On the efficiency and consitency of covariance localisation in the - - PowerPoint PPT Presentation

On the efficiency and consitency of covariance localisation in the EnKF

Alban Farchi and Marc Bocquet

CEREA, laboratoire commun École des Ponts ParisTech and EDF R&D, Université Paris-Est, Champs-sur-Marne, France

Tuesday, 4 June 2018

International EnKF workshop

Alban Farchi Localisation of the EnSRF 4 June 2019 1 / 26

Contents

1

Rank deficiency of the EnKF

2

The LEnSRF with augmented ensembles

3

A new perturbation update scheme

4

References

Alban Farchi Localisation of the EnSRF 4 June 2019 2 / 26

Rank deficiency of the EnKF The Kalman filter

The ensemble Kalman filter analysis

◮ We focus on the Kalman filter analysis step:

K = BHT R + HBHT−1 , (1) xa = xb + K (y − Hxb) , (2) Pa = (I − KH) B, (3)

◮ In the EnKF, the statistics are carried through by the ensemble xi, i = 1 . . . Ne

• :

xb = x, (4) B = XXT, (5)

where x is the ensemble mean and X is the normalised anomaly matrix.

Alban Farchi Localisation of the EnSRF 4 June 2019 3 / 26

Rank deficiency of the EnKF Rank deficiency and localisation

Rank deficiency of the EnKF

◮ The matrix XXT has rank limited by Ne − 1, too small to accurately represent B in a high-dimensional system (Nx ≫ Ne). ◮ When Ne is too small, XXT is characterised by large sampling errors, which often take the form of spurious correlations at long distance. ◮ A common solution is to use localisation: either make the analysis local (domain localisation) or use a localised B (covariance localisation).

domain localisation covariance localisation relies on a collection of local analyses relies on a localised background error B embarrassingly parallel no obvious parallelisation of the perturbation update ad hoc treatment of non-local observations rigourous treatment of non-local observations Alban Farchi Localisation of the EnSRF 4 June 2019 4 / 26

Rank deficiency of the EnKF Rank deficiency and localisation

Covariance localisation in the deterministic EnKF

The (exact) LEnSRF Regularise XXT with a localisation matrix ρ:

B = ρ ◦ XXT . (6)

Update the perturbation as:

T = I + BHTR−1H, (7) Xa = T−1/2X. (8)

◮ First focus: efficient implementation (accuracy/speed) of the perturbation update. ◮ Second focus: consistency (how well does XaXT

a match Pa) of the perturbation

update.

Alban Farchi Localisation of the EnSRF 4 June 2019 5 / 26

The LEnSRF with augmented ensembles

Contents

1

Rank deficiency of the EnKF

2

The LEnSRF with augmented ensembles The augmented ensemble Update the perturbations Numerical illustration: the Lorenz 1996 model Numerical illustration: the multi-layer Lorenz 1996 model

3

A new perturbation update scheme

4

References

Alban Farchi Localisation of the EnSRF 4 June 2019 6 / 26

The LEnSRF with augmented ensembles The augmented ensemble

Using an augmented ensemble to represent B

◮ In the analysis step, we choose to compute Ne > Ne perturbations X such that B = ρ ◦ XXT ≈ X

• XT. There are two methods.

The modulation method Suppose that there is a matrix W with Nm columns such that ρ ≈ WWT. Let X be the matrix with NmNe columns:

• XjNe+i

n

= [W]j

n [X]i n .

(9) The random svd method Compute the truncated svd B = ρ ◦ XXT ≈ UΣUT with Nm columns, and form the matrix

• X = UΣ1/2.

Alban Farchi Localisation of the EnSRF 4 June 2019 7 / 26

The LEnSRF with augmented ensembles The augmented ensemble

The modulation method

1:

Compute the modulation product X = W∆X

• XjNe+i

n

= [W]j

n [X]i n .

(10)

◮ This is a mix between a Schur product (for the state variable index n) and a tensor product (for the ensemble indices i and j). [Buehner, 2005] ◮ The modulation product is based on a factorisation property shown by [Lorenc, 2003] and is currently used for covariance localisation [e.g., Bishop et al., 2017], including in

• perational centers [e.g., Arbogast et al., 2017].

Alban Farchi Localisation of the EnSRF 4 June 2019 8 / 26

The LEnSRF with augmented ensembles The augmented ensemble

The random svd method

1:

Compute the truncated svd B = ρ ◦ XXT ≈ UΣUT with Nm columns.

2:

Form the matrix X = UΣ1/2.

◮ The matrix B is sparse, which means that efficient (and parallelisable) methods can be used to compute the truncated svd, e.g., the random svd algorithm of [Halko et al., 2011].

[Farchi and Bocquet, 2019]

Alban Farchi Localisation of the EnSRF 4 June 2019 9 / 26

The LEnSRF with augmented ensembles Update the perturbations

Compute the updated perturbations

◮ How to obtain Ne updated members using the Ne analysis perturbations Xa of the augmented ensemble? ◮ A solution is to use the augmented ensemble to compute an approximate left-transform update of the (non-augmented) ensemble as follows:

Xa = I + BHTR−1H−1/2 X, (11) Xa = I + X XTHTR−1H−1/2 X, (12) Xa =

• I −

X

• I +

YTR−1 Y + I + YTR−1 Y1/2−1 YTR−1 H

• X,

(13)

where Y = H X. ◮ The linear algebra is performed in the augmented ensemble space (i.e., using Ne × Ne matrices) using the formula derived by [Bocquet, 2016], later used by [Bishop et al., 2017] under the name gain form of the ETKF.

Alban Farchi Localisation of the EnSRF 4 June 2019 10 / 26

The LEnSRF with augmented ensembles Numerical illustration: the Lorenz 1996 model

The Lorenz 1996 model

◮ We use the L96 model with Nx = 400 variables:

dxn dt = (xn+1 − xn−2) xn−1 − xn + 8, n = 1 . . . Nx. (14)

◮ The observations are given every ∆t = 0.05 by

y = x + v, v ∼ N (0, I) . (15)

◮ The localisation matrix is constructed using the Gaspari–Cohn function, assuming that the grid points are equally distributed in space. ◮ All algorithms use an ensemble of Ne = 10 members. ◮ The runs are 2 × 104∆t long and our criterion is the time-average analysis RMSE.

Alban Farchi Localisation of the EnSRF 4 June 2019 11 / 26

The LEnSRF with augmented ensembles Numerical illustration: the Lorenz 1996 model

Results with the L96 model

50 100 150 200 250 300 Augmented ensemble size Ne 0.20 0.25 0.30 0.35 0.40 0.45 RMSE LETKF random svd modulation 0.20 0.22 0.24 0.26 0.28 0.30 RMSE 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 Wall-clock time

• ×103s
• random svd, 1 thread

random svd, 24 threads modulation

◮ Both methods can yield similar RMSE scores as the LETKF. ◮ The modulation method requires a larger augmented ensemble size to yield similar RMSE scores as the random svd method. ◮ For a given level of RMSE score, the random svd method is faster than the modulation method.

[Farchi and Bocquet, 2019]

Alban Farchi Localisation of the EnSRF 4 June 2019 12 / 26

The LEnSRF with augmented ensembles Numerical illustration: the multi-layer Lorenz 1996 model

A multilayer extension of the L96 model

◮ We introduce the mL96 model, that consists of Pz = 32 coupled layers of the L96 model with Ph = 40 variables:

dx(z, h) dt = x(z, h+1) − x(z, h−2)

• x(z, h−1) − x(z, h) + Fz

+ δ{z>0}

• x(z−1, h) − x(z, h)
• + δ{z<Pz}
• x(z+1, h) − x(z, h)
• .

(16)

◮ The forcing term linearly decreases from F1 = 8 to F32 = 4.

Alban Farchi Localisation of the EnSRF 4 June 2019 13 / 26

The LEnSRF with augmented ensembles Numerical illustration: the multi-layer Lorenz 1996 model

Satellite observations for the mL96 model

◮ Each column is observed independently via: yc, h =

Pz

• z=1

[Ω]c, z xz, h + vc, h, vc, h ∼ N (0, 1) , (17) where Ω is a weighting matrix with Nc = 8 channels that is designed to mimic satellite radiances. ◮ The 8 × 40 observations are available every ∆t = 0.05. ◮ The runs are 104∆t long. ◮ All algorithms use an ensemble of Ne = 8 members.

0.0 0.1 0.2 0.3 0.4 Weight 1 8 16 24 32 Vertical layer index z

Covariance localisation (with augmented ensembles) is used only in the vertical

• direction. Domain localisation (LETKF-like) is used in the horizontal direction.

Alban Farchi Localisation of the EnSRF 4 June 2019 14 / 26

The LEnSRF with augmented ensembles Numerical illustration: the multi-layer Lorenz 1996 model

Results with the mL96 model

◮ Using covariance localisation in the vertical direction yields better RMSE scores than the LETKF. ◮ The modulation method requires a larger augmented ensemble size to yield similar RMSE scores as the random svd method. ◮ Both methods benefit from the parallelisation of the local analyses, but the parallelisation potential of the random svd method is not fully exploited because of our limited computational platform.

[Farchi and Bocquet, 2019]

1 5 9 13 17 21 0.1 0.2 0.3 0.4 0.5 RMSE ad hoc LETKF

• trunc. svd,

Ne = 8

• trunc. svd,

Ne = 32

• trunc. svd,

Ne = 64 modulation, Ne = 8 modulation, Ne = 32 modulation, Ne = 128 1 5 9 13 17 21 Horizontal localisation radius rh 101 102 103 Wall-clock time [s]

Alban Farchi Localisation of the EnSRF 4 June 2019 15 / 26

The LEnSRF with augmented ensembles Numerical illustration: the multi-layer Lorenz 1996 model

Conclusions related to the augmented ensembles

◮ We have shown how to use augmented ensembles to implement covariance localisation in the deterministic EnKF. ◮ We have proposed an alternative to the modulation method to construct the augmented ensemble, based on the random svd algorithm. ◮ We have compared both methods using the L96 model and a multilayer extension of the L96 model with satellite-like observations. ◮ We have shown that the augmented ensemble size need to be smaller with the random svd method than with the modulation method. ◮ We have seen that using domain localisation in the horizontal and covariance localisation in the vertical seems to be an adequate approach to assimilate satellite radiances in a spacially extended model.

Alban Farchi Localisation of the EnSRF 4 June 2019 16 / 26

A new perturbation update scheme

Contents

1

Rank deficiency of the EnKF

2

The LEnSRF with augmented ensembles

3

A new perturbation update scheme Improving the consistency of the perturbation update Numerical cost of computing the gradient and the cost function Numerical illustrations: the Lorenz 1996 and the Kuramoto–Sivashinsky models

4

References

Alban Farchi Localisation of the EnSRF 4 June 2019 17 / 26

A new perturbation update scheme Improving the consistency of the perturbation update

Failed deterministic sampling

◮ The perturbation update of the LEnSRF:

T = I + BHTR−1H, (18) Xa = T−1/2X. (19)

◮ A first source of inconsistency lies in the fact that B = XXT. ◮ A potential fix: apply the left-transform to X, defined as the Ne dominant modes of B = ρ ◦ XXT , for which B ≈ X XT. ◮ Empirically, this fix systematically makes the EnKF diverge. Intuitively, this comes from a double application of localisation.

Alban Farchi Localisation of the EnSRF 4 June 2019 18 / 26

A new perturbation update scheme Improving the consistency of the perturbation update

A new approach

◮ For the LEnSRF update, instead of using the left-transform it could be more consistent with how the perturbations are defined to look for a low-rank perturbation matrix Xa such that

Pa ≈ ρ ◦ XaXT

a

• .

(20)

A solution of Eq. (20) trades the accuracy of the representation of the long range covariances (which may eventually be discarded at the next cycle) for a potentially better accuracy of the short range covariances. ◮ In this context, the goal is to solve the optimisation problem Xa = arg min

rank(X)≤Ne−1

L (X) , with L (X) = ln

• ρ ◦

XXT − Pa

• F ,

(21) where ∗F is the Frobenius matrix norm.

[Bocquet and Farchi, 2019]

Alban Farchi Localisation of the EnSRF 4 June 2019 19 / 26

A new perturbation update scheme Improving the consistency of the perturbation update

A new approach

◮ This problem is similar to a weighted low-rank approximation (WLRA) problem. We expect that it has no tractable solution as opposed to the unweighted case for which the Eckart-Young solution holds. ◮ However, the gradient of the cost function can be computed: ∇L (X) = 2 ∆−2

F (ρ ◦ ∆) X,

(22) with ∆ = ρ ◦ XXT − Pa. ◮ Therefore, a local solution can be found using a gradient-based minimisation algorithm.

Alban Farchi Localisation of the EnSRF 4 June 2019 20 / 26

A new perturbation update scheme Numerical cost of computing the gradient and the cost function

Computing the gradient and the cost function

◮ We use a mode expansion with Ne modes for the analysis error covariance matrix: Pa ≈ Xa XT

a .

◮ To compute L and ∇L, we use the (classical) formula ρ ◦ XXT · v =

Ne

• i=1

X(i) ◦ ρ · X(i) ◦ v , (23) for any matrix X of size Nx × Ne and any vector v of size Nx, and where X(i) is the i-th column of X. ◮ Therefore, the numerical cost of computing L and ∇L is: O

• Ne

NeNxNb

• if ρ is banded with bandwidth Nb;

O

• Ne

NeNxlnNx

• if ρ is homogeneous.

◮ The numerical cost can be reduced by a factor Ne or Ne if Eq. (23) is computed in parallel.

Alban Farchi Localisation of the EnSRF 4 June 2019 21 / 26

A new perturbation update scheme Numerical illustrations: the L96 and the KS models

The Lorenz 1996 and the Kuramoto–Sivashinsky models

◮ We use the L96 model with Nx = 40 variables, in the same configuration as in the previous test series. ◮ As a complement, we use the KS model:

∂tu = −u∂xu − ∂2

xu − ∂4 xu.

(24)

Equation (24) is solved over the domain x ∈ [0, 32π] using a pseudo-spectral integration with Nx = 128. For this model, the observations are given every ∆t = 1 by

y = x + v, v ∼ N (0, I) . (25)

◮ For both models, the localisation matrix is constructed using the Gaspari–Cohn function, assuming that the grid points are equally distributed in space. ◮ The runs are 2 × 104∆t long and our criterion is the time-average analysis RMSE. ◮ We compare: the LETKF algorithm; the LEnSRF algorithm (exact left transform update); the new LEnSRF algorithm (new perturbation update scheme).

Alban Farchi Localisation of the EnSRF 4 June 2019 22 / 26

A new perturbation update scheme Numerical illustrations: the L96 and the KS models

Accuracy for the L96 model and KS models

4 6 8 10 12 14 16 18 20

Ensemble size Ne

0.18 0.20 0.22 0.24 0.26 0.28 0.30

Average RMSE L96 LEnSRF LETKF new LEnSRF

4 6 8 10 12 14 16 18 20

Ensemble size Ne

1.00 1.02 1.04 1.06 1.08 1.10

Optimal multiplicative inflation L96 LEnSRF LETKF new LEnSRF

4 6 8 10 12 14 16 18 20

Ensemble size Ne

5 10 15 20 25 30

Optimal localisation length L96 LEnSRF LETKF new LEnSRF

4 6 8 10 12 14 16

Ensemble size Ne

0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18

Average RMSE KS LEnSRF LETKF new LEnSRF

4 6 8 10 12 14 16

Ensemble size Ne

1.00 1.02 1.04 1.06 1.08 1.10 1.12

Optimal multiplicative inflation KS LEnSRF LETKF new LEnSRF

4 6 8 10 12 14 16

Ensemble size Ne

10 15 20 25 30 35 40 45 50 55

Optimal localisation length KS LEnSRF LETKF new LEnSRF

Comparison of the LETKF, the LEnSRF, and the LEnSRF with the new update scheme, applied to L96 (top) and to KS (bottom). The RMSE, optimal localisation and optimal inflation are plotted as functions of the ensemble size Ne. [Bocquet and Farchi, 2019] Alban Farchi Localisation of the EnSRF 4 June 2019 23 / 26

A new perturbation update scheme Numerical illustrations: the L96 and the KS models

Robustness of the new scheme

1.0 1.005 1.01 1.015 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 1.1

Multiplicative inflation

0.30 0.35 0.40 0.45 0.50 0.55

Average RMSE L96 Ne =4 LEnSRF LETKF new LEnSRF

1.0 1.005 1.01 1.015 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 1.1

Multiplicative inflation

0.20 0.25 0.30 0.35 0.40 0.45

Average RMSE L96 Ne =8 LEnSRF LETKF new LEnSRF

1.0 1.005 1.01 1.015 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 1.1

Multiplicative inflation

0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25 0.26

Average RMSE L96 Ne =16 LEnSRF LETKF new LEnSRF

1 . 1 . 5 1 . 1 1 . 1 5 1 . 2 1 . 4 1 . 6 1 . 8 1 . 1 1 . 1 2 1 . 1 4 1 . 1 6

Multiplicative inflation

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Average RMSE KS Ne =4 LEnSRF LETKF new LEnSRF

1.0 1.005 1.01 1.015 1.02 1.04 1.06 1.08 1.1 1.12 1.14 1.16

Multiplicative inflation

0.15 0.20 0.25 0.30 0.35

Average RMSE KS Ne =8 LEnSRF LETKF new LEnSRF

1.0 1.005 1.01 1.015 1.02 1.04 1.06 1.08 1.1 1.12 1.14 1.16

Multiplicative inflation

0.12 0.14 0.16 0.18 0.20

Average RMSE KS Ne =16 LEnSRF LETKF new LEnSRF

Time-averaged RMSE as a function of the multiplicative inflation, the localisation length being tuned so as to minimise the RMSE for L96 (top) and KS (bottom), for Ne = 4, 8, 16. [Bocquet and Farchi, 2019] Alban Farchi Localisation of the EnSRF 4 June 2019 24 / 26

A new perturbation update scheme Numerical illustrations: the L96 and the KS models

Conclusions related to the new perturbation update scheme

◮ The updated perturbations in the local EnKFs based on covariance localisation (in particular the LEnSRF) are not the main modes of Pa (= LETKF). ◮ We have proposed a perturbation update scheme potentially more consistent such that the perturbations X are related to the error covariance matrix by P ≈ ρ ◦ XXT throughout the EnKF: Xa = arg min

rank(X)≤Ne−1

L (X) , with L (X) = ln

• ρ ◦

XXT − Pa

• F ,

(26) ◮ We have compared numerically the new LEnSRF to the LETKF and to the LEnSRF (with the exact left-transform update), using the L96 and the KS models. ◮ We have shown that for both models, the requirement for residual multiplicative inflation is much weaker with the new LEnSRF than with both the LETKF and the LEnSRF: much weaker imbalance of the new update scheme? ◮ Moreover, there is an accuracy improvement of up to 6% in the analysis RMSE in mildly nonlinear conditions, which is significant in these very well tuned configurations.

Alban Farchi Localisation of the EnSRF 4 June 2019 25 / 26

References

References

Bocquet, M. and A. Farchi (2019). “On the consistency of the local ensemble square root Kalman filter perturbation update”. In: Tellus A. In press. Farchi, A. and M. Bocquet (2019). “On the Efficiency of Covariance Localisation of the Ensemble Kalman Filter Using Augmented Ensembles”. In: Front. App. M. Stat. 5, p. 3. Alban Farchi Localisation of the EnSRF 4 June 2019 26 / 26