arXiv:1605.08875v1 [math.ST] 28 May 2016

Computer Science Technical
Report CSTR-2/2016
May 31, 2016
Elias D. Niño, Adrian Sandu and
Xinwei Deng
An Ensemble Kalman Filter
Implementation Based on Modified
Cholesky Decomposition for Inverse
Covariance Matrix Estimation
Computational Science Laboratory
Computer Science Department
Virginia Polytechnic Institute and State University
Blacksburg, VA 24060
Phone: (540)-231-2193
Fax: (540)-231-6075
Email: sandu@cs.vt.edu
Web: http://csl.cs.vt.edu

Innovative Computational Solutions

An Ensemble Kalman Filter Implementation
Based on Modified Cholesky Decomposition
for Inverse Covariance Matrix Estimation
Elias D. Nino
Computational Science Laboratory
Department of Computer Science
Virginia Tech
Blacksburg, VA 24060, USA, E-mail: enino@vt.edu
Adrian Sandu
Computational Science Laboratory
Department of Computer Science
Virginia Tech
Blacksburg, VA 24060, USA, E-mail: asandu7@vt.edu
Xinwei Deng
Department of Statistics
Virginia Tech
Blacksburg, VA 24060, USA, E-mail: xdeng@vt.edu
May 31, 2016
Abstract
This paper develops an efficient implementation of the ensemble
Kalman filter based on a modified Cholesky decomposition for inverse
covariance matrix estimation. This implementation is named EnKFMC. Background errors corresponding to distant model components
with respect to some radius of influence are assumed to be conditionally independent. This allows to obtain sparse estimators of the inverse

1

background error covariance matrix. The computational effort of the
proposed method is discussed and different formulations based on various matrix identities are provided. Furthermore, an asymptotic proof
of convergence with regard to the ensemble size is presented. In order
to assess the performance and the accuracy of the proposed method,
experiments are performed making use of the Atmospheric General
Circulation Model SPEEDY. The results are compared against those
obtained using the local ensemble transform Kalman filter (LETKF).
Tests are performed for dense observations (100% and 50% of the
model components are observed) as well as for sparse observations
(only 12%, 6%, and 4% of model components are observed). The results reveal that the use of modified Cholesky for inverse covariance
matrix estimation can reduce the impact of spurious correlations during the assimilation cycle, i.e., the results of the proposed method are
of better quality than those obtained via the LETKF in terms of root
mean square error.
Keywords: Modified Cholesky decomposition, background error covariance estimation, spurious correlations, ensemble Kalman filter.

1

Introduction

The goal of sequential data assimilation is to estimate the true state of a
dynamical system xtrue ∈ n×1 using information from numerical models,
priors, and observations. A numerical model captures (with some approximation) the physical laws of the system and evolves its state forward in time
[CJAS10b]:

R

xk = Mtk−1 →tk (xk−1 ) ∈

Rn×1,

for x ∈

Rn×1,

(1)

where n is the dimension of the model state, k denotes time index, and M
can represent, for example, the dynamics of the ocean and/or atmosphere.
A prior estimation xbk ∈ n×1 of xtrue
is available, and the prior error ν is
k
usually assumed to be normally distributed:

R

xbk − xtrue = ν k ∼ N (0, Bk ) ∈

R

Rn×1,

(2)

where Bk ∈ n×n is the background error covariance matrix. Noisy observations (measurements) of the true state yk ∈ m×1 are taken, and the

R

2

observation errors  are usually assumed to be normally distributed:

yk − H xtrue
= k ∼ N (0, Rk ) ∈ m×1 ,
k

R

R

(3)

R

where m is the number of observed components, H : n×1 → m×1 is the
observation operator, and Rk ∈ m×m is the data error covariance matrix.
Making use of Bayesian statistics and matrix identities, the assimilation
of the observation (3) is performed as follows:

R

Rn×1,


−1 

xak = xbk + Bk · HTk · Hk · Bk · HTk + Rk
· yk − H(xbk ) ∈
h
i

−1
Ak = I − Bk · HTk · Rk + Hk · Bk · HTk
· Hk · Bk ∈ n×n ,

R

(4)

R

where Hk ≈ H0 (xbk ) ∈ m×n is a linear approximation of the observational
operator, and xak ∈ n×1 is the analysis (posterior) state.
According to equation (4) the elements of Bk determine how the information about the observed model components contained in the innovations
yk − H(xbk ) ∈ m×1 is distributed to properly adjust all model components,
including the unobserved ones. Thus, the successful assimilation of the observation (3) will rely, in part, on how well the background error statistics
are approximated.
In the context of ensemble based methods, an ensemble of model realizations
h
i
b[1]
b[2]
b[N ]
Xbk = xk , xk , . . . , xk ens ∈ n×Nens ,
(5)

R

R

R

is used in order to estimate the unknown moments of the background error
distribution:
xbk

Nens
1 X
b[i]
=
·
x ∈
Nens i=1 k

Bk ≈ Pb =

1
Nens − 1

Rn×1,

· Ubk · Ubk

T

(6a)

∈

Rn×n,

b[i]

(6b)

R

where Nens is the number of ensemble members, xk ∈ n×1 is the i-th
ensemble member, xbk ∈ n×1 is the background ensemble mean, Pbk is the
background ensemble covariance matrix, and Uk ∈ n×Nens is the matrix of
member deviations:

R

R

Ubk = Xbk − xbk · 1TNens ∈
3

Rn×N

ens

.

(6c)

One attractive feature of Pbk is its flow-dependency which allows to approximate the background error correlations based on the dynamics of the numerical model (1). However, in operational data assimilation, the number
of model components is much larger than the number of model realizations
n  Nens and therefore Pbk is rank-deficient. Spurious correlations (e.g., correlations between distant model components in space) can degenerate the
quality of the analysis corrections. One of the most succesful EnKF formulations is the local ensemble transform Kalman filter (LETKF) in which the
impact of spurious analysis corrections is avoided by making use of local domain analyses. In this context, every model component is surrounded by a
box of a prescribed radius, and then the assimilation is performed within every local box. In this case the background error correlations are provided by
the local ensemble covariance matrix. The local analyses are mapped back
onto the global domain to obtain the global analysis state. Nevertheless,
when sparse observational networks are considered many boxes can contain
no observations, in which case the local analyses coincide with the background. The local box sizes can be increased in order to include observations
within the local domains, in which case local analysis corrections can be impacted by spurious correlations. Moreover, in practice, the size of local boxes
can be still larger than the number of ensemble members and therefore, the
local sample covariance matrix can be rank-deficient.
In order to address the above issues this paper proposes a better estimation of the inverse background error covariance matrix B−1 obtained via
a modified Cholesky decomposition. By imposing conditional independence
between errors in remote model components we obtain sparse approximations
of B−1 .
This paper is organized as follows. In Section 2 ensemble based methods
and the modified Cholesky decomposition are introduced. Section 3 discusses
the proposed ensemble Kalman filter based on a modified Cholesky decomposition for inverse covariance matrix estimation; a theoretical convergence
of the estimator in the context of data assimilation as well as its computational effort are discussed. Section 4 presents numerical experiments using
the Atmospheric General Circulation Model SPEEDY; the results of the new
filter are compared against those obtained by the local ensemble transform
Kalman filter. Future work is discussed in Section 4.4 and conclusions are
drawn in Section 5.

4

2

Background

The ensemble Kalman filter is a sequential Monte Carlo method for state
and parameter estimation of non-linear models such as those found in atmospheric and oceanic sciences [SO08, Eve03, GMC+ 06]. The EnKF popularity
is due to its basic theoretical formulation and its relative ease of implementation [Eve03]. Given the background ensemble (5) EnKF builds the analysis
ensemble as follows:


(7a)
Xa = Xb + Pb · HT · R + H · Pb · HT · ∆ ∈ n×Nens ,

R

where:
∆ = Ys − H(Xb ) ∈

Rm×N

ens

,

(7b)

R

and the matrix of perturbed observations Ys ∈ m×Nens is:


Ys = y + [1] , y + [2] , . . . , y + [Nens ] ∈ m×Nens ,

R

[i] ∼ N (0, R) ,

1 ≤ i ≤ Nens .

(7c)

For ease of notation we have omitted the time index superscripts.
The use of perturbed observations (7c) during the assimilation provides
asymptotically correct analysis-error covariance estimates for large ensemble
sizes and makes the formulation of the EnKF statistically consistent [WT02].
However, it also has been shown that the inclusion of perturbed observations
introduces sampling errors in the assimilation [JFW14, And12].
One of the important problems faced by current ensemble based methods is that spurious correlations between distant components in the physical
space lead to spurious analysis corrections. Better approximations of the
background error covariance matrix are proposed in the literature in order to
alleviate this problem. A traditional approximation of B is the Hollingworth
and Lonnberg method [HL86] in which the difference between observations
and background states are treated as a combination of background and observations errors. However, this method provides statistics of background
errors in observation space, and requires dense observing networks (not the
case in practice). Another method has been proposed by Benedetti and
Fisher [BF07] based on forecast differences in which the spatial correlations
of background errors are assumed to be similar at 24 and 48 hours forecasts.
This method can be efficiently implemented in practice, however, it does not
5

perform well in data-sparse regions, and the statistics provided are a mixture
of analysis and background errors. Another way to reduce the impact of spurious correlations is based on adaptive modeling [Ler07]. In this context, the
model learns and changes with regard to the data collected (i.e., parameters
values and model structures). This allows to calibrate, in time, the error
subspace rank (i.e., number of empirical orthogonal functions used in the
assimilation process), the tapering parameter (i.e., local domain sizes), and
the ensemble size, among others. Yet another method based on error subspace statistical estimation is proposed in [LR99]. This approach develops an
evolving error subspace, of variable size, that targets the processes where the
dominant errors occur. Then, the dominant errors are minimized in order
to estimate the best model state trajectory with regard to the observations.
We proposed approximations based on autoregressive error models [CSCC07]
and using hybrid subspace techniques.[CJAS10a].
Covariance matrix localization artificially reduces correlations between
distant model components via a Schur product with a localization matrix
Π ∈ n×n :
b b = Π ◦ Pb ∈ n×n
P
(8)

R

R

R

b b ∈ n×n in the EnKF analysis equation (7a).
and then Pb is replaced by P
The entries of Π decrease with the distance between model components depending on the radius of influence ζ:


π (mi , mj )
, for 1 ≤ i ≤ j ≤ n ,
(9)
{Π}i,j = exp −
f (ζ)
where π (mi , mj ) represents the physical distance between the model components mi and mj while, f (ζ) is a function of ζ (e.g., f (ζ) = 2 · ζ 2 ). The
exponential decay allows to reduce the impact of innovations between distant model components. The use of covariance matrix localization alleviates
the impact of sampling errors. However, the explicit computation of Π (and
even Pb ) is prohibitive owing to numerical model dimensions. Thus, domain localization methods [Bue11, Kep00] are commonly used in the context
of operational data assimilation. One of the best EnKF implementations
based on domain localization is the local ensemble transform Kalman filter
(LETKF) [OHS+ 04b]. In the LETKF the analysis increments are computed
in the space spanned by the ensemble perturbations Ub defined in (6c). An
approximation of the analysis covariance matrix in this space reads:


b a = (Nens − 1) · I + QT · R−1 · Q −1 ∈ Nens ×Nens ,
P
(10a)

R

6

R

where Q = H · Ub ∈ m×Nens and I is the identity matrix consistent with
the dimension. The analysis increments in the subspace are:


b a · QT · R−1 · y − H(xb ) ∈ Nens ×1 ,
(10b)
αa = P

R

from which an estimation of the analysis mean in the model space can be
obtained:
xa = xb + Ub · αa ∈

Rn×1.

Finally, the analysis ensemble reads:
h
i1/2
b
a
a
T
a
b
∈
x = x · 1Nens + U · (Nens − 1) · P

(10c)

Rn×N

ens

.

(10d)

The domain localization in the LETKF is performed as follows: each
model component is surrounded by a local box of radius ζ. Within each
local domain the analysis equations (10) are applied, and therefore a local
analysis component is obtained. All local analysis components are mapped
back onto the model space to obtain the global analysis state. Local boxes
for different radii are shown in Figure 1. The local sample covariance matrix
(6b) is utilized as the covariance estimator of the local B. This can perform
well when small radii ζ are considered during the assimilation step. However,
for large values of ζ, the analysis corrections can be impacted by spurious
correlations since the local sample covariance matrix can be rank deficient.
Consequently, the local analysis increments can perform poorly.

(a) ζ = 1

(b) ζ = 2

(c) ζ = 3

Figure 1: Local domains for different radii of influence ζ. The red dot is the
model component to be assimilated, blue components are within the scope
of ζ, and black model components are unused during the local assimilation
process.
7

There is an opportunity to reduce the impact of sampling errors by improving the background error covariance estimation. We achieve this by
making use of the modified Cholesky decomposition for inverse covariance
matrix estimation [BL08]. Consider a sample of Nens Gaussian random vectors:


S = s[1] , s[2] , . . . , s[Nens ] ∈ n×Nens ,

R

with statistical moments:
s[j] ∼ N (0n , Q) , for 1 ≤ j ≤ Nens ,

R

R

where s[j] ∈ n×1 denotes the j-th sample. Denote by x[i] ∈ Nens ×1 the
vector holding the i-th component across all the samples (the i-th row of
S, transposed). The modified Cholesky decomposition arises from regressing
each component on his predecessors according to some component ordering:
x[i] =

i−1
X

x[j] · βi,j + ε[i] ∈

RN

ens ×1

,

2 ≤ i ≤ n,

(11)

j=1

where x[j] is the j-th model component which precedes x[i] for 1 ≤ j ≤ i − 1,
ε[1] = x[1] , and ε[i] ∈ Nens ×1 is the error in the i-th component regression
for i ≥ 2. Likewise, the coefficients βi,j in (11) can be computed by solving
the optimization problem:

R

β [i] = arg min x[i] − Z[i] · β
β

2
2

(12)

where
Z[i] =
β [i] =

R(i−1)×N
[βi,1 , βi,2 , . . . , βi,i−1 ]T ∈ R(i−1)×1 .
 [1] [2]
T
x , x , . . . , x[i−1] ∈

ens

,

2 ≤ i ≤ n,

The regression coefficients form the lower triangular matrix


 −βi,j for 1 ≤ j < i,

b
1 for j = i,
T
=
1 ≤ i ≤ n,
i,j


0 for j > i,

8

(13a)


b
b ∈ Rn×n . The
where T
denotes the (i, j)-th component of matrix T
i,j
empirical variances cd
ov of the residuals ε[i] form the diagonal matrix:
b = diag cd
D
ov(ε[i] ) = diag
1≤i≤n

where

n o
b
D

1,1

1≤i≤n

Nens
X


1



Nens − 1

!
2
ε[i] j

∈ Rn×n .

(13b)

j=1


= cd
ov x[1] . Then an estimate of Q−1 can be computed as

follows:
b −1 = T
bT · D
b −1 · T
b∈
Q

Rn×n,

(14a)

or, by basic matrix algebra identities the estimate of Q reads:
b =T
b −1 · D
b ·T
b −T ∈
Q

Rn×n .

(14b)

b −1 is strictly related to the structure of T.
b This
Note that the structure of Q
−1
can be exploited in order to obtain sparse estimators of Q by imposing that
b are zero. This is important for high dimensional probability
some entries of T
b or Q
b −1 is prohibitive. The
distributions where the explicit computation of Q
b can be justified as follows: when two components are
zero components in T
b −1 is zero. In the
conditionally independent their corresponding entry in Q
context of data assimilation, the conditional independence of background
errors between different model components can be achieved by making use of
domain localization. We can consider zero correlations between background
errors corresponding to model components located at distances that exceed
a radius of influence ζ. In the next section we present an ensemble Kalman
filter implementation based on modified Cholesky decomposition for inverse
covariance matrix estimation.

3

Ensemble Kalman Filter Based On Modified Cholesky Decomposition

In this section we discuss the new ensemble Kalman filter based on modified
Cholesky decomposition for inverse covariance matrix estimation ( EnKFMC).

9

3.1

Estimation of the inverse background covariance

The columns of matrix (6c)


Ub = ub[1] , ub[2] , . . . , ub[Nens ] ∈

Rn×N

ens

can be seen as samples of the (approximately normal) distribution:
xb[j] − xb = ub[j] ∼ N (0, B) , for 1 ≤ j ≤ Nens ,

R

and therefore, if we let x[i] ∈ Nens ×1 in (11) to be the vector formed by
the i-th row of matrix (6c), for 1 ≤ i ≤ n, according to equations (14), an
estimate of the inverse background error covariance matrix reads:
b −1 = T
bT · D
b −1 · T
b∈
B−1 ≈ B

Rn×n,

(15a)

and therefore:
b =T
b −1 · D
b ·T
b −T ∈
B≈B

Rn×n .

(15b)

b −1 depends on that of T.
b If we
As we mentioned before, the structure of B
assume that the correlations between model components are local, and there
are no correlations outside a radius of influence ζ, we obtain lower-triangular
b Consequently, the resulting B
b −1 will also be sparse,
sparse estimators of T.
b will be localized. Since the regression (11) is performed only on the
and B
predecessors of each model component, an ordering (labeling) must be set
b Since we work with
on the model components prior the computation of T.
gridded models we consider column-major and row-major orders. They are
illustrated in Figure 2 for a two-dimensional domain. Figure 3 shows the local
domain and the predecessors of the model component 6 when column-major
order is utilized.

(a) Column-major order

(b) Row-major order

Figure 2: Row-major and column-major ordering for a 4 × 4 domain. The
total number of model components is 16.
10

(a) In blue, local box for the
model component 6 when ζ = 2.

(b) In blue, predecessors of the
model component 6 for ζ = 2.

Figure 3: Local model components (local box) and local predecessors for the
model component 6 when ζ = 2. Column-major ordering is utilized to label
the model components.
b −1 proceeds as follows:
The estimation of B

R

1. Form the matrix Z[i] ∈ pi ×Nens with the predecessors of the i-th model
component:

T
Z[i] = x[q(i,1)] , x[q(i,2)] , . . . , x[q(i,pi )] ∈ pi ×Nens ,
(16)

R

where x[e] is the e-th row of matrix (6c), pi is the number of predecessors
of component i, and 1 ≤ q(i, j) ≤ n is the index (row of matrix (6c))
of the j-th predecessor of the i-th model component.
2. For the i-th model components the regression coefficients are obtained
as follows:
x[i] =

pi
X

βi,j · x[q(i,j)] + ε[i] ∈

RN

ens ×1

.

j=1

Rp ×1 by solving

For 2 ≤ i ≤ n, compute β [i] = [βi,1 , βi,2 , . . . , βi,pi ] ∈
the optimization problem (12) with Z[i] given by (16).
3. Build the matrices

b
T
= −βi,j
i,q(i,j)

for 1 ≤ j ≤ pi , 1 < i ≤ n ;


b
T

i

i,i

= 1,

b according to equation (13b). Note that the number of non-zero
and D
b equals the number of predecessors pi .
elements in the i-th row of T
11

Note that the solution of the optimization problem (12) can be obtained
as follows:

−1
β [i] = Z[i] · Z[i] T
· Z[i] · x[i]
(17)
and since the ensemble size can be smaller than the number of model components, Z[i] · Z[i] T ∈ pi ×pi can be rank deficient. To overcome this situation,
regularization of the zero singular values of Z[i] · Z[i] T can be used. One
possibility is Tikhonov regularization [GHO99, KHE07, Mai94]:
n
o
2
2
2
β [i] = arg min x[i] − Z[i] · β 2 + λ · β 2
(18)

R

β

R

where λ ∈ . In our context the best choice for λ relies on prior knowledge
of the background and the observational errors [Neu88]. Another approach
to regularization is to use a truncated singular value decomposition (SVD)
of Z[i] :
T

Z[i] = UZ[i] · ΣZ[i] · VZ[i] ∈

R

R

Rp ×N
i

ens

,

where UZ[i] ∈ pi ×pi and VZ[i] ∈ Nens ×Nens are the right and the left singular
vectors of Z[i] , respectively. Likewise, ΣZ[i] ∈ pi ×Nens is a diagonal matrix
whose diagonal entries are the singular values of Z[i] in descending order.
The solution of (12) can be computed as follows [JB00, Huf91, Han90]:
β [i] =

R

ki
X
1
Z T
Z
· uj [i] · vj [i] · x[i]
τ
j=1 j

with

τj
τmax

≥ σr ,

(19)

where τj is the j-th singular value with corresponding right and left singular
Z
Z
vectors uj [i] ∈ pi ×1 and vj [i] ∈ pi ×1 , respectively, σr ∈ (0, 1) is a predefined
threshold, and τmax = max {τ1 , τ2 , . . . , τNens −1 }. Since small singular values
are more sensitive to the noise in x[i] , the threshold τj > τmax · σr seeks to
neglect their contributions.

R

3.2

R

Formulation of EnKF-MC

b −1 is estimated, the EnKF based on modified Cholesky decomposition
Once B
(EnKF-MC) computes the analysis using Kalman’s formula:
xa = xb + Æ · HT · R−1 · ∆ ∈
12

Rn×N

ens

,

(20a)

where Æ ∈

Rn×n is the estimated analysis covariance matrix
h
i−1
b −1 + HT · R−1 · H
Æ= B
,

R

and ∆ ∈ m×Nens is the innovation matrix on the perturbed observations
given in (7b).
Computationally-friendlier alternatives to (20c) can be obtained by making use of elementary matrix identities:
h
i
a
−1
b
T
−1
s
b
x = Æ · B · x + H · R · Y ∈ n×Nens ,
(20b)


b −1 · D
b 1/2 · VTb · R + V b · VTb −1 · ∆,
xa = xb + T
(20c)
B
B
B
b −1 · D
b 1/2 ∈ n×m ,
Vb = H · T

R

B

R

where Ys are the perturbed observations. The formulation (20c) is wellknown as the EnKF dual formulation, (20b) is known as the EnKF primal
formulation, and the equation (20a) is the incremental form of the primal
formulation. In the next subsection, we discuss the computational effort of
the EnKF-MC implementations (20).

3.3

Computational effort of EnKF-MC implementations

The computational cost of the different EnKF-MC implementations depend,
in general, on the model state dimension n, the number of observed components m, the radius of influence ζ, and the ensemble size Nens . Typically
[TAB+ 03] the data error covariance matrix R has a simple structure (e.g.,
block diagonal), the ensemble size is much smaller than the model dimension (n  Nens ), and the observation operator H is sparse or can be applied
efficiently. We analyze the computational effort of the formulation (20a); similar analyses can be carried out for the other formulations. The incremental
formulation can be written as follows:
xa = xb + δXa ,

R

where the analysis increments δXa ∈ n×Nens are given by the solution of
the linear system:
h
i
b −1 + RH · RT · δXa = ∆H .
B
H
13

R

R

with RH = HT · R−1/2 ∈ n×m , ∆H = HT · R−1 · ∆ ∈ n×Nens , and ∆ is
given in (7b). This linear system can be solved making use of the iterative
Sherman Morrison formula [NRSA14] as follows:
1. Compute:
(0)[i]
WZ
(0)[j]

WU
[i]

=

h

=

h

b T · D−1 · T
b
T

i−1

b T · D−1 · T
b
T

i−1

R

[j]

[i]

(21a)

[j]

(21b)

· ∆H , for 1 ≤ i ≤ Nens ,
· RH , for 1 ≤ j ≤ m .

R

where ∆H ∈ n×1 and RH ∈ n×1 denote the i and j columns of
b is a sparse unitary lower
matrices ∆H and RH , respectively. Since T
triangular matrix, the direct solution of the linear system (21) can be
obtained by making use of forward and backward substitutions. Hence,
this step can be performed with:
O (nnz · n · Nens + nnz · n · m)

(22)

long computations, where nnz denotes the maximum number of nonb this is
zero elements across all rows of T,
nnz = max {p1 , p2 , . . . , pn }
where pi is the number of predecessors of model component i, for 1 ≤
i ≤ n.
2. For 1 ≤ i ≤ m compute:
i−1
h
1
(i−1)[i]
(i−1)[i]
[i] T
(i)
·
W
,
·
W
,
with
γ
=
1
+
R
U
U
H
γ (i)
h
i
(i−1)[j]
[i] T
(i−1)[j]
= WZ
− h(i) · RH · WZ
, for 1 ≤ j ≤ Nens ,
h
i
(i−1)[k]
[i] T
(i−1)[k]
= WU
− h(i) · RH · WU
, for i + 1 ≤ k ≤ m .

h(i) =
(i)[j]

WZ

(i)[k]

WU

Note that, at each step, h(i) can be computed with n long computations, while WZ and WU can be obtained with n · Nens and n · m
long computations, respectively. This leads to the next bound for the
number of long computations:

O m · n + m · n · Nens + m2 · n .
14

Hence, the computational effort involved during the assimilation step of formulation (20a) can be bounded by:

O m · n + m · n · Nens + m2 · n + nnz · n · Nens + nnz · n · m ,
which is linear with respect to the number of model components. For dense
observational networks, when local observational operators can be approximated, domain decomposition can be exploited in order to reduce the computational effort during the assimilation cycle. This can be done as follows:
1. The domain is split in certain number of sub-domains (typically matching a given number of processors).
2. Background error correlations are estimated locally.
3. The assimilation is performed on each local domain.
4. The analysis sub-domains are mapped back onto the model domain
from which the global analysis state is obtained.
Figure 4 shows the global domain splitting for different sub-domain sizes. In
Figure 4c the boundary information needed during the assimilation step for
two particular sub-domains is shown in dashed blue lines. Note that each
sub-domain can be assimilated independently. Note that we only use domain
decomposition in order to reduce the computational effort of the proposed
implementation (and its derivations) and not in order to reduce the impact
of spurious correlations.

15

(a) Number of sub-domains 12

(b) Number of sub-domains 80

(c) Number of sub-domains 16.

Figure 4: Global domain splitting in different sub-domain sizes. Blue local
boxes reflects the boundary information utilized in order to perform local
data assimilation.

3.4

Convergence of the covariance inverse estimator

b −1 estimator in the context
In this section we prove the convergence of the B
of data assimilation.
Comment 1 (Sparse Cholesky factors and localization). The modified Cholesky
decomposition for inverse covariance matrix estimation can be seen as a form
of covariance matrix localization method in which the resulting matrix approximates the inverse of a localized ensemble covariance matrix. This process
is implicit in the resulting estimator when only a local neighborhood for each
16

model component is utilized in order to perform the local regression and to
b and D.
b Figure 5 shows an example for the Lorenz 96 [Lor05]:
estimate T

(x2 − xn−1 ) · xn − x1 + F
for j = 1
dxj 
= (xj+1 − xj−2 ) · xj−1 − xj + F for 2 ≤ j ≤ n − 1
(23)

dt

(x1 − xn−2 ) · xn−1 − xn + F
for j = n
where F is usually set to 8 to exhibit chaotic behavior and the number of
model components is n = 40. We assume B to be a sample covariance matrix
based on 105 samples, the localized ensemble covariance matrix Pb and the
b −1 are based on 80 samples. The radius of influence is ζ = 7.
estimator B
The similarities among the different Cholesky factors is evident. Even more,
along the main diagonal, the correlations decay with respect to the distance
of the model components. This is reflected in the resulting estimator of B−1
for each case. Definition 1 of covariance matrices relies on this assumpation
b
for the Cholesky factors T and T.

17

(a) Exact B−1 ≈ Pb

−1

for Nens = 105

(b) T, B−1 = TT · D−1 · T

cb −1
(c) Localized ensemble estimate P

cb −1 = TT · D · T
(d) TL , P
L
L
L

b −1
(e) Cholesky estimate B

b B
b −1 = T
bT · D
b −1 · T
b
(f) T,

Figure 5: Decay of correlations in the Cholesky factors for different approximations of B−1 .
18

We consider a two-dimensional square domain with s×s grid points. Our
proof below can be extended immediately to non-square domains, as well as
to three-dimensional domains. In our domain each space point is described
by two indices (i, j), a zonal component i and a meridional component j, for
1 ≤ i, j ≤ s. A particular case for s = 4 is shown in Figure 6a. We make
use of row-major order in order to map model grid components to the one
dimensional “index space”:
k = f (i, j) = (j − 1) · s + i,

for 1 ≤ k ≤ n.

where here, n = s2 . For a particular grid component (i, j), the resulting
b −1 . The results of labeling each model
k = f (i, j) denotes the row index in B
component in this manner can be seen in Figure 2b.

(a) Grid components (i, j)

(b) Index space f (i, j)

Figure 6: Grid distribution of model components and corresponding index
b −1 .
terms in B
To start our proof, the inverse of the (exact) background error covariance
b −1 can be written as
matrix B−1 and of the its estimator B
h
iT
h
i
−1
−1
b
b
b
b
B = I − C · D · I − C ∈ n×n
(24a)

R

and
B−1 = [I − C]T · D−1 · [I − C] ∈
b = I−T
b∈
respectively, where C
b are diagonal matrices:
and D

Rn×n,

(24b)

Rn×n and C = I − T ∈ Rn×n. Moreover, D


D = diag d21 , d22 , . . . , d2n
n
o
b = diag db2 , db2 , . . . , db2
D
1
2
n
19

where {D}i,i = d2i and

R

n o
b
D

i,i
{j}

= db2i , for 1 ≤ i ≤ n. In what follows we

denote by b
c{j} ∈ n×1 and c ∈
C, respectively, for 1 ≤ j ≤ n.

Rn×1 the j-th columns of matrices Cb and

Definition 1 (Class of matrices under consideration.). We consider the class
of covariance matrices matrices with correlations decreasing quickly:
(
U −1 (ε0 , C, α) =

B : 0 < ε0 ≤ λmin (B) ≤ λmax (B) ≤ ε−1
0 ,
max
k

n
X

(25a)
)

γk,` · {T}k,` ≤ C · ζ −α , for ζ ≤ s − 1

`=1

where B−1 = TT D−1 T, α is the decay rate (related to the dynamics of the
numerical model),


0 j − ζ ≤ q ≤ j − 1 and i − ζ ≤ p ≤ i + ζ
,
γk(i,j) ,`(p,q) =
0 q = j and i − ζ ≤ p ≤ i


1 otherwise
and the grid components (i, j) and (p, q), for 1 ≤ i, j, p, q ≤ s are related to
the (k(i,j) , `(p,q) ) matrix entry by k(i,j) = f (i, j) and ` = f (p, q).
Comment 2. The factors γk,` for the grid component (i, j) in Definition (1)
are zero inside the scope of ζ.
Theorem 2 (Error in the covariance inverse estimation). Uniformly for B ∈

−1/2(α+1)
U −1 (ε0 , C, α), if ζ ≈ Nens −1 · log n
and Nens −1 · log n = o(1),

α(α+1)/2 !
log(n)
b −1 − B−1
B
=O
(25b)
∞
Nens
where

·

∞

denotes the infinity norm (matrix or vector)

Comment 3. The factors γk,` in Theorem (2) are zero for the predecessors
of the grid component (i, j) inside the scope of ζ.
In order to prove Theorem 2, we need the following result.
20

Lemma 3. Under the conditions of Theorem 2, uniformly on U −1


 {j}
max b
c − c{j} ∞ : 1 ≤ j ≤ n = O Nens −1/2 log1/2 n , (26a)
o

n
α/(2(α+1)) 
, (26b)
max db2j − d2j : 1 ≤ j ≤ n = O Nens −1 log n
and
C

= D−1

∞

∞

= O (1) .

(26c)

The proof of Lemma 3 is based on the following results of Bickel and
Levina in [BL08].
Lemma 4. [[BL08, Lemma A.2]] Let ν [k] ∼ N (0, B) and λmax (B) ≤ ε−1
0 <
∞, for 1 ≤ k ≤ Nens . Then, if {B}i,j denotes the (i, j)-th component of B,
for 1 ≤ i ≤ j ≤ n,
"N
#
ens h
i
X
 [k]
 [k]
Prob
ν i · ν j − {B}i,j ≥ Nens · ν
(27)
k=1


≤ C1 · exp −C2 · Nens · ν 2 ,

for |ν| ≤ δ, where ν [k] i is the i-th component of the sample ν [k] , for 1 ≤
k ≤ Nens , and 1 ≤ i ≤ n. Likewise, C1 , C2 and δ depend on ε0 only.
Proof of Lemma 3. In what follows we denote by cov and cd
ov denote the
true and the empirical
covariances, respectively. In the context of EnKF we

b
have that cov U = B.
Recall that
cd
ov U

b



1

b

=P =

b

Nens − 1

bT

·U ·U

=

1
Nens − 1

·

Nens
X

T

ub[k] · ub[k] ,

k=1

and therefore



cd
ov Ub


i,j

=

1
Nens − 1

·

Nens
X


ub[k]

 b[k]
·
u
i

j

.

k=1





For ν > 0, ν [k] i · ν [k] j −{B}i,j ≥ Nens ·ν implies ν [k] i · ν [k]
(Nens − 1) · ν, and therefore by Lemma 4 we have:




cov Ub − cd
ov Ub ∞ = O Nens −1/2 · log1/2 n ,
21

j

−{B}i,j ≥
(28a)



since the entries of cov Ub − cd
ov Ub can be bounded by:


b

cov U



b

− cd
ov U



≤ Nens

i,j

−1

·

Nens
X


 b[k]
·
u
i

ub[k]

j

− {B}i,j .

k=1

Lemma 4 ensures that:
"
Prob max Nens −1 ·
i,j

Nens
X


#
ub[k]

· ub[k]


i

j

− {B}i,j ≥ ν

k=1
2


≤ C1 · n · exp −C2 · Nens · ν 2 ,
for |ν| ≤ δ. Let ν =



log n2
Nens ·C2

1/2

· M , for M arbitrary.

Since Z[i] stores the columns of Ub corresponding to the predecessors of
model component i, an immediate consequence of (28a) is




(28b)
max cov Z[i] − cd
ov Z[i] ∞ = O Nens −1/2 · log1/2 n .
i

Also,
B−1

∞

= cov Ub

−1
∞

≤ ε−1
0 .

According to equation (17),
c[i]

j

 [i]
b
c

j



o
n
−1
· Z[i] · x[i] ,
cov Z[i]
j
n
o
−1
= cd
ov Z[i]
· Z[i] · x[i] ,

=

j

therefore:
 [i]
c[i] k − b
c k
k
n
o
n
o
−1
−1
= max cov Z[i]
· Z[i] · x[i] − cd
ov Z[i]
· Z[i] · x[i]
(29)
k
k
k
nh
i
o
−1
−1
= max
cov Z[i]
− cd
ov Z[i]
· Z[i] · x[i]
k
k


1/2
−1/2
= O Nens
· log n
(30)
max



22

from which (26a) follows. Note that:
x[i] =

n
X

 [i]
γ̃i,j · b
c

j

· x[j] + b
ε[i]

j=1



⇔ cd
ov x[i] = cd
ov

n
X

!
γ̃i,j

 [i]
· b
c

γ̃i,j

 [i]
· b
c

j

· x[j] + b
ε[i]

j=1



⇔ cd
ov x[i] = cd
ov

n
X

!
j

· x[j]

 
+ cd
ov b
ε[i]

j=1
n
X

⇔ db2i = cd
ov x[i] − cd
ov


!
 [i]
c
· b

γ̃i,j

· x[j]

j

,

j=1

and similarly
d2i = cov x[i] − cov


n
X

!
γ̃i,j · c[i]


j

· x[j]

.

j=1

The claim (26b) and the first part of (26c) follow from (28a), (28b) and
(30). Since
db2i − d2i

≤
+



cov x[i] − cd
ov x[i]
cd
ov

n
X

γ̃i,j

h
c[i]
· b

γ̃i,j

 [i]
· b
c

− c[i]

i



j

j

!
· x[j]

(31)

j=1

+

cd
ov

n
X

!
j

· x[j]

− cov

n
X

!
γ̃i,j

 [i]
c
· b

j

· x[j]

j=1

j=1

where γ̃i,j = 1 − γi,j . By Lemma 4 the maximum over i of the first term is:




1/2
−1/2
max cov x[i] − cd
ov x[i] = O Nens
· log n .
i

23

The second term can be bounded as follows:
n
X

≤

2
·
γ̃i,j

j=1
n
X

2
γ̃i,j
·

h
b
c[i]

h
b
c[i]

j

j

i2


− c[i] j ·[
cov x[j]

i2


− c[i] j · cd
ov x[j]

j=1

 [i]
≤ max b
c
k



k

− c

−1

2

2

[i]
k

· max cd
ov x[i]



i

·

n
X

2
γ̃i,j

j=1



= O ζ · Nens · log n

α/2·(α+1) 
= O Nens −1 · log n

1/2·(α+1)
−1
by (26a) and B ≤ ε−1
· log n
and even
0 . Recall that ζ = Nens
more, note that:
n
X


ζ2
1
[ζ + 1]2
=
+ ζ + = O ζ2 .
=
2
2
2

2
γ̃i,j

j=1

The third term can be bounded similarly. Thus (26b) follows. Furthermore,
!
!
n
n h
i2
X
X

[i]
b
≥ ε0 ,
d2i = cov x[i] −
γ̃i,j · b
c[i] j · x[j] ≥ ε0 · 1 +
cj
j=1

i=1

and the lemma follows.
We now are ready to prove Theorem 2.
Proof of Theorem 2. We need only check that:


1/2
−1/2
−1
−1
b
B −B
= O Nens
· log (n)
∞

(32a)

and
B−1 − Φζ B−1


∞

= O ζ −α



where the entries of Φζ (B−1 ) are given by:



Φζ B−1 k,` = δk,` · B−1 k,` , for 1 ≤ k, ` ≤ n
24

(32b)

(32c)

where k = f (i, j) and ` = f (q, p) for 1 ≤ i, j, p, q ≤ s, and
(
1 j − ζ ≤ q ≤ j + ζ and i − ζ ≤ p ≤ i + ζ
δk,` =
0 otherwise
We first prove (32a). By definition,
b −1 − B−1 = T
bT · D
b −1 · T
b − TT · D−1 · T.
B

(33)

Applying the standard inequality:
bT · D
b −1 · T
bT
TT · D−1 · T − T

bT · D
b
TT − T
b · T
bT
D−D
b · T
b ·
T−T

≤
+
+

b
· T
b
· T

b
D
b T − TT
· T
b T − TT
· T

b · D−D
b
T
b · T−T
b
D
bT · D − D
b · T
b −T
T
b · T−T
b · TT − T
bT
D−D

+
+
+
+

all previous terms can be bounded making use of Lemma 3 and therefore,
(32a) follows. Likewise, for (32b), we need to note that for any matrix M,
M · MT − Φζ (M) · Φζ (M)T

∞

≤ 2· M
+

∞

· Φζ (M) − M−1

Φζ (M) − M

∞

2
∞

and by letting M = TT · D−1/2 , the theorem follows from Definition 1.

4

Numerical Experiments

In this section we study the performance of the proposed EnKF-MC implementation. The experiments are performed using the atmospheric general circulation model SPEEDY [Mol03, KMB06]. SPEEDY is a hydrostatic, spectral coordinate, spectral transform model in the vorticity-divergence form,
with semi-implicit treatment of gravity waves. The number of layers in the
SPEEDY model is 8 and the T-63 model resolution (192 × 96 grids) is used
for the horizontal space discretization of each layer. Four model variables
25

are part of the assimilation process: the temperature (K), the zonal and
the meridional wind components (m/s), and the specific humidity (g/kg).
The total number of model components is n = 589, 824. The number of
ensemble members is Nens = 94 for all the scenarios. The model state space
is approximately 6,274 times larger than the number of ensemble members
(n  Nens ).
Starting with the state of the system xref
−3 at time t−3 , the model solution
ref
x−3 is propagated in time over one year:

ref
xref
−2 = Mt−3 →t−2 x−3 .
The reference solution xref
−2 is used to build a perturbed background solution:



b
ref
b
b
ref
2
b−2 = x−2 + −2 , −2 ∼ N 0n , diag (0.05 {x−2 }i )
x
.
(34)
i

The perturbed background solution is propagated over another year to obtain
the background solution at time t−1 :

bb−2 .
xb−1 = Mt−2 →t−1 x
(35)
This model propagation attenuates the random noise introduced in (34) and
makes the background state (35) consistent with the physics of the SPEEDY
model. Then, the background state (35) is utilized in order to build an
ensemble of perturbed background states:



b[i]
b
b
b
b
2
b−1 = x−1 + −1 , −1 ∼ N 0n , diag (0.05 {x−1 }i )
x
, 1 ≤ i ≤ Nens ,(36)
i

from which, after three months of model propagation, the initial ensemble is
obtained at time t0 :


b[i]
bb[i]
x0 = Mt−1 →t0 x
.
−1
Again, the model propagation of the perturbed ensemble ensures that the
ensemble members are consistent with the physics of the numerical model.
The experiments are performed over a period of 24 days, where observations are taken every 2 days (M = 12). At time k synthetic observations are
built as follows:

2
yk = Hk · xref
k ∼ N (0m , Rk ) , Rk = diagi (0.01 {Hk xref
.
k + k ,
k }i )
26

The observation operators Hk are fixed throughout the time interval. We
perform experiments with several operators characterized by different proportions p of observed components from the model state xref
k (m ≈ p · n).
We consider four different values for p: 0.50, 0.12, 0.06 and 0.04 which represent 50%, 12 %, 6 % and 4 % of the total number of model components,
respectively. Some of the observational networks used during the experiments are shown in Figure 7 with their corresponding percentage of observed
components from the model state.
The analyses of the EnKF-MC are compared against those obtained making use of the LETKF implementation proposed by Hunt et al in [OHS+ 04b,
OHS+ 04a, OHS+ 08] . The analysis accuracy is measured by the root mean
square error (RMSE)
v
u
M
u1 X
 ref
T 

a
xk − xak · xref
(37)
RMSE = t ·
k − xk
M k=1

R

R

where xref ∈ n×1 and xak ∈ n×1 are the reference and the analysis solutions
at time k, respectively, and M is the number of assimilation times.
b −1 is σr = 0.10.
The threshold used in (19) during the computation of B
During the assimilation steps, the data error covariance matrices Rk are
used (no representativeness errors are involved during the assimilations) and
therefore. The different EnKF implementations are performed making use of
FORTRAN and specialized libraries such as BLAS and LAPACK are used
in order to perform the algebraic computations.

(a) p = 12%

(b) p = 4%

Figure 7: Observational networks for different values of p. Dark dots denote
the location of the observed components. The observed model variables are
the zonal and the meridional wind components, the specific humidity, and
the temperature.
27

4.1

Results with dense observation networks

We first consider dense observational networks in which 100% and 50% of
the model components are observed. We vary the radius of influence ζ from
1 to 5 grid points.
Figure 8 shows the RMSE values for the LETKF and EnKF-MC analyses
for different values of ζ for the specific humidity when 50% of model components are observed. When the radius of influence is increased the quality of
the LETKF results degrades due to spurious correlations. This is expected
since the local estimation of correlations in the context of LETKF is the
sample covariance matrix. For instance, for a radius of influence of 1, the
total number of local components for each local box is 36 which matches the
dimension of the local background error distribution. Now, when we compare it against the ensemble size (96 ensemble members), sufficient degrees
of freedom (95 degrees of freedom) are available in order to estimate the local
background error distribution onto the ensemble space, and consequently all
directions of the local probability error distribution are accounted during the
estimation and posterior assimilation. On the other hand, when the radius
of influence is 5, the local box sizes have dimension 484 (model components)
which is approximately 5 times larger than the ensemble size. Thus, when
the analysis increments are computed onto the ensemble space, just part
of the local background error distribution is accounted during the assimilation. Consequently, the larger the local box, the more local background error
information cannot be represented in the ensemble space.
Figure 8 shows that EnKF-MC analyses improve with increasing radius
of influence ζ. Since a dense observational network is considered during the
assimilation, when the radius of influence is increased, a better estimation
of the state of the system is obtained by the EnKF-MC. This can be seen
clearly in Figure 9, where the RMSE values within the assimilation window
are shown for the LETKF and the EnKF-MC solutions for the specific humidity variable and different values of ζ and p. The quality of the EnKF-MC
analysis for ζ = 5 is better than that of the LETKF with ζ = 1. Likewise,
when a full observational network is considered (p = 100%), the proposed implementation outperforms the LETKF implementation. EnKF-MC is able to
exploit the large amount of information contained in dense observational networks by properly estimating the local background error correlations. The
RMSE values for all model variables and different values for ζ and p are
summarized in Table 1.
28

(a) ζ = 1.

(b) ζ = 2.

(c) ζ = 3.

(d) ζ = 4.

(e) ζ = 5.

Figure 8: RMSE of specific humidity analyses with a dense observational
network. When the radius of influence ζ is increased the performance of
LETKF degrades.
29

Root Mean Square Error (RMSE)

Root Mean Square Error (RMSE)

Specific Humidity (g/Kg)

12

EnKF-MC
LETKF

10

8

6

4

1

2

3
Radius (ζ)

4

5

(a) p = 100%.

Specific Humidity (g/Kg)

12

EnKF-MC
LETKF

10

8

6

4

1

2

3
Radius (ζ)

4

5

(b) p = 50%.

Figure 9: Analysis RMSE for the specific humidity variable. The RMSE
values of the assimilation window are shown for different values of ζ and
percentage of observed components p. When the local domain sizes are increased the accuracy of the LETKF analysis degrades, while the accuracy of
EnKF-MC analysis improves.

30

Variable (units)

ζ
1
2

Zonal Wind Component (u), (m/s)

3
4
5
1
2

Meridional Wind Component (v) (m/s)

3
4
5
1
2

Temperature (K)

3
4
5
1
2

Specific Humidity (g/Kg)

3
4
5

p
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%
100%
50%

EnKF-MC
6.012 × 101
4.264 × 102
6.078 × 101
2.255 × 102
6.080 × 101
2.341 × 102
6.088 × 101
2.418 × 102
6.092 × 101
2.673 × 102
3.031 × 101
2.632 × 102
3.046 × 101
1.641 × 102
3.047 × 101
1.964 × 102
3.052 × 101
2.084 × 102
3.054 × 101
2.428 × 102
9.404 × 102
6.644 × 102
9.416 × 102
6.129 × 102
9.425 × 102
5.815 × 102
9.432 × 102
5.585 × 102
9.432 × 102
5.500 × 102
1.733 × 101
8.680 × 101
1.712 × 101
8.204 × 101
1.705 × 101
8.089 × 101
1.699 × 101
7.525 × 101
1.694 × 101
7.642 × 101

LETKF
6.394 × 101
9.825 × 102
6.820 × 101
1.330 × 103
7.969 × 101
1.124 × 103
9.687 × 101
1.072 × 103
1.190 × 102
1.017 × 103
6.418 × 101
3.247 × 102
6.597 × 101
4.138 × 102
7.565 × 101
4.418 × 102
9.332 × 101
4.832 × 102
1.151 × 102
5.029 × 102
5.078 × 102
7.059 × 102
4.112 × 102
1.138 × 103
3.447 × 102
1.389 × 103
2.939 × 102
1.355 × 103
2.554 × 102
1.104 × 103
5.427 × 101
7.602 × 101
5.669 × 101
1.045 × 102
6.630 × 101
1.298 × 102
7.344 × 101
1.431 × 102
7.617 × 101
1.458 × 102

Table 1: RMSE values for the EnKF-MC and the LETKF analyses with the
SPEEDY model and for different values for ζ and p. Dense observational
networks are considered in this experimental setting.
31

4.2

Results with sparse observation networks

For sparse observational networks, in general, the results obtained by the
EnKF-MC are more accurate than those obtained by the LETKF, as reported in the Tables 2 and 3. We vary the values of ζ from 1 to 5. Three
sparse observational networks with p = 12%, 6%, and 4%, respectively are
considered.
Figure 10 shows the RMSE values of the specific humidity analyses for
different radii of influence and 4% of the model components being observed.
The best performance of the LETKF analyses is obtained when the radius
of influence is set to 2. Note that for ζ = 1 the LETKF performs poorly,
which is expected since during the assimilation most of model components
will not have observations in their local boxes. For ζ ≥ 3 the effects of
spurious correlations degrade the quality of the LETKF analysis. On the
other hand, the background error correlations estimated by the modified
Cholesky decomposition allows the EnKF-MC formulation to obtain good
analyses even for largest radius of influence ζ = 5.
Figure 11 shows the RMSE values of the LETKF and the EnKF-MC implementations for different radii of influences and two sparse observational
networks. Clearly, when the radius of influence is increased, in the LETKF
context, the analysis corrections are impacted by spurious correlations. On
the other hand, the quality of the results in the EnKF-MC case is considerably better. When data errors components are uncorrelated ζ can be seen as
a free parameter and the choice can be based on the “optimal performance
of the filter”. For the largest radius of influence ζ = 5 the RMSE values
of the ENKF-MC and the LETKF implementations differ by one order of
magnitude.
Figure 12 reports the RMSE values for the zonal and the meridional
wind component analyses, and for different values of p and ζ. As can be
b can reduce the impact of
seen, the estimation of background errors via B
spurious correlations; the RMSE values of the EnKF-MC analyses remain
small at all assimilation times, from which we infer that the background
error correlations are properly estimated. On the other hand, the impact
of spurious correlations is evident in the context of LETKF. Since most of
the model components are unobserved, the background error correlations
drive the quality of the analysis, and spurious correlations lead to a poor
performance of the filter at many assimilation times.

32

(a) ζ = 1.

(b) ζ = 2.

(c) ζ = 3.

(d) ζ = 4.

(e) ζ = 5.

Figure 10: RMSE of specific humidity analyses with a sparse observational
network (p ∼ 4%) and different values of ζ.

33

Root Mean Square Error (RMSE)

Root Mean Square Error (RMSE)

Specific Humidity (g/Kg)

15

EnKF-MC
LETKF

14
13
12
11
10

1

2

3
Radius (ζ)

4

5

(a) p = 6%.

Specific Humidity (g/Kg)

15

EnKF-MC
LETKF

14
13
12
11
10

1

2

3
Radius (ζ)

4

5

(b) p = 4%.

Figure 11: Analysis RMSE for the specific humidity variable with sparse
observation networks. RMSE values are shown for different values of ζ and
percentage of observed components p.
Figures 13 and 14 provide snapshots of the meridional and the zonal wind
components, respectively, at the first assimilation time. For this particular
case the percentage of observed model components is p = 4%. At this step,
only the initial observation has been assimilated in order to compute the
analysis corrections by the EnKF-MC and the LETKF methods. The background solution contains erroneous waves for the zonal and the meridional
wind components. For instance, for the u model variable, such waves are
clearly present near the poles. After the first assimilation step, the LETKF
analysis solution dissipates the erroneous waves but, the numerical values of
the wind components are slightly greater than those of the reference solutions.
This numerical difference increases at later times due to the highly-nonlinear
dynamics of SPEEDY, as can bee seen in Figure 12. On the other hand, the
EnKF-MC implementation recovers the reference shape, and the analysis
values of the numerical model components are close to that of the reference
solution. This shows again that the use of the modified Cholesky decomposition as the estimator of the background error correlations can mitigate the
impact of spurious error correlations.

34

Variable (units)

ζ
1

2

Zonal Wind Component (u), (m/s)

3

4

5

1

2

Meridional Wind Component (v) (m/s)

3

4

5

p
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%

EnKF-MC
5.514 × 102
6.972 × 102
9.393 × 102
4.187 × 102
6.090 × 102
7.853 × 102
4.388 × 102
6.146 × 102
7.438 × 102
4.323 × 102
5.990 × 102
7.124 × 102
4.456 × 102
6.106 × 102
7.160 × 102
3.540 × 102
5.165 × 102
7.770 × 102
3.009 × 102
4.605 × 102
6.217 × 102
3.172 × 102
4.735 × 102
6.014 × 102
3.399 × 102
4.812 × 102
5.913 × 102
3.626 × 102
5.107 × 102
6.122 × 102

LETKF
5.471 × 102
1.168 × 103
1.737 × 103
1.275 × 103
7.591 × 102
8.569 × 102
1.661 × 103
1.237 × 103
9.997 × 102
1.752 × 103
1.608 × 103
1.258 × 103
1.862 × 103
1.983 × 103
1.602 × 103
4.496 × 102
1.158 × 103
1.749 × 103
7.285 × 102
5.520 × 102
7.420 × 102
9.510 × 102
8.334 × 102
7.455 × 102
1.048 × 103
1.146 × 103
9.026 × 102
1.101 × 103
1.575 × 103
1.102 × 103

Table 2: RMSE values of the wind-components for the EnKF-MC and
LETKF making use of the SPEEDY model.

35

Variable (units)

ζ
1

2

Temperature (K)

3

4

5

1

2

Specific Humidity (g/Kg)

3

4

5

p
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%
12%
6%
4%

EnKF-MC
6.054 × 102
5.692 × 102
6.522 × 102
5.680 × 102
5.193 × 102
5.299 × 102
5.279 × 102
4.982 × 102
4.926 × 102
5.023 × 102
4.757 × 102
4.766 × 102
4.898 × 102
4.644 × 102
4.684 × 102
9.862 × 101
1.133 × 102
1.405 × 102
1.029 × 102
1.146 × 102
1.270 × 102
1.068 × 102
1.205 × 102
1.317 × 102
1.065 × 102
1.246 × 102
1.324 × 102
1.089 × 102
1.301 × 102
1.373 × 102

LETKF
6.033 × 102
6.704 × 102
8.073 × 102
6.693 × 102
5.556 × 102
5.529 × 102
1.217 × 103
6.458 × 102
6.073 × 102
1.817 × 103
1.030 × 103
7.464 × 102
1.600 × 103
1.473 × 103
1.172 × 103
9.026 × 101
1.449 × 102
1.941 × 102
1.125 × 102
1.137 × 102
1.321 × 102
1.341 × 102
1.418 × 102
1.458 × 102
1.640 × 102
1.652 × 102
1.739 × 102
2.078 × 102
1.950 × 102
2.068 × 102

Table 3: RMSE values for the EnKF-MC and LETKF making use of the
SPEEDY model.

36

Zonal Wind Component (U), (m/s)
Background
EnKF-MC
LETKF

4000
3500
3000
2500
2000
1500
1000
500
0

0

5

10
15
Time (days)

20

Meridional Wind Component (V) (m/s)

4500
Root Mean Square Error (RMSE)

Root Mean Square Error (RMSE)

4500

3500
3000
2500
2000
1500
1000
500
0

25

(a) ζ = 3 and p = 12%

Background
EnKF-MC
LETKF

4000

0

5

10
15
Time (days)

20

25

(b) ζ = 4 and p = 12%

Figure 12: RMSE of the LETKF and EnKF-MC implementations for different model variables, radii of influence and observational networks.

37

(a) Reference

(b) Background

(c) EnKF-MC

(d) LETKF

Figure 13: Snapshots of the reference solution, background state, and analysis fields from the EnKF-MC and LETKF for the fifth layer of the meridional
wind component (v).

38

(a) Reference

(b) Background

(c) EnKF-MC

(d) LETKF

Figure 14: Snapshots of the reference solution, background state, and analysis fields from the EnKF-MC and LETKF for the second layer of the zonal
wind component (u).

4.3

Statistics of the ensemble

In this section, we briefly discuss the spread of the ensemble making use
of rank histograms. Of course, we do not claim this to be a verification
procedure but, it provides useful insights about the dispersion of the members
and the level of uncertainty about the ensemble mean. The plots are based on
the 5-th numerical layer of the atmosphere. We collect information across all
model variables and the plots are shown in figures 15, 18, 16, and 17. Based

39

on the results, the proposed implementation seems to be lesser sensitive to
the intrinsic need of inflation than the LETKF formulation. For instance,
after the assimilation, the ensemble members from the EnKF-MC are spread
almost uniformly across different observation times. On the other hand, the
spread in the context of the LETKF is impacted by the constant inflation
factor used during the experiments (1.04) In practice, the inflation factor is
set up according to historical information and/or heuristically with regard to
some properties of the dynamics of the numerical model. This implies that,
the dispersion of the LETKF members after the analysis will rely in how-well
we estimate the optimal inflation factor for such filter. In operational data
assimilation, an answer to this question can be hard to find. We think that
inflation methodologies such as adaptive inflation can lead to better spread
of the ensemble members in the context of the LETKF. For the proposed
method, based on the experimental results, such methodology is not needed.

(a) EnKF-MC

(b) LETKF

Figure 15: Rank-histograms for the Specific Humidity model variable. The
information is collected from the 5-th model layer.

40

(a) EnKF-MC

(b) LETKF

Figure 16: Rank-histograms for the Zonal Wind Component model variable.
The information is collected from the 5-th model layer.

(a) EnKF-MC

(b) LETKF

Figure 17: Rank-histograms for the Meridional Wind Component model variable. The information is collected from the 5-th model layer.

41

(a) EnKF-MC

(b) LETKF

Figure 18: Rank-histograms for the Temperature model variable. The information is collected from the 5-th model layer.

4.4

The impact of SVD truncation threshold

An important question arising from this research is the number of singular
values/vectors to be used in (19). To study this question we use the same
experimental setting and the sparse observational network where only 4%
of the model components are observed. We apply EnKF-MC algorithm and
truncate the summation (19) based on different thresholds σr .
The results are reported in Figure 19. Different thresholds lead to different levels of accuracy for the EnKF-MC analyses. There is no unique value
of σr that provides the best ensemble trajectory in general; for instance, the
best performance at the beginning of the assimilation window is obtained for
σr = 0.05, but, at the end the best solution is obtained with σr = 0.2. This
indicates that the results can be improved when σr is dynamically and optimally chosen. Note that, on average, the results obtained by the EnKF-MC
with σr ∈ {0.15, 0.20, 0.25} are much better than those when σr = 0.10 (and
therefore much better than the results obtained by the LETKF). In Figure
20 snapshots of the specific humidity for different σr are shown. It can be
seen that the spurious errors can be quickly decreased when σr is chosen
accordingly.
In order to understand the optimal truncation level note that the sum-

42

mation (19) can be written as follows:
β [i] =

Nens
X

Z

αj · uj [i] ,

(38)

j=1


1
1
Z T
Z T 
e[i] + θ [i]
· vj [i] · x[i] = · vj [i] · x
τj
τj
1
1
Z T
Z T
e[i] + · vj [i] · θ [i]
· vj [i] · x
=
τ
τ
{z
} |j
{z
}
|j

αj =

Error

Uncorrupted data

e[i] is the perfect data (x[i] = x
e[i] + θ [i] ). The components with small
where x
singular values τj will amplify the error more. The threshold should be large
e[i] , but small enough in order to
enough to include useful information from x
prune out the components with large error amplification. We expect that
model components with large variances will need more basis vectors from
(19) than those with lesser variance. An upper bound for the number of
basis vectors (and therefore the threshold σr ) can be obtained by inspection
of the values αj in (38). Figure 21 shows the weights αj for different singular
values for the 500-th model component of the SPEEDY model. The large
zig-zag behaviors are evidence of error amplifications and therefore, we can
truncate the summation (38) before this pattern starts to take place in the
values of αj .

43

1600

σr=0.05
σ =0.10
r

Root Mean Square Error (RMSE)

Root Mean Square Error (RMSE)

1400

σr=0.15

1200

σr=0.20
σr=0.25

1000

σr=0.30

800

600

400
0

5

10
15
Time (days)

20

(a) Meridional wind component (m/s)

400

20

5

10
15
Time (days)

20

25

σr=0.10

σr=0.30

10
15
Time (days)

600

σ =0.05

σr=0.25

5

800

220

500

300
0

σr=0.30

1000

240

σr=0.20

600

σr=0.25

σr=0.10
σr=0.15

700

σr=0.20

1200

σ =0.05
r

800

r

σr=0.15

(b) Zonal wind component (m/s)
Root Mean Square Error (RMSE)

Root Mean Square Error (RMSE)

900

σ =0.10
1400

400
0

25

σr=0.05

r

σr=0.15

200

σr=0.20
σr=0.25

180

σr=0.30

160
140
120
100
80
0

25

(c) Temperature (K)

5

10
15
Time (days)

20

25

(d) Specific humidity (g/kg)

Figure 19: RMSE for the SPEEDY analyses obtained using different SVD
truncation levels based on the σr values.

44

(a) Reference

(b) σr = 0.05

(c) σr = 0.10

(d) σr = 0.15

(e) σr = 0.20

45

(f) σr = 0.30

Figure 20: Snapshots at the final assimilation time (day 22) of the EnKF-MC
analysis making use of different thresholds σr for ζ = 5 and p = 4%.

Meridional wind component

Zonal wind component

15

20

10

10

j

α

αj

5
0

0
−10

−5
−10
0

20

40
60
Singular vectors j

80

−20
0

100

(a) Meridional wind component (m/s)

20

40
60
Singular vectors j

80

(b) Zonal wind component (m/s)

Temperature

Specific humidity

6

0.4

4

0.2
0

2
αj

αj

−0.2
0

−0.4
−2

−0.6

−4
−6
0

−0.8
20

40
60
Singular vectors j

80

−1
0

100

(c) Temperature (K)

20

40
60
Singular vectors j

80

100

(d) Specific humidity (g/kg)

Figure 21: The effect of θ on the weights αj for some model component i of
the SPEEDY model when ζ = 5 and p = 4%.

5

Conclusions

This paper develops an efficient implementation of the ensemble Kalman filter, named EnKF-MC, that is based on a modified Cholesky decomposition
to estimate the inverse background covariance matrix. This new approach
has several advantages over classical formulations. First, a predefined spar46

sity structure can be built into the factors of the inverse covariance. This
reflects the fact that if two distant model components are uncorrelated then
the corresponding entry in the inverse covariance matrix is zero; the only
nonzero entries in the Cholesky factors correspond to components of the
model that are located in each other’s proximity. Therefore, imposing a
sparsity structure on the inverse background covariance matrix is a form of
covariance localization. Second, the formulation allows for a rigorous theoretical analysis; we prove the convergence of the covariance estimator for a
number of ensemble members that is proportional to the logarithm of the
number of states of the model therefore, when Nens ≈ log n, the background
error correlations can be well-estimated making use of the modified Cholesky
decomposition.
We discuss different implementations of the new EnKF-MC, and asses
their computational effort. We show that domain decomposition can be used
in order to decrease even more the computational effort of the proposed implementation. Numerical experiments are carried out using the Atmospheric
General Circulation Model SPEEDY reveal that the analyses obtained by
EnKF-MC are better than those of the LETKF in the root mean square
sense when sparse observations are used in the analysis. For dense observation grids the EnKF-MC solutions are improved when the radius of influence increases, while the opposite holds true for LETKF analyses. (We
stress the fact that these conclusions are true for our implementation of the
basic LETKF; other implementations may incorporate advances that could
make the filter perform considerably better). The use of modified Cholesky
decomposition can mitigate the impact of spurious correlation during the
assimilation of observations.

Acknowledgements
This work was supported in part by awards NSF CCF–1218454, AFOSR
FA9550–12–1–0293–DEF, and by the Computational Science Laboratory at
Virginia Tech.

47

References
[And12] Jeffrey L. Anderson. Localization and Sampling Error Correction
in Ensemble Kalman Filter Data Assimilation. Monthly Weather
Review, 140(7):2359–2371, 2012.
[BF07] Angela Benedetti and Michael Fisher. Background Error Statistics for Aerosols. Quarterly Journal of the Royal Meteorological
Society, 133(623):391–405, 2007.
[BL08] Peter J. Bickel and Elizaveta Levina. Regularized estimation of
large covariance matrices. Ann. Statist., 36(1):199–227, 02 2008.
[Bue11] Mark Buehner. Evaluation of a Spatial/Spectral Covariance Localization Approach for Atmospheric Data Assimilation. Monthly
Weather Review, 140(2):617–636, 2011.
[CJAS10a] H. Cheng, M. Jardak, M. Alexe, and A. Sandu. A hybrid approach to estimating error covariances in variational data assimilation. Tellus A, 62A(3):288–297, 2010.
[CJAS10b] Haiyan Cheng, Mohamed Jardak, Mihai Alexe, and Adrian
Sandu. A Hybrid Approach to Estimating Error Covariances
in Variational Data Assimilation. Tellus A, 62(3):288–297, 2010.
[CSCC07] E.M. Constantinescu, A. Sandu, T. Chai, and G.R. Carmichael.
Autoregressive models of background errors for chemical data assimilation. Journal of Geophysical Research, 112(D12309), 2007.
[Eve03] Geir Evensen. The Ensemble Kalman Filter: Theoretical
Formulation and Practical Implementation. Ocean Dynamics,
53(4):343–367, 2003.
[GHO99] Gene H. Golub, Per Christian Hansen, and Dianne P. O’Leary.
Tikhonov regularization and total least squares. SIAM Journal
on Matrix Analysis and Applications, 21(1):185–194, 1999.
[GMC+ 06] S. Gillijns, O.B. Mendoza, J. Chandrasekar, B. L R De Moor,
D.S. Bernstein, and A Ridley. What is the Ensemble Kalman
Filter and How Well Does It Work? In American Control Conference, 2006, pages 6 pp.–, June 2006.
48

[Han90] Per Christian Hansen. Truncated singular value decomposition
solutions to discrete ill-posed problems with ill-determined numerical rank. SIAM Journal on Scientific and Statistical Computing, 11(3):503–518, 1990.
[HL86] A. Hollingsworth and P. Lonnberg. The Statistical Structure
of Short-range Forecast Errors as Determined from Radiosonde
Data. Part I: The Wind Field. Tellus A, 38A(2):111–136, 1986.
[Huf91] Sabine Van Huffel. Iterative algorithms for computing the singular subspace of a matrix associated with its smallest singular
values. Linear Algebra and its Applications, 154:675 – 709, 1991.
[JB00] Eric P. Jiang and Michael W. Berry. Solving total least-squares
problems in information retrieval. Linear Algebra and its Applications, 316(1–3):137 – 156, 2000. Special Issue: Conference
celebrating the 60th birthday of Robert J. Plemmons.
[JFW14] Poterjoy Jonathan, Zhang Fuqing, and Yonghui Weng. The Effects of Sampling Errors on the EnKF Assimilation of Inner-Core
Hurricane Observations. Monthly Weather Review, 142(4):1609–
1630, 2014.
[Kep00] Christian L. Keppenne. Data Assimilation into a PrimitiveEquation Model with a Parallel Ensemble Kalman Filter.
Monthly Weather Review, 128(6):1971–1981, 2000.
[KHE07] Misha E. Kilmer, Per Christian Hansen, and Malena I. Español.
A projectionbased approach to generalform tikhonov regularization. SIAM Journal on Scientific Computing, 29(1):315–330,
2007.
[KMB06] Fred Kucharski, Franco Molteni, and Annalisa Bracco. Decadal
interactions between the western tropical pacific and the north
atlantic oscillation. Climate Dynamics, 26(1):79–91, 2006.
[Ler07] P. F. J. Lermusiaux. Adaptive modeling, adaptive data assimilation and adaptive sampling. Physica D Nonlinear Phenomena,
230:172–196, June 2007.

49

[Lor05] Edward N. Lorenz. Designing chaotic models. Journal of the
Atmospheric Sciences, 62(5):1574–1587, 2005.
[LR99] P. F. J. Lermusiaux and A. R. Robinson. Data assimilation via
error subspace statistical estimation. part i: Theory and schemes,
1999. Accessed: 08-29-2015.
[Mai94] B. A. Mair. Tikhonov regularization for finitely and infinitely
smoothing operators. SIAM Journal on Mathematical Analysis,
25(1):135–147, 1994.
[Mol03] F. Molteni. Atmospheric simulations using a gcm with simplified
physical parametrizations. i: model climatology and variability
in multi-decadal experiments. Climate Dynamics, 20(2-3):175–
191, 2003.
[Neu88] Andreas Neubauer. An a posteriori parameter choice for
tikhonov regularization in hilbert scales leading to optimal
convergence rates. SIAM Journal on Numerical Analysis,
25(6):1313–1326, 1988.
[NRSA14] EliasD. Nino Ruiz, Adrian Sandu, and Jeffrey Anderson. An
Efficient Implementation of the Ensemble Kalman Filter Based
on an Iterative Sherman–Morrison Formula. Statistics and Computing, pages 1–17, 2014.
[OHS+ 04a] Edward Ott, Brian R. Hunt, Istvan Szunyogh, Aleksey V. Zimin,
Eric J. Kostelich, Matteo Corazza, Eugenia Kalnay, D. J. Patil,
and James A. Yorke. A Local Ensemble Kalman Filter for Atmospheric Data Assimilation. Tellus A, 56(5):415–428, 2004.
[OHS+ 04b] Edward Ott, Brian R. Hunt, Istvan Szunyogh, Aleksey V. Zimin,
Eric J. Kostelich, Matteo Corazza, Eugenia Kalnay, D. J. Patil,
and James A. Yorke. A local ensemble kalman filter for atmospheric data assimilation. Tellus A, 56(5):415–428, 2004.
[OHS+ 08] Edward Ott, Brian Hunt, Istvan Szunyogh, Aleksey V Zimin,
Eic J. Kostelich, Matteo Corazza, Eugenia Kalnay, D. J. Patil,
and James A. Yorke. A Local Ensemble Transform Kalman Filter
Data Assimilation System for the NCEP Global Model. Tellus
A, 60(1):113–130, 2008.
50

[SO08] Pavel Sakov and Peter R. Oke. A Deterministic Formulation
of the Ensemble Kalman Filter: An Alternative to Ensemble
Square Root Filters. Tellus A, 60(2):361–371, 2008.
[TAB+ 03] Michael K. Tippett, Jeffrey L. Anderson, Craig H. Bishop,
Thomas M. Hamill, and Jeffrey S. Whitaker. Ensemble square
root filters. Monthly Weather Review, 131(7):1485–1490, Jul
2003.
[WT02] Jeffrey S. Whitaker and M. Hamill Thomas. Ensemble Data
Assimilation without Perturbed Observations. Monthly Weather
Review, 16(3):1913–1924, 2002.

51

