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

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

$\begin{array}{l}{\text{Z}}^f=\frac{\left[\left(\left({\text{x}}_1^f\right)-\overline{\left({\text{x}}_{}^f\right)}\right),\left(\left({\text{x}}_2^f\right)-\overline{\left({\text{x}}_{}^f\right)}\right),\dots ,\left(\left({\text{x}}_K^f\right)-\overline{\left({\text{x}}_{}^f\right)}\right)\right]}{\sqrt{K-1}},\text{ and}\\ \tilde{H}{\text{Z}}^f=\frac{{\text{R}}^{-1/2}\left[\left(H\left({\text{x}}_1^f\right)-\overline{H\left({\text{x}}_{}^f\right)}\right),\left(H\left({\text{x}}_2^f\right)-\overline{H\left({\text{x}}_{}^f\right)}\right),\dots ,\left(H\left({\text{x}}_K^f\right)-\overline{H\left({\text{x}}_{}^f\right)}\right)\right]}{\sqrt{K-1}}& (1)\end{array}$

where K is the total number of ensemble members in the ensemble forecast and where the n-vector ${\text{x}}_i^f$ is the ith member of the prior ensemble forecast and where the p-vector $H\left({\text{x}}_i^f\right)$ 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

$\begin{array}{ll}\tilde{H}{\text{Z}}^f{\left(\tilde{H}{\text{Z}}^f\right)}^T=\text{E}{\Gamma }_{pxp}{\text{E}}^T& (2)\end{array}$

is less than the K × K eigen decomposition

$\begin{array}{ll}{\left(\tilde{H}{\text{Z}}^f\right)}^T\tilde{H}{\text{Z}}^f=\text{C}{\Gamma }_{K\times K}{\text{C}}^T.& (3)\end{array}$

In (2), E is a pxp eigenvector matrix for which $\text{E}{\text{E}}^T={\text{E}}^T\text{E}={\text{I}}_{pxp}$, Γ_pxp is a pxp diagonal matrix of eigenvalues. In (3), C is a K × K orthonormal matrix of eigenvectors $(\text{C}{\text{C}}^T={\text{C}}^T\text{C}={\text{I}}_{K\times 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 $\tilde{H}{\text{Z}}^f={\text{E}}_{pxp}{\Gamma }_{pxK}^{1/2}{\text{C}}_{K\times K}^T$ where ${\Gamma }_{pxK}^{1/2}=\left[{\Gamma }_{pxp}^{1/2}{\text{0}}_{px(K-p)}\right]$ where 0_px(K−p) is a px(K−p) matrix of zeros. However, since the columns of C associated with zero eigenvalues cannot contribute to products of the matrix $\tilde{H}{\text{Z}}^f$ with other vectors, it is more efficient to work with the concise svd given by $\tilde{H}{\text{Z}}^f={\text{E}}_{pxp}{\Gamma }_{pxp}^{1/2}{\left({\text{L}}_{Kxp}^{}\right)}^T$ where ${\text{L}}_{Kxp}^{}$ lists the p columns of ${\text{C}}_{K\times K}^{}$ having non-zero eigenvalues. Posselt and Bishop [4, 5] note that ${\text{L}}_{Kxp}^{}$ is given by

$\begin{array}{ll}{\left({\text{L}}_{Kxp}^{}\right)}^T={\Gamma }_{pxp}^{-1/2}{\text{E}}^T\tilde{H}{\text{Z}}^f\text{ or }{\text{L}}_{Kxp}^{}={\left(\tilde{H}{\text{Z}}^f\right)}^T\text{E}{\Gamma }_{pxp}^{-1/2}& (4)\end{array}$

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 $\tilde{H}$, if

$\begin{array}{ll}{\text{Z}}^a=\left\{\text{I}-{\text{Z}}^f{\text{L}}_{Kxp}^{}\left[{\text{I}}_{pxp}-{({\Gamma }_{pxp}+\text{I})}^{-1/2}\right]{\Gamma }_{pxp}^{-1/2}{\text{E}}^T\tilde{H}\right\}{\text{Z}}^f& (5)\end{array}$

(see [5], Equation A10) then

$\begin{array}{l}{\text{P}}^a={\text{Z}}^a{\text{Z}}^{aT}\\ ={\text{Z}}^f{\text{Z}}^{fT}-{\text{Z}}^f{\left(\tilde{H}{\text{Z}}^f\right)}^T{\left(\left(\tilde{H}{\text{Z}}^f\right){\left(\tilde{H}{\text{Z}}^f\right)}^T+{\text{I}}_{pxp}\right)}^{-1}\left(\tilde{H}{\text{Z}}^f\right){\text{Z}}^{fT}\\ ={\text{P}}^f-{\text{P}}^f{\tilde{H}}^T{\left(\tilde{H}{\text{P}}^f{\tilde{H}}^T+{\text{I}}_{pxp}\right)}^{-1}\tilde{H}{\text{P}}^f.& (6)\end{array}$

The analysis perturbations are given by ${\text{X}}^a={\text{Z}}^a\sqrt{K-1}$, hence,

$\begin{array}{l}{\text{X}}^a=\left\{{\text{I}}_{KxK}-{\text{Z}}^f{\text{L}}_{Kxp}\left[{\text{I}}_{pxp}-{({\Gamma }_{pxp}+\text{I})}^{-1/2}\right]{\Gamma }_{pxp}^{-1/2}{\text{E}}^T\tilde{H}\right\}{\text{X}}^f,\\ \text{where }{\text{X}}^f={\text{Z}}^f\sqrt{K-1}& (7)\end{array}$

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,

$\begin{array}{l}{X}_{Bouquet}^a=\{I_{K\times K}-Z^f[I_{K\times K}+{(\tilde{H}{Z}^f)}^T(\tilde{H}{Z}^f)\\ +{(I_{K\times K}+{(\tilde{H}{Z}^f)}^T(\tilde{H}{Z}^f))}^{1/2}{]}^{-1}{(\tilde{H}{Z}^f)}^T\tilde{H}\}{X}^f.& (8)\end{array}$

A fundamental difference between (7) and (8) is that while Bouquet multiplies Z^f by the KxK matrix ${\left[{\text{I}}_{K\times K}+{\left(\tilde{H}{\text{Z}}^f\right)}^T\left(\tilde{H}{\text{Z}}^f\right)+{\left({\text{I}}_{K\times K}+{\left(\tilde{H}{\text{Z}}^f\right)}^T\left(\tilde{H}{\text{Z}}^f\right)\right)}^{1/2}\right]}^{-1}$ Posselt and Bishop multiply it by the Kxp matrix ${\text{L}}_{Kxp}^{}\left[{\text{I}}_{pxp}-{({\Gamma }_{pxp}+\text{I})}^{-1/2}\right]$. 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), ${\text{L}}_{Kxp}^{}$ becomes identical to the K × K matrix ${\text{C}}_{K\times K}^{}$ and $\tilde{H}{\text{Z}}^f={\text{E}}_{pxK}{\Gamma }_{K\times K}^{1/2}{\text{C}}_{K\times K}^T$ becomes the concise svd of $\tilde{H}{\text{Z}}^f$. In this case, E_pxK is efficiently given by $\tilde{H}{\text{Z}}^f{\text{C}}_{K\times K}^{}{\Gamma }_{K\times K}^{-1/2}={\text{E}}_{pxK}$ and (7) becomes

$\begin{array}{l}{X}^a=\{I-Z^f{C}_{K\times K}\left[I_{K\times K}-{\left({\Gamma }_{K\times K}+I\right)}^{-1/2}\right]\\ {\Gamma }_{K\times K}^{-1}{C}_{K\times K}^T{\left(\tilde{H}{Z}^f\right)}^T\tilde{H}\}{X}^f& (9)\end{array}$

Dividing Equation (9) by $\sqrt{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]Γ⁻¹ as follows

$\begin{array}{l}\begin{array}{l}{\Gamma }^{-1}\left[\text{I}-{(\Gamma +\text{I})}^{-1/2}\right]\\ ={\Gamma }^{-1}\left[\text{I}-{(\Gamma +\text{I})}^{-1/2}\right]\left[\text{I}+{(\Gamma +\text{I})}^{-1/2}\right]\\ {\left[\text{I}+{(\Gamma +\text{I})}^{-1/2}\right]}^{-1}\\ ={\Gamma }^{-1}\left[\text{I}-{(\Gamma +\text{I})}^{-1}\right]{\left[\text{I}+{(\Gamma +\text{I})}^{-1/2}\right]}^{-1}\\ ={\Gamma }^{-1}\left[(\Gamma +\text{I}){(\Gamma +\text{I})}^{-1}-{(\Gamma +\text{I})}^{-1}\right]\\ {\left[\text{I}+{(\Gamma +\text{I})}^{-1/2}\right]}^{-1}\\ \end{array}\\ \begin{array}{l}={\Gamma }^{-1}\left[\Gamma {(\Gamma +\text{I})}^{-1}\right]{\left[\text{I}+{(\Gamma +\text{I})}^{-1/2}\right]}^{-1}\\ ={\Gamma }^{-1}\Gamma {\left[(\Gamma +\text{I})\left[\text{I}+{(\Gamma +\text{I})}^{-1/2}\right]\right]}^{-1}\\ ={\left[(\Gamma +\text{I})+{(\Gamma +\text{I})}^{1/2}\right]}^{-1}\end{array}& (10)\end{array}$

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

$\begin{array}{l}{\text{X}}^a=\left\{\text{I}-{\text{Z}}^f\text{C}{\left[(\Gamma +\text{I})+{(\Gamma +\text{I})}^{1/2}\right]}^{-1}{\text{C}}_{}^T{\left(\tilde{H}{\text{Z}}^f\right)}^T\tilde{H}\right\}{\text{X}}^f\\ =\left\{\text{I}-{\text{Z}}^f{\left[{\text{C}}_{}(\Gamma +\text{I}){\text{C}}_{}^T+{\text{C}}_{}{(\Gamma +\text{I})}^{1/2}{\text{C}}_{}^T\right]}^{-1}{\left(\tilde{H}{\text{Z}}^f\right)}^T\tilde{H}\right\}{\text{X}}^f\\ =\left\{\text{I}-{\text{Z}}^f{\left[\left(\text{C}\Gamma {\text{C}}_{}^T+\text{I}\right)+{\text{C}}_{}{(\Gamma +\text{I})}^{1/2}{\text{C}}_{}^T\right]}^{-1}{\left(\tilde{H}{\text{Z}}^f\right)}^T\tilde{H}\right\}{\text{X}}^f& (11)\end{array}$

(iii) But

$\begin{array}{ll}{\left[\left({\text{C}}_{}\Gamma {\text{C}}_{}^T+\text{I}\right)\right]}^{1/2}={\left[{\text{C}}_{}(\Gamma +\text{I}){\text{C}}_{}^T\right]}^{1/2}=\text{C}{(\Gamma +\text{I})}^{1/2}{\text{C}}_{}^T& (12)\end{array}$

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

$\begin{array}{ll}\begin{array}{l}{X}^a=\left\{I-Z^f{\left[\left(C\Gamma C_{}^T+I\right)+{\left(C\Gamma C_{}^T+I\right)}^{1/2}\right]}^{-1}{\left(\tilde{H}{Z}^f\right)}^T\tilde{H}\right\}{X}^f\hfill \\ =\{I-Z^f[\left({\left(\tilde{H}{Z}^f\right)}^T\tilde{H}{Z}^f+I\right)\hfill \\ {+{\left({\left(\tilde{H}{Z}^f\right)}^T\tilde{H}{Z}^f+I\right)}^{1/2}]}^{-1}{\left(\tilde{H}{Z}^f\right)}^T\tilde{H}\}{X}^f\hfill \end{array}& (13)\end{array}$

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 ${\gamma }_i^{max}$ and ${\gamma }_i^{min}$ to, respectively, denote the maximum and minimum of the eigenvalues listed in the eigenvalue matrix Γ_pxp. The condition number of ${(\tilde{H}\text{Z})}^T\tilde{H}\text{Z}$ is $\kappa \left[{(\tilde{H}\text{Z})}^T\tilde{H}\text{Z}\right]=\frac{{\gamma }_i^{max}}{{\gamma }_i^{min}}$. Because ${\gamma }_i^{min}$ can be zero, $\kappa \left[{(\tilde{H}\text{Z})}^T\tilde{H}\text{Z}\right]$ can be infinite. In contrast, $\kappa \left[\left({\left(\tilde{H}\text{Z}\right)}^T\tilde{H}\text{Z}+\text{I}\right)+{\left({\left(\tilde{H}\text{Z}\right)}^T\tilde{H}\text{Z}+\text{I}\right)}^{1/2}\right]=\frac{\left({\gamma }_i^{max}+1\right)+{\left({\gamma }_i^{max}+1\right)}^{1/2}}{\left({\gamma }_i^{min}+1\right)+{\left({\gamma }_i^{min}+1\right)}^{1/2}}$ is bounded above by $\left[({\gamma }_i^{max}+1)+{({\gamma }_i^{max}+1)}^{1/2}\right]/2$. However, note that the matrix $\left[{(\tilde{H}\text{Z})}^T\tilde{H}\text{Z}+\alpha \text{I}\right]$ has the eigenvalue decomposition

$\begin{array}{ll}\left[{(\tilde{H}\text{Z})}^T\tilde{H}\text{Z}+\alpha \text{I}\right]=\text{C}\Lambda {\text{C}}^T=\text{C}(\Gamma +\alpha \text{I}){\text{C}}^T& (14)\end{array}$

and hence has $\kappa \left[{(\tilde{H}\text{Z})}^T\tilde{H}\text{Z}+\alpha \text{I}\right]=\frac{{\gamma }_i^{max}+\alpha }{{\gamma }_i^{min}+\alpha }$ which is bounded above by $\frac{{\gamma }_i^{max}}{\alpha }+1$ where α is a positive scalar. Hence, α can be chosen to create a matrix that is better conditioned than $\left[\left({\left(\tilde{H}\text{Z}\right)}^T\tilde{H}\text{Z}+\text{I}\right)+{\left({\left(\tilde{H}\text{Z}\right)}^T\tilde{H}\text{Z}+\text{I}\right)}^{1/2}\right]$. Once the eigen decomposition CΛC^T 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