World-class research. Ultimate impact.
More on impact ›

General Commentary ARTICLE

Front. Appl. Math. Stat., 17 March 2020 | https://doi.org/10.3389/fams.2020.00002

Commentary: On the Efficiency of Covariance Localisation of the Ensemble Kalman Filter Using Augmented Ensembles

  • 1School of Earth Sciences and Centre of Excellence for Climate Extremes, University of Melbourne, Parkville, VIC, Australia
  • 2NOAA Earth System Research Lab, Boulder, CO, United States
  • 3School of Atmospheric Sciences, Nanjing University, Nanjing, China

A Commentary on
On the Efficiency of Covariance Localisation of the Ensemble Kalman Filter Using Augmented Ensembles

by Farchi, A., and Bocquet, M. (2019). On the Efficiency of Covariance Localisation of the Ensemble Kalman Filter Using Augmented Ensembles. Front. Appl. Math. Stat. 5:3. doi: 10.3389/fams.2019.00003

Introduction

In discussing Equation (39) (Equation (25) of Bocquet [1]), Farchi and Bouquet [2] state that “This perturbation update has been rediscovered by Bishop et al. [3] and included in their gain ETKF (GETKF) algorithm. However, the update formula used in the GETKF is prone to numerical cancellation errors as opposed to Equation (39)”. Here, we note:

  (i) The predecessor of the GETKF eigenvalue form of the modified gain matrix equation appeared in Posselt and Bishop [4, 5]—before Bocquet [1].

 (ii) The spectral shift theorem reduces the differences in the numerical cancellation errors referred to by Farchi and Bouquet.

(iii) The eigenvalue form enables Wang et al.'s [6] corrections for ensemble rank deficiency.

(iv) A proof of the equivalence of the eigenvalue form and Bouquet's form.

On page 12, Farchi and Bouquet [2] also state that “Such an extension had been discussed by Bishop et al. [3] but without numerical illustration.” This is incorrect. Lei et al. [7] used the GETKF to show that model space ensemble covariance localization provided satellite data assimilation (DA) performance comparable to 3DEnsVar.

To be specific about the forms of the modified gain matrix, let

  Zf=[((x1f)-(xf)¯),((x2f)-(xf)¯),,((xKf)-(xf)¯)]K-1, andH~Zf=R-1/2[(H(x1f)-H(xf)¯),(H(x2f)-H(xf)¯),,(H(xKf)-H(xf)¯)]K-1    (1)

where K is the total number of ensemble members in the ensemble forecast and where the n-vector xif is the ith member of the prior ensemble forecast and where the p-vector H(xif) is the ith member of the prior ensemble forecast of the p-vector y of p observations. When p<K, the numerical cost of the pxp eigen decomposition

H~Zf(H~Zf)T=EΓpxpET    (2)

is less than the K × K eigen decomposition

(H~Zf)TH~Zf=CΓK×KCT.    (3)

In (2), E is a pxp eigenvector matrix for which EET=ETE=Ipxp, Γpxp is a pxp diagonal matrix of eigenvalues. In (3), C is a K × K orthonormal matrix of eigenvectors (CCT=CTC=IK×K) and ΓKxK is a K × K diagonal matrix of eigenvalues. At least K-p of the eigenvalues in ΓK × K will be equal to zero in the case of K>p. Equation's (2) and (3) are directly connected to the verbose singular value decomposition H~Zf=EpxpΓpxK1/2CK×KT where ΓpxK1/2=[Γpxp1/2 0px(K-p)] where 0px(Kp) is a px(Kp) matrix of zeros. However, since the columns of C associated with zero eigenvalues cannot contribute to products of the matrix H~Zf with other vectors, it is more efficient to work with the concise svd given by H~Zf=EpxpΓpxp1/2(LKxp)T where LKxp lists the p columns of CK×K having non-zero eigenvalues. Posselt and Bishop [4, 5] note that LKxp is given by

(LKxp)T=Γpxp-1/2ETH~Zf or LKxp=(H~Zf)TEΓpxp-1/2    (4)

and hence can be computed without performing an eigen decomposition of the larger K × K matrix in (3). Posselt and Bishop [4, 5] prove that for a linear observation operator H~, if

Za={I-ZfLKxp[Ipxp-(Γpxp+I)-1/2]Γpxp-1/2ETH~}Zf    (5)

(see [5], Equation A10) then

Pa=ZaZaT     =ZfZfT-Zf(H~Zf)T((H~Zf)(H~Zf)T+Ipxp)-1(H~Zf)ZfT     =Pf-PfH~T(H~PfH~T+Ipxp)-1H~Pf.    (6)

The analysis perturbations are given by Xa=ZaK-1, hence,

Xa={IKxK-ZfLKxp[Ipxp-(Γpxp+I)-1/2]Γpxp-1/2ETH~}Xf,where Xf=ZfK-1    (7)

is the perturbation update equation implied by Posselt and Bishop [4, 5].

In the above notation, and when propagation of small amplitude ensemble perturbations by the non-linear model is replaced by the propagation of raw ensemble perturbations by the non-linear model (i.e., no tangent linear model approximation is made), Bocquet's Equation (25) [1] for the ensemble perturbation update takes the form,

XBouqueta={IK×KZf[IK×K+(H˜Zf)T(H˜Zf)+(IK×K+(H˜Zf)T(H˜Zf))1/2]1(H˜Zf)TH˜}Xf.    (8)

A fundamental difference between (7) and (8) is that while Bouquet multiplies Zf by the KxK matrix [IK×K+(H~Zf)T(H~Zf)+(IK×K+(H~Zf)T(H~Zf))1/2]-1 Posselt and Bishop multiply it by the Kxp matrix LKxp[Ipxp-(Γpxp+I)-1/2]. When K>p, Posselt and Bishop's form only requires the eigenvector decomposition of a pxp matrix, whereas Bouquet's form requires the inversion of a larger K × K matrix. However, when p>K, the eigen decomposition (3) is cheaper than (2), LKxp becomes identical to the K × K matrix CK×K and H~Zf=EpxKΓK×K1/2CK×KT becomes the concise svd of H~Zf. In this case, EpxK is efficiently given by H~ZfCK×KΓK×K-1/2=EpxK and (7) becomes

Xa={IZfCK×K[IK×K(ΓK×K+I)1/2]ΓK×K1CK×KT(H˜Zf)TH˜}Xf    (9)

Dividing Equation (9) by K-1 recovers Equation (24) of Bishop et al. [3].

The above shows that Bishop et al.'s Equation (24) [3] was not “rediscovered” from Bocquet's [1] form as implied by Farchi and Bocquet [2]. It is an extension of Posselt and Bishop's [4, 5] eigenvalue form to the case of K>p. Equation (9) is just an eigenvalue form of the modified gain matrix of Whitaker and Hamill's [8] Ensemble Square Root Filter.

Equivalence of (9) and (8)

Bocquet's Equation (25) [1] can be derived from (9) with the following steps:

  (i) Drop the dimension subscripts and manipulate [I − (Γ + I)−1/2−1 as follows

Γ-1[I-(Γ+I)-1/2]=Γ-1[I-(Γ+I)-1/2][I+(Γ+I)-1/2][I+(Γ+I)-1/2]-1=Γ-1[I-(Γ+I)-1][I+(Γ+I)-1/2]-1=Γ-1[(Γ+I)(Γ+I)-1-(Γ+I)-1][I+(Γ+I)-1/2]-1=Γ-1[Γ(Γ+I)-1][I+(Γ+I)-1/2]-1=Γ-1Γ[(Γ+I)[I+(Γ+I)-1/2]]-1=[(Γ+I)+(Γ+I)1/2]-1    (10)

 (ii) Use (10) in (9) to give

    Xa={I-ZfC[(Γ+I)+(Γ+I)1/2]-1CT(H~Zf)TH~}Xf={I-Zf[C(Γ+I)CT+C(Γ+I)1/2CT]-1(H~Zf)TH~}Xf={I-Zf[(CΓCT+I)+C(Γ+I)1/2CT]-1(H~Zf)TH~}Xf    (11)

(iii) But

[(CΓCT+I)]1/2=[C(Γ+I)CT]1/2=C(Γ+I)1/2CT    (12)

(iv) Using (12) and (3) in (11) gives

Xa={IZf[(CΓCT+I)+(CΓCT+I)1/2]1(H˜Zf)TH˜}Xf={IZf[((H˜Zf)TH˜Zf+I)+((H˜Zf)TH˜Zf+I)1/2]1(H˜Zf)TH˜}Xf    (13)

Equation (13) is equivalent to (8) and Bocquet's Equation (25) [1].

Numerical Issues, Condition Numbers, and Understanding

Numerical cancellation errors increase when the condition number of the matrix increases. Let us define the scalars γimax and γimin to, respectively, denote the maximum and minimum of the eigenvalues listed in the eigenvalue matrix Γpxp. The condition number of (H~Z)TH~Z is κ[(H~Z)TH~Z]=γimaxγimin. Because γimin can be zero, κ[(H~Z)TH~Z] can be infinite. In contrast, κ[((H~Z)TH~Z+I)+((H~Z)TH~Z+I)1/2]=(γimax+1)+(γimax+1)1/2(γimin+1)+(γimin+1)1/2 is bounded above by [(γimax+1)+(γimax+1)1/2]/2. However, note that the matrix [(H~Z)TH~Z+αI] has the eigenvalue decomposition

[(H~Z)TH~Z+αI]=CΛCT=C(Γ+αI)CT    (14)

and hence has κ[(H~Z)TH~Z+αI]=γimax+αγimin+α which is bounded above by γimaxα+1 where α is a positive scalar. Hence, α can be chosen to create a matrix that is better conditioned than [((H~Z)TH~Z+I)+((H~Z)TH~Z+I)1/2]. Once the eigen decomposition CΛCT of (14) has been obtained, one obtains the eigenvalues required by the GETKF or ETKF using Γ = Λ − αI. Thus, condition number differences between the Bouquet and eigenvalue form are easily eliminated.

The eigenvalue form lends understanding to the performance of DA schemes in much the same way that Empirical Orthogonal Functions lend understanding to climate variability. Wang et al. [6] used this understanding to correct gross aspects of the eigenvalue overestimation that occurs when the size of the ensemble is much smaller than the rank of the true observation space forecast error covariance matrix.

Discussion

Bocquet [1] and Farchi and Bocquet [2] may have overlooked Posselt and Bishop's [4, 5] work because:

(i)   It is difficult to find all relevant literature to one's own work.

(ii)  Bishop et al. [3] did not cite Posselt and Bishop [4, 5].

(iii) The equivalence of Posselt and Bishop's [4, 5] form and Bocquet's [1] form is not obvious.

Similarly, Bishop et al. [3] overlooked Bocquet's [1] work because of (i) and (iii). This note serves to clarify the origins and uses of modified gain matrices used in ensemble DA.

Author Contributions

CB led the study and wrote the text. JW and LL carefully reviewed the text.

Funding

CB acknowledges support from the Australian Research Council's Centre of Excellence in Climate Extremes (CE170100023). LL acknowledges joint sponsorship by the National Key R&D Program of China through grant 2017YFC1501603 and the National Natural Science Foundation of China through grant 41675052. JW acknowledges support through the NOAA High-Impact Weather Prediction Project (HIWPP) under award NA14OAR4830123, and the NOAA/NWS Next-Generation Global Prediction System (NGGPS) project.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

1. Bocquet M. Localization and the iterative ensemble Kalman smoother. Q J R Meteorol Soc. (2016) 142:1075–89. doi: 10.1002/qj.2711

CrossRef Full Text | Google Scholar

2. Farchi A, Bocquet M. On the efficiency of covariance localisation of the ensemble Kalman filter using augmented ensembles. Front Appl Math Stat. (2019) 5:3. doi: 10.3389/fams.2019.00003

CrossRef Full Text | Google Scholar

3. Bishop CH, Whitaker JS, Lei L. Gain form of the ensemble transform Kalman filter and its relevance to satellite data assimilation with model space ensemble covariance localization. Mon Wea Rev. (2017) 145:4575–92. doi: 10.1175/MWR-D-17-0102.1

CrossRef Full Text | Google Scholar

4. Posselt DJ, Bishop CH. Nonlinear parameter estimation: comparison of an ensemble Kalman smoother with a Markov chain Monte Carlo algorithm. Mon Wea Rev. (2012) 140:1957–74. doi: 10.1175/MWR-D-11-00242.1

CrossRef Full Text | Google Scholar

5. Posselt DJ. Bishop CH. Corrigendum. Mon Wea Rev. (2014) 142:1382. doi: 10.1175/MWR-D-13-00342.1

CrossRef Full Text

6. Wang X, Hamill TM, Whitaker JS, Bishop CH. A comparison of hybrid ensemble transform Kalman filter–optimum interpolation and ensemble square root filter analysis schemes. Mon Wea Rev. (2007) 135:1055–76. doi: 10.1175/MWR3307.1

CrossRef Full Text | Google Scholar

7. Lei L, Whitaker JS, Bishop C. Improving assimilation of radiance observations by implementing model space localization in an ensemble Kalman filter. J Adv Model Earth Syst. (2018) 10:3221–32. doi: 10.1029/2018MS001468

CrossRef Full Text | Google Scholar

8. Whitaker JS, Hamill TM. Ensemble data assimilation without perturbed observations. Mon Wea Rev. (2002) 130:1913–24. doi: 10.1175/1520-0493(2002)130<1913:EDAWPO>2.0.CO;2

CrossRef Full Text | Google Scholar

Keywords: ensemble Kalman filter, modified Kalman gain, eigenvalue form of Kalman gain, GETKF, ETKF

Citation: Bishop CH, Whitaker JS and Lei L (2020) Commentary: On the Efficiency of Covariance Localisation of the Ensemble Kalman Filter Using Augmented Ensembles. Front. Appl. Math. Stat. 6:2. doi: 10.3389/fams.2020.00002

Received: 26 November 2019; Accepted: 29 January 2020;
Published: 17 March 2020.

Edited by:

Raluca Eftimie, University of Dundee, United Kingdom

Reviewed by:

Xin Tong, National University of Singapore, Singapore

Copyright © 2020 Bishop, Whitaker and Lei. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Craig H. Bishop, craig.bishop@unimelb.edu.au