A note on the formulation of the Ensemble Adjustment Kalman Filter
AA note on the formulation of the Ensemble Adjustment KalmanFilter
Ian GroomsDepartment of Applied Mathematics, University of Colorado, Boulder, Colorado, USA80309
Abstract
The ensemble adjustment Kalman filter (EAKF; Anderson, 2001) is one of the earliest ensemblesquare root filters. This note clarifies the correct formulation of the EAKF, which depends on a carefultreatment of an eigen-decomposition of one of the matrices involved in the formulation.
Ensemble Kalman filters (EnKFs) are widely used for high-dimensional, nonlinear, non-Gaussian data as-similation. EnKFs use an ensemble to approximate the forecast mean and covariance, and then update thesein a manner consistent with the classical Kalman Filter formulas (Evensen, 2009). The ensemble adjustmentKalman filter (EAKF; Anderson, 2001) is one of the earliest EnKFs. The EAKF provides a deterministicupdate of the forecast ensemble, unlike the perturbed-observation version of the EnKF (Burgers et al., 1998,Houtekamer and Mitchell, 1998); as such it is a member of the class of ensemble square root filters (Tippettet al., 2003). The presentation of the EAKF in Anderson (2001) and Tippett et al. (2003) suffers from alack of clarity potentially leaving the practitioner confused as to the correct implementation of the EAKF aswell as confused about whether the EAKF is in fact consistent with the Kalman filter update formula. Thisnote clarifies the correct formulation of the EAKF, and presents a novel detail required to ensure consistencyof the update formula with the Kalman filter covariance update. No clarification is needed in the serialassimilation algorithm of Anderson (2003), which bypasses the matrix algebra of the EAKF.We follow the notation from Tippett et al. (2003). Denote the forecast ensemble by x f , · · · , x fm where m is the ensemble size. The ensemble forecast covariance matrix is P f = Z f Z fT where Z f is the scaledensemble pertubation matrix Z f = 1 √ m − (cid:104) x f − µ f , · · · , x fm − µ f (cid:105) (1)and µ f is the ensemble forecast mean. Denote the observation matrix by H and the observation errorcovariance matrix by R .The key to the EAKF is the formula for updating the ensemble perturbations in such a way that theyare consistent with the KF covariance update formula P a = ( I − KH ) P f (2)where K is the Kalman gain matrix K = P f H T (cid:16) HP f H T + R (cid:17) − . (3)Defining V = (cid:16) HZ f (cid:17) T and inserting P f = Z f Z fT , the KF update formula becomes P a = Z f (cid:20) I − V (cid:16) V T V + R (cid:17) − V T (cid:21) Z fT . (4)1 a r X i v : . [ s t a t . M E ] J un he Woodbury matrix identity can be applied to write this as P a = Z f (cid:104) I + VR − V T (cid:105) − Z fT . (5)The EAKF uses an adjustment matrix A to update the scaled ensemble perturbation matrix as Z a = AZ f in such a way that P a = Z a Z aT exactly satisfies Eqn. (5). The EAKF adjustment matrix presented byTippett et al. (2003) is A = Z f C ( I + Γ ) − / G − F T (6)where Z f = FGU T is the singular value decomposition (SVD) of Z f and VR − V T = CΓC T is theeigenvalue decomposition of VR − V T (cf. Eq. (19) of Tippett et al. (2003)).The use of the notation G − in (6) implies that the expression uses the version of the SVD that producesa square, invertible matrix G . If the state dimension is n and the rank of Z f is r then Z f is n × m , and theSVD produces G of size r × r , F of size n × r , and U of size m × r . This is precisely the interpretation givenin a footnote by Tippett et al. (2003). In this interpretation the only possible way for (6) to be meaningful isfor r = m , because otherwise the matrix ( I + Γ ) − / , which is m × m , can’t be multiplied by G − , which is r × r . But the rank of Z f is r ≤ min { n, m − } , so (6) is a mathematical impossibility, and the practitioneris left uncertain the correct formulation of the EAKF adjustment matrix.The situation can be fixed by using the version of the SVD that produces G of size r × m , F of size n × r ,and U of size m × m , and by replacing G − in (6) by the Moore-Penrose pseudoinverse G † . Here G † is an m × r diagonal matrix whose diagonal elements are 1 /σ , . . . , /σ r where σ i are the singular values of Z f that appear on the diagonal of G . With this substitution the ensemble adjustment matrix becomes A = Z f C ( I + Γ ) − / G † F T . (7)This formulation is implied by the discussion of rank deficiency in appendix A of Anderson (2001), thoughthe consistency of this adjustment matrix with the KF covariance update formula is not fully demonstratedwhen G is not invertible.We next show that this update formula leads to an analysis ensemble covariance matrix that exactlysatisfies (2), as long as the matrices C and Γ are constructed carefully. We begin by simplifying theexpression for the scaled analysis ensemble perturbation matrix Z a = Z f C ( I + Γ ) − / G † F T FGU T = Z f C ( I + Γ ) − / G † GU T . (8)The matrix G † G in this expression is an m × m diagonal matrix whose first r diagonal entries are 1, andwhose remaining entries are all 0. The product ( I + Γ ) − / G † G is a diagonal matrix whose first r entriesare the same as those of ( I + Γ ) − / , and whose remaining entries are all 0. Denoting this product by( I + Γ ) − / r we plug in to the update formula to find P a = Z a Z aT = Z f C ( I + Γ ) − r C T Z fT . (9)The correct expression differs only in the r subscript, meaning that in this formulation of EAKF the last m − r elements on the diagonal of ( I + Γ ) − r are zero when they should not be. This suggests that theEAKF produces a posterior covariance that is under-dispersed compared to the correct posterior covariance.Nevertheless, this version of the EAKF still does produce the correct posterior covariance if implementedcarefully, as we now show.First note that VR − V T = Z fT H T R − HZ f , so the null space of Z f is equal to the null space of VR − V T . This implies that there are m − r columns of C (the columns of C are orthonormal eigenvectorsof VR − V T ) that are in the null space of Z f . Multiplying any one of these columns by Z f will yield , sothe matrix Z f C that appears in (7) will have m − r columns equal to . If the eigenvalue decomposition of VR − V T is constructed so that these columns of C are the last m − r columns of C , then we have Z f C ( I + Γ ) − / r = Z f C ( I + Γ ) − / (10)2ecause the m − r trailing zeros on the diagonal of ( I + Γ ) − / r end up multiplying columns of Z f C thatare already . If, on the other hand, the m − r eigenvectors in the null space of VR − V T are not arrangedas the last m − r columns of C , then the m − r trailing zeros on the diagonal of ( I + Γ ) − / r end upmultiplying nonzero columns of Z f C , causing erroneous cancellations and leading to an under-dispersedposterior ensemble. If the columns of C are arranged correctly, then we have P a = Z a Z aT = Z t C ( I + Γ ) − r C T Z fT = Z t C ( I + Γ ) − C T Z fT = Z f (cid:104) I + VR − V T (cid:105) − Z fT (11)i.e. the EAKF update exactly satisfies (2).The presentation of the EAKF in appendix A of Anderson (2001) avoids this potential problem by usingthe language ‘singular value decomposition’ for the eigenvalue decomposition of VR − V T . Since VR − V T is a symmetric non-negative definite matrix, the SVD is the same as the eigenvalue decomposition. Thecrucial difference is that the convention for the SVD is to always arrange vectors corresponding to the nullspace as the last columns of the eigenvector matrix C . Eigenvalue solvers often do not follow this convention,and scatter the null space eigenvectors randomly in the columns of C .In summary, the correct EAKF update formula is Z a = AZ f , A = Z f C ( I + Γ ) − / G † F T (12)where Z f = FGU T is the singular value decomposition (SVD) of Z f and VR − V T = CΓC T is theeigenvalue decomposition of VR − V T , with the specific requirement that eigenvectors of VR − V T thatare in the null space of VR − V T must be arranged as the final columns of C . The purpose of this note isto (i) fully demonstrate the consistency of the EAKF of Anderson (2001), (ii) clarify the correct treatmentof the pseudoinverse in the discussion of Tippett et al. (2003), and (iii) point out a potential pitfall in theimplementation of the EAKF that could lead to inconsistency and an under-dispersed analysis ensemble ifan eigenvalue decomposition is not treated carefully. References
J. L. Anderson. An ensemble adjustment Kalman filter for data assimilation.
Mon. Wea. Rev. , 129(12):2884–2903, 2001.J. L. Anderson. A local least squares framework for ensemble filtering.
Mon. Wea. Rev. , 131(4):634–642,2003.G. Burgers, P. Jan van Leeuwen, and G. Evensen. Analysis scheme in the ensemble Kalman filter.
Mon.Wea. Rev. , 126(6):1719–1724, 1998.G. Evensen.
Data Assimilation: The Ensemble Kalman Filter . Springer, 2009.P. L. Houtekamer and H. L. Mitchell. Data assimilation using an ensemble Kalman filter technique.
Mon.Wea. Rev. , 126(3):796–811, 1998.M. K. Tippett, J. L. Anderson, C. H. Bishop, T. M. Hamill, and J. S. Whitaker. Ensemble square rootfilters.