Performance analysis of local ensemble Kalman filter
arXiv:1705.10598v2 [math.PR] 21 Mar 2018

Xin T. Tong

∗

March 22, 2018

Abstract
Ensemble Kalman filter (EnKF) is an important data assimilation method for high
dimensional geophysical systems. Efficient implementation of EnKF in practice often
involves the localization technique, which updates each component using only information within a local radius. This paper rigorously analyzes the local EnKF (LEnKF) for
linear systems, and shows that the filter error can be dominated by the ensemble covariance, as long as 1) the sample size exceeds the logarithmic of state dimension and a
constant that depends only on the local radius; 2) the forecast covariance matrix admits
a stable localized structure. In particular, this indicates that with small system and
observation noises, the filter error will be accurate in long time even if the initialization
is not. The analysis also reveals an intrinsic inconsistency caused by the localization
technique, and a stable localized structure is necessary to control this inconsistency.
While this structure is usually taken for granted for the operation of LEnKF, it can also
be rigorously proved for linear systems with sparse local observations and weak local
interactions. These theoretical results are also validated by numerical implementation
of LEnKF on a simple stochastic turbulence in two dynamical regimes.

1

Introduction

Data assimilation is a sequential procedure, in which observations of a dynamical system
are incorporated to improve the forecasts of that system. In many of its most important
geoscience and engineering applications, the main challenge comes from the high dimensionality of the system. For contemporary atmospheric models, the dimension can reach
d ∼ 108 , and the classical particle filter is no longer feasible [1, 2]. The ensemble Kalman
filter (EnKF) was invented by meteorologists [3, 4, 5] to resolve this issue. By sampling the
forecast uncertainty with a small ensemble, and then employing Kalman filter procedures
to the empirical distribution, EnKF can often capture the major uncertainty and produce
accurate predictions. The simplicity and efficiency of EnKF have made it a popular choice
for weather forecasting and oil reservoir management [6, 7].
One fundamental technique employed by EnKF is localization [8, 4, 9, 10, 11]. In most
geophysical applications, each component [X]i of the state variable X holds information of
one spatial location. There is a natural distance d(i, j) between two components. In most
∗

National University of Singapore, mattxin@nus.edu.sg

1

physical systems, the covariance between [X]i and [X]j is formed by information propagation
in space, intuitively its strength decays with the distance d(i, j). In particular, when d(i, j)
exceeds a threshold L, the covariance is approximately zero. This is a special sparse and
localized structure that can be exploited in the EnKF operation. In particular, the forecast
covariance can be artificially enforced as zero if d(i, j) > L. In other words, there is no need
to sample these covariance terms, and indeed sampling from them leads to higher errors [4].
Such modification significantly reduces the sampling difficulty and the associated sample
size. This is crucial for EnKF operation, since often only a few hundred samples can be
generated in practice. Various versions of localized EnKF (LEnKF) are derived based on
this principle, and there is ample numerical evidence showing their performance is robust
against the growth of dimension [4, 9, 10, 11, 12, 13, 14, 15, 16, 17]. Moreover, there is a
growing interest in applying the same technique to the classical particle filters [18, 19, 20].
While there is a consensus on the importance of the localization technique for EnKF,
currently there is no rigorous explanation of its success. This paper contributes to this issue
by showing that in the long run, the LEnKF can reach its estimated performance for linear
systems, if the ensemble size K exceeds DL log d, and the ensemble covariance matrix admits
a stable localized structure of radius L. The constant DL above depends on the radius L
but not on d.
Showing the necessary sampling size has only logarithmic dependence on d is our major
interest. In the simpler scenario of sampling a static covariance matrix, [21] shows that the
necessary sample size scales with DL log d. Generalizing this result to the setting of EnKF
is highly nontrivial, since the target covariance matrix evolves constantly in time, and the
sampling error at one time step has a nonlinear impact on future iterations. By analyzing the
filter forecast error evolution, and compare it with the filter covariance evolution, we show
the filter error covariance can be dominated by the ensemble covariance with high probability.
In other words, the LEnKF can reach its estimated performance.
One important corollary
√
is that if the system and observation noise are of scale ǫ, then the error covariance scales
as ǫ, which indicates that LEnKF can be accurate regardless of the initial condition. Such
property is often termed as accuracy for practical filters or observers [22, 23, 24].
Interestingly, our analysis also captures an intrinsic inconsistency caused by the localization technique. Generally speaking, the localization technique can be applied to the ensemble
covariance matrix, but not the ensemble. However, the Kalman update is applied to the ensemble, but not to the localized ensemble covariance matrix. As these two operations do not
commute, an inconsistency emerges, which we will call the localization inconsistency. This
phenomenon has been mentioned in [9, 25]. Moreover, [15] numerically examines its role
with serial observation processing, and shows that it may lead to significant filter error. In
correspondence to these findings, one crucial step in our analysis is showing that the localization inconsistency is controllable, if the forecast covariance matrix indeed has a localized
structure.
While most applications of LEnKF assume the underlying covariance matrices are localized, rigorous justification of this assumption is sorely missing in the literature. A recent
work [26] considers applying a projection to the continuous time Kalman-Bucy filter, and
shows that if the projection is a small perturbation on the covariance matrix, its impact on
the filter process is also small. It is shown through an example that if the filter system can
be decoupled into independent local parts, a projection similar to the LEnKF localization
2

procedure can be made. Unfortunately, in most practical problems, all spatial dimensions
are coupled with local interactions, and it is very difficult to show that the localization
procedure is a small perturbation.
This paper partially investigates the theoretical gaps mentioned above. We show that
for linear systems with weak local interactions and sparse local observations, the localized
structure is stable for the LEnKF ensemble covariance. Weak local interaction is an intuitive
requirement, else fast information propagation will form strong covariances between far away
locations. Sparse local observation, on the other hand, is assumed to simplify the assimilation
formulas.
In rough words, our main results consist of the following statements.
1. To sample a localized covariance matrix correctly, the necessary sample size scales with
DL log d (Theorem 2.1). This reveals the sampling advantage gained by applying the
localization procedure.
2. While localization improves the sampling, it creates an inconsistency in the assimilation
steps. For the LEnKF ensemble covariance to capture the filter error covariance with
DL log d samples, the localization inconsistency needs to be small (Theorem 2.4).
3. One way to guarantee a small localization inconsistency, is to have a stable localized
structure in the forecast ensemble covariance matrix (Proposition 2.3).
4. The LEnKF forecast covariance has a stable localized structure, if the underlying
linear system has weak interactions and sparse local observations. (Theorem 2.5). So
by points 2 and 3, we know that LEnKF has good forecast skills, since its ensemble
covariance captures the true filter error covariance.
5. The results above scale linearly with the variance of the noises. So when applying
LEnKF to a linear system with small system and observation noises, its long time
performance is accurate (Theorem 2.7).
Section 2 will provide the setup of our problem, and present the precise statements of the
main results. The implication of these results on the issue of localized radius is discussed in
Section 2.6.
Section 3 verifies the theoretical results by implementing LEnKF on a stochastically
forced dissipative advection equation [6]. One stable and one unstable dynamical regimes
are tested. In both of them, LEnKF have shown robust forecast skill with only K = 10
ensemble members, while the dimension varies between 10 and 1000. Moreover the localized
covariance structure and the accuracy with small noises can also be verified for LEnKF in
both regimes.
Section 4 investigates the covariance sampling problem of LEnKF, and proves Theorem
2.1. Section 5 analyzes the localization inconsistency and filter error evolution. It contains
the proofs of Theorem 2.4 and Proposition 2.3. Section 6 studies the localized structure
of linear systems with weak local interactions and sparse observations, and shows that the
small noise scaling can be applied to our results. Section 7 concludes this paper and discusses
some interesting extensions.

3

2
2.1

Main Results
Problem Setup

Since its invention, the ensemble Kalman filter (EnKF) has been modified constantly for two
decades, and its formulation has become rather sophisticated today. In this subsection we
briefly review some of the key modifications, in particular the localization techniques.
The following notations will be used throughout the paper. For two vectors a and b,
kak denotes the l2 norm of a, a ⊗ b denotes the matrix abT . Square bracket with subscripts
indicates a component or entry of an object. So [a]i is the i-th component of vector a. In
particular, we use ei to denote the i-th standard basis vector, i.e. [ei ]j = 1i=j .
Given a matrix A, [A]i,j is the (i, j)-th entry of A. The l2 operator norm is denoted
by kAk
P = inf{c : kAvk ≤ ckvk, ∀v}. The l∞ operator norm is denoted by kAk1 =
maxi j |[A]i,j |. The maximum absolute entry is denoted by kAk∞ = maxi,j |[A]i,j |. We
also use Im to denote the m × m dimensional identity matrix. Given two matrices A and D,
their Schur (Hadamard) product can be defined by entry wise product
[A ◦ D]i,j = [A]i,j [D]i,j .
For two real symmetric matrices A and B, A  B indicates that B−A is positive semidefinite.
Ensemble Kalman Filter
In this paper, we consider a linear system in Rd with partial observations,
Xn+1 = An Xn + bn + ξn , ξn+1 ∼ N (0, Σn ),
Yn+1 = HXn+1 + ζn , ζn+1 ∼ N (0, σo2Iq ).

(2.1)

Throughout our discussion, we assume the matrices An , Σn are bounded:
kAn k ≤ MA ,

mΣ Id  Σn  MΣ Id .

The time-inhomogeneous generality can be used to model intermittent dynamical systems [6,
27]. We assume that the observations are made at q < d distinct locations {o1 , o2 , · · · , oq } ⊂
{1, · · · , d}. This can be modelled by letting
[H]k,j = 1j=ok ,

1 ≤ k ≤ q, 1 ≤ j ≤ d.

(2.2)

Note that the operator norm kHk = 1.
It is well known that the optimal estimate of Xn given historical observations Y1 , . . . , Yn
is provided by the Kalman filter [28], assuming X0 is Gaussian distributed. Unfortunately,
direct implementation of the Kalman filter involves a stepwise computation complexity of
O(d2q). When the state dimension d is high, the Kalman filter is not computationally
feasible.
The ensemble Kalman filter (EnKF) is invented by meteorologists [5] to reduce the computation complexity. K samples of (2.1) are updated using the Kalman filter rules, and their
ensemble mean and covariance are employed to estimate the signal Xn . In specific, suppose
4

(k)

the posterior ensemble for Xn is denoted by {Xn }k=1,...,K . The forecast ensemble of Xn+1
is first generated by propagating the linear system in (2.1):
b (k) = An X (k) + bn + ξ (k) ,
X
n+1
n
n+1

(k)

ξn+1 ∼ N (0, Σn ).

b n+1 , C
bn+1 ), where the mean
The EnKF then estimates Xn+1 with a prior distribution N (X
and covariance are obtained by the forecast ensemble:
K
1 X b (k)
b
X ,
X n+1 =
K k=1 n+1

b (k) := X
b (k) − X
b n+1 ,
∆X
n+1
n+1

K
1 X b (k)
b (k) .
b
∆Xn+1 ⊗ ∆X
Cn+1 =
n+1
K k=1

Applying the Bayes’ formula to the prior distribution and the linear observation Yn+1, a
target Gaussian posterior distribution for Xn+1 can be obtained. There are several ways to
update the forecast ensemble so its statistics approximate the target ones. Here we consider
the standard EnKF in [5, 6] with artificial perturbations:
(k)
(k)
(k)
e n+1 H)X
bn+1
e n+1 Yn+1 − K
b n+1 ζn+1
Xn+1 = (I − K
+K
.

(2.3)

(k)
e n+1 = C
bn+1 H T (σo2 Iq + H C
bn+1H T )−1 . The ζn+1
The Kalman gain matrix is given by K
are
independent noises sampled from N (0, σo2Iq ).
The computation complexity of EnKF is roughly O(K 2 d), assuming An and Σn are sparse
[29]. In practice, the ensemble size K is often less than a few hundred, so the operational
speed is significantly improved. On the other hand, with the sample size K much smaller than
bn+1 often produces spurious correlations
the state space dimension d, the sample covariance C
[30, 5]. Spurious correlations may seriously reduce the filter accuracy, since the Kalman
filter operation hinges heavily on the correctness of covariance estimation. The localization
techniques are often employed to resolve such problems.

Localization techniques
In most geophysical applications, each dimension index i ∈ {1, . . . , d} corresponds to a
spatial location. For simplicity, we assume different indices correspond to different spatial
locations. Let d(i, j) be the spatial distance between the locations i and j specify, then d is
also a distance on the index set {1, . . . , d}. In other words,
• d(i, j) = 0 if and only if i = j;
• d(i, j) = d(j, i);
• d(i, j) + d(j, k) ≥ d(i, k).
For a simple example, one can correspond index i with the integer i, then d(i, j) = |i − j|
clearly defines a distance.
For most geophysical problems that can be modeled by a (stochastic) partial differential
equation, the covariance between two locations is caused by the propagation of information
through local interactions. Information often is also dissipated during its propagation, so
its impact gets less significant when it reaches far-away locations. This leads to a localized
5

covariance structure. In other words, there is a decreasing function φ : [0, ∞) 7→ [0, 1],
φ(0) = 1 such that
[Cn ]i,j ∝ φ(d(i, j)).
In geophysical applications, a localization radius l is often defined, so φ(x) = 0 for x > l.
Consequentially, it is natural to model the localization function as
[Dl ]i,j = φ(d(i, j)).

(2.4)

In particular, the widely used Gaspari-Cohn matrix [31] is of this form with




x
x
φ(x) = 1 +
exp −
1x≤l ,
(2.5)
cl
cl
p
where the radius is often picked with l = 10/3cl or 2cl [32]. Another simple localization
matrix corresponds to the cutoff or heavyside function φ(x) = 1x≤l , and we denote it by
Dlcut . In other words
[Dlcut ]i,j = 1d(i,j)≤l .
(2.6)
As a remark, while (2.5) is more useful in practice, (2.6) is much simpler for theoretical
analysis and interpretation. Most of our analysis results in below only apply to (2.6), except Theorem 2.1. It will be very interesting to generalize the analysis framework here for
localization functions like (2.5).
The notion of localization radius is closely related to the bandwidth of a matrix [33]. For
a matrix A, we define its bandwidth as:
l := inf{x ≥ 0 : [A]i,j = 0 if d(i, j) > x}.

(2.7)

The bandwidth roughly captures how fast different components interact with each other. If
A has bandwidth l, each component interacts with at most Bl components when product
with A, where the volume constant Bl is defined by
Bl = max #{j : d(i, j) ≤ l}.
i

(2.8)

A localized covariance structure is extremely useful for EnKF. It indicates only covariances between nearby indices are worth sampling. By ignoring the far apart covariances,
the necessary sampling size can be significantly reduced. To apply this idea, the localization
technique modifies the Kalman gain matrix in (2.3), and ensures the assimilation updates
from far away observation is insignificant. There are two main types of localization methods
in the literature, domain localization and covariance localization [14]. This paper discusses
only the former, while similar analysis should in principal applies to the latter as well.
With domain localization, the i-th component is updated using only observations of
indices within distance l, which are elements of Ii = {j : d(i, j) ≤ l}. Let PIi be the
projection matrix of a Rd vector to its components on Ii , note that it is diagonal so it is
bi
b
symmetric. Then C
n+1 := PIi Cn+1 PIi contains the local covariance relevant to the i-th
component. The corresponding Kalman gain is
i
bi H T (σ 2 Iq + H C
b i H T )−1 ,
Kn+1
=C
n+1
o
n+1

6

(2.9)

i
and the i-th component is updated using the i-th row of (2.9), namely ei eTi Kn+1
. Again
ei is the i-th standard basis vector of Rd . The final Kalman gain matrix patches all rows
together
d
X
i
b
Kn+1 =
ei eTi Kn+1
.
(2.10)
i=1

i
b n+1 H is of bandwidth l as
Since each Kn+1
has nonzero entries only with indies in Ii × Ii , K
well. The proof in Proposition 2.3 below verifies this. Therefore, each component is updated
using observations of distance at most l from it.

Localized EnKF with covariance inflation
Other than spurious correlations, a small sampling size also jeopardizes the EnKF operation,
as the forecast covariance is often undervalued [34, 35, 23]. In order to resolve this issue, the
covariance needs to be inflated with a fixed ratio r > 1. [23] has shown these modification
are pivotal to EnKF performance. We also incorporate this idea in our LEnKF.
(k)
In summary, the localized EnKF (LEnKF) updates an posterior ensemble {Xn , k =
P
(k)
(k)
(k)
1, · · · , K} of its mean X n = K1 K
and spread ∆Xn = Xn − X n through the
k=1 Xn
b n+1 given by (2.9) and (2.10):
following steps with K
b n+1 = An X n + bn ,
X

bn+1 =
∆X
(k)

√

K
1 X b (k)
(k)
bn+1
b
∆Xn+1 ⊗ ∆X
,
Cn+1 =
K
k=1
(k)
∆Xn+1 =

(k)

r(An ∆Xn(k) + ξn+1),

(k)

ξn+1 ∼ N (0, Σn ),

b n+1 H)X
b n+1 + K
b n+1 Yn+1 ,
X n+1 = (I − K

b n+1 H)∆X
b (k) + K
b n+1 ζ (k) ,
(I − K
n+1
n+1

(2.11)

(k)

ζn+1 ∼ N (0, σo2 Iq ).

The posterior covariance matrix can be obtained through the spread
Cn+1

K
1 X
(k)
(k)
=
∆Xn+1 ⊗ ∆Xn+1 .
K k=1

Note here we update the mean and ensemble spread, the ∆ terms,
is different
P (k) separately.
P This
(k)
1
1
ξn+1 and K
ζn+1 are ignored
from the standard EnKF, since the average noise terms K
P
(k)
for simplicity. Also the sum of the ensemble spread,
∆Xn , may not be zero. On the
other hand, these differences are small by the law of large numbers. The proofs can also be
generalized to admit these terms, but the discussion will be notationally complicated.
One classical property of the Kalman filter is that the filter covariances and the Kalman
gain matrices are predetermined with no dependence on the realization of system (2.1). This
is inherited by the LEnKF (2.11), the covariances and Kalman gain depend only on the
(k) (k)
sample noise ξn , ζn realizations, but not on (Xn , Yn ).
To illustrate, consider the filtration generated by sample noise realization,
b , ξt , ζ , t = 1, . . . , n, k = 1, . . . , K}.
FnS = σ{∆X
0
t−1
(k)

(k)

(k)

7

(2.12)

Using induction, it is easy to verify the ensemble spread, ensemble covariance and Kalman
gain, are all FnS adapted:
b (k) , ∆X (k) , C
bn , Cn−1 , K
bn ∈ F S.
∆X
n
n−1
n

W
S
The corresponding conditional expectation is denoted by EFnS . We will use F∞
= FnS to
denote the σ-field for all ensemble spread information.
The other randomness of EnKF comes from the realization of system (2.1). We can
S
average out this part of randomness by conditioning on F∞
, which we will denote as ES .
This is useful when comparing the filter error and sample covariance. The natural filtration
generated by all random outcome at time n is
(k)
b0(k) , ξt , ζt−1 , ξt(k), ζt−1
Fn = σ{X0 , X
, t = 1, . . . , n, k = 1, · · · , K}.

We will denote the conditional expectation with Fn as En .

2.2

Sampling errors of localized forecast covariance

Since EnKF relies on the ensemble forecast covariance matrix to assimilate new observations,
its performance depends on the accuracy of the sampling procedure. The sampling procedure
updates the forecast matrix from time n to n + 1.
bn , based on the Kalman update rule, the inflated
Given the forecast ensemble covariance C
bn ), with the posterior Riccati map
target forecast covariance at n + 1 is given by rRn (C
bn ) := An (I − K
b n H)C
bn (I − K
b n H)T AT + σ 2 An K
b nK
b T AT + Σn .
Rn (C
n
o
n n

(2.13)

P
bn+1 = 1
b (k) ⊗ ∆X
b (k) is generated by the
The real ensemble forecast covariance C
∆X
n+1
n+1
K
ensemble spread
√
√
√
b (k) = rAn (I − K
b n H)∆X
b (k) + rAn K
b n ζ (k) + rξ (k) .
∆X
(2.14)
n+1
n
n
n+1

bn+1 over ζn(k) and ξ (k) matches Rn (C
bn ), that
It is straight forward to verify the average of C
n+1
bn+1 = Rn (C
bn ).
is, En C
bn+1 − rRn (C
bn )k, it is necessary to have a
In order to control the sampling error kC
sufficiently large K. Unfortunately, the size of K would need to grow linearly with d [21].
(k)
(k)
bn = K
b n = 0, Σn = Id , r = 1, then ∆X
bn+1
As a simple example, let C
= ξn+1 are i.i.d.
p
bn+1 k = 1 + d/K with
samples from N (0, Id), and the target sample matrix is Id . Yet kC
high probability by the Bai-Yin’s law [36]. In practical settings, K ≪ d, so the sample
covariance is unlikely to be correct.
As discussed in Section 2.1, the main idea of localization is that we assume the target
bn ) is localized, so it suffices to consider Rn (C
bn ) ◦ DL, which can be sampled
covariance Rn (C
bn+1 ◦ DL . Here DL can be any matrix of form (2.4), where its radius L does not
by C
need to match l used in (2.9). In fact, we will mostly use DL = DLcut (2.6) with L ≥ 4l
in our discussion. One important advantage gained by localization is that, in order for
bn+1 − Rn (C
bn )) ◦ DL k to be small, the
the covariance sampling to be accurate, that is k(C
necessary sample size scales only with DL log d, instead of d, where DL is some constant
8

that only depends on L. This phenomenon was discovered in statistics [21], assuming the
samples are generated from one fixed distribution.
But in EnKF, the conditional mean of
√
(k)
b
b n H)∆X
bn(k) . A generalization of [21]
each sample is different, i.e. En ∆Xn+1 = rAn (I − K
is our first result:

Theorem 2.1. For any fixed group of ak ∈ Rd , k = 1, . . . , K, and K i.i.d. samples zk ∼
N (0, Σz ). Consider the sample covariances
Z=

K
1 X
(ak + zk ) ⊗ (ak + zk ),
K k=1

Let

Σa =

1/2

K
1 X
ak ⊗ ak .
K k=1

1/2

σa,z = max{[Σz ]i,i , [Σa ]i,i [Σz ]j,j }.
i,j

Z concentrates around its mean in the following two ways, where c is an absolute constant:
a) Schur product with a symmetric matrix DL . For any t ≥ 0


P(k(Z − EZ) ◦ DL k ≥ kDL k1 σa,z t) ≤ 8 exp 2 log d − cK min{t, t2 } .

Recall that kDL k1 := maxi

Pd

j=1 |[DL ]i,j |,

which is often independent of d.

b) Entry-wise. Consider kZ − EZk∞ = maxi,j |[Z − EZ]i,j |, then for any t ≥ 0

P(kZ − EZk∞ ≥ σa,z t) ≤ 8 exp 2 log d − cK min{t, t2 } .
In application to LEnKF, we will let
√
b n H)∆X
bn(k) ,
ak = rAn (I − K

zk =

√

b n ζn(k) +
rAn K

√

(k)

rξn+1 ,

bn+1 ◦ DL concentrates around rRn (C
bn ) ◦ DL . The exact
and Theorem 2.1 shows that C
statement is given below by Corollary 5.4. The result in [21] is equivalent to the special case
where ak ≡ 0. Fortunately, the generalization is not difficult and is in Section 4.

2.3

Localization inconsistency with localized covariance

While the localization technique makes the covariance sampling much easier, they also introduce additional errors. The fundamental reason is that the localization techniques are
applied to the covariance matrices, but cannot be applied to the ensemble members themselves. On the other hand, the analysis update is applied to the ensemble but not to the
covariance. This leads to a matrix inconsistency [9, 25, 15].
b n −Xn . At this moment,
To illustrate, we look at the forecast filter error at time n, ên = X
the sample noise realization of FnS is available, so it is natural to consider the conditional
covariance of the forecast filter error :
EFnS ên ⊗ ên = ES ên ⊗ ên .
9

The identity holds because the sample noises after time n are independent of ên ∈ Fn .
Suppose this covariance is captured by the localized ensemble covariance, in other words
bn ◦ DL . Based on the LEnKF formulation (2.11), the filter errors after the
ES ên ⊗ ên = C
next assimilation step and forecast step are:
bn − K
b n (H X
b n − HXn − ζn ) − Xn = (I − K
b n H)ên + K
b n ζn ,
en = X n − Xn = X

b n+1 − Xn+1 = An (X n − Xn ) − ξn+1 = An (I − K
b n H)ên + An K
b n ζn − ξn+1. (2.15)
ên+1 = X

b n ∈ F S , ζn and ξn+1 are independent of F S , the new forecast error
Since the Kalman gain K
n
∞
covariance is
b n H)(ES ên ⊗ ên )(I − K
b n H)T + σ 2 K
b bT T
ES ên+1 ⊗ ên+1 = An [(I − K
o n Kn ]An + Σn ,
b n H)[C
bn ◦ DL ](I − K
b n H)T + σ 2 K
bnK
b T ]AT + Σn =: R′ (C
bn ).
= An [(I − K
o

n

n

n

(2.16)

On the other hand, the ensemble covariance is generated by the update in (2.14). With no
bn+1 ◦ DL is near its average
inflation, r = 1, Theorem 2.1 indicates C
bn ) ◦ DL = [An [(I − K
b n H)C
bn (I − K
b n H)T + σ 2 K
b bT T
Rn (C
o n Kn ]An + Σn ] ◦ DL .

(2.17)

bn ) is defined by (2.13).
Recall the posterior Riccati map Rn (C
The difference between (2.16) and (2.17) can be interpreted as the inconsistency caused
by commuting the localization and Kalman covariance update. In order for the ensemble
covariance to capture the error covariance, it is necessary for this difference to be small. This
is an issue not governed by the sampling scheme, but governed by the localization operation.
As discussed in the introduction, the major motivation behind localization techniques is
that the covariance is localized. We formalize this notion through the following definition.
Definition 2.2. Given a decreasing function Φ : R+ 7→ [0, 1] with Φ(0) = 1, we say the
bn follows an (Mn , Φ, L)-localized structure, if
forecast covariance sequence C
(
d(i, j) ≤ L;
bn ]i,j | ≤ Mn Φ(d(i, j))
|[C
(2.18)
Mn Φ(L)
d(i, j) > L.
The decay function Φ and L need not coincide with the φ and l used in Kalman gain
localization (2.4). This flexibility is useful when we try to verify the localized structure.
Intuitively, in order for localization techniques to be effective, we need Φ(x) to be near zero
when x is large. This holds true for most localized covariance structures, such as the Gaspari
Cohn matrix (2.5), and also the function Φ(x) = λxA with a certain λA < 1, which will appear
below in Theorem 2.5 for linear systems.
One interesting phenomenon, is that if the forecast covariance is already localized, then
the localization inconsistency is in general small:
bn
Proposition 2.3. Suppose kAn k ≤ MA , An and Σn are of bandwidth less than l, and C
L
follows an (Mn , Φ, L)-localized structure, then the localization inconsistency with DL = Dcut
and L ≥ 4l, given by
∆loc = (2.16) − (2.17),
10

has nonzero entries only around the localization boundary:
[∆loc ]i,j = 0

if |d(i, j) − L| > 2l.

Moreover, it is bounded by
k∆loc k ≤ Mn MA2 (1 + σo−2 Bl Mn )2 Bl2 BL,l Φ(L − 2l).

(2.19)

BL,l is a volume constant BL,l = maxi #{j : |d(i, j) − L| ≤ 2l}, and Bl is given by (2.8).
Note that if Φ(L − 2l) is close to zero, the right side is very small.

2.4

Main result: LEnKF performance

There are different ways to quantify the performance of EnKF. One approach is to compare
EnKF with its large ensemble limit, which is the Kalman filter, and estimate the convergence
rate [?, 37, 38, 39]. Moreover, advanced sampling techniques, such as multilevel Monte Carlo,
can be applied to the EnKF procedures, and speed up the convergence [?, ?]. However, these
results have not investigated the dependence of sample size K on the underlying dimension,
thus they are not helpful in explaining the advantages of the localization procedures. Moreover, the large ensemble limit for LEnKF is not necessarily the optimal, since the localization
techniques may violate the Bayes’ formula.
A more practical approach looks for qualitative EnKF properties, where the necessary
sample size K scales with quantities much less than d [40, 41, 42, 43], for example a low
effective dimension [23]. One central issue of EnKF is that, unlike Kalman filter, it estimates
the forecast uncertainty by the ensemble covariance, which can be faulty. Since the forecast
covariance matrix plays a pivotal role in the EnKF operation, it is important to ask if the
ensemble covariance captures the real filter error covariance.
In our particular case, we are interested in finding a bound for filter error covariance
bn . Note that the
ES ên ⊗ ên . We will compare it with the filter ensemble covariance C
S
conditioning ES is with respect to the sample noise filtration F∞
given in (2.12), moreover
S
b
note that Cn ∈ F∞ . Therefore the comparison is legitimate. By showing ES ên ⊗ ên is
bn with large probability, we demonstrate that the LEnKF
dominated by a proper inflation of C
reaches its estimated performance. In order to achieve that, we need the localized structure
to be stable as well.
Theorem 2.4. Suppose the forecast ensemble covariance follows a stable (Mn , Φ, L)-localized
structure, and the sample size K exceeds DL log d with a constant DL that depends on L, the
LEnKF (2.11) reaches its estimated performance in the long time average. In specific, for
any δ > 0, suppose the following conditions hold
1) In the signal-observation system (2.1), An and Σn are of bandwidth l, moreover
kAn k ≤ MA ,

mΣ Id  Σn  MΣ Id ,

MA2 ≥ mΣ .

b0 + ρId ) for some r0 and ρ that
2) Suppose the initial error satisfies ES ê0 ⊗ ê0  r0 (C
0 < r0 ,

0 < ρ < ( 21 −

1
) min{MA2 /mΣ , σo2 }.
2r

This can always be achieved by picking a larger r0 .
11

bn follows a (Mn , Φ, L)-localized structure as in Definition 2.2.
3) The forecast covariance C
Moreover, the localized structure is stable, so there are constants B0 , D0 and M0 so that
T
1 X
1
b0 k + D0 ) + M0 .
Mn ≤ (B0 EkC
E
T n=1
T

(2.20)

4) The localized structure Φ and radius L satisfy
L ≥ 4l,

−1
Φ(L − 2l) ≤ δ 3 BL,l
MA−2 Bl−6 .

The volume constants are given by Proposition 2.3.
5) The sample size K > Γ(rBl δ −1 , d), with
2

, 18x
log d},
Γ(x, d) = max{9x2 , 24x
c
c

(2.21)

and the absolute constant c is given by Theorem 2.1.
Then for any 1 < r∗ < r, the filter error covariance is dominated by the filter covariance
bn ◦ DL + ρId )
ES ên ⊗ ên  r∗ (C
cut

with high 1 − O(δ) probability in long time average
1−

T −1
1X
bn ◦ DL + ρId ))
P(ES ên ⊗ ên  r∗ (C
cut
T n=0

≤

2.5

b0 k + D0 )
δ(B0 kC
r0
2r 1/3
+
(ρ−1 Bl2 MA2 + (ρσ
1/3 )
o)
T log r∗
T log r∗

δ  −1 2 2
−1
2r 1/3
r 1/3 2/3
+
(ρ Bl MA + (ρσo )1/3 )M0 + ρ MΣ + 2 ρ1/3 σo
.
log r∗

Weak local interaction with sparse observations

By Theorem 2.4, the stability of localized structure is a necessary condition for the LEnKF
to reach its estimated performance. While in practice this condition is often assumed to
be true to motivate the localization technique, and one can check it while the algorithm is
running, it is interesting to find some sufficient a-priori conditions of system (2.1), so that
(2.20) holds. Unfortunately, rigorous investigations in this direction is sorely missing. Here
we provide a stability analysis in a simple setting.
The origin of localized covariances is intuitively clear. In most physical systems, the
covariance between [X]i and [X]j comes from information propagation in space. So if the
propagation is weak and decays at the same time, there will be a localized covariance. For
our linear models, the information propagation is carried by local interactions, described by
the off diagonal terms of An . To enforce its weakness, we assume that there is a λA < 1,
such that
)
( d
X
−d(i,k)
|[An ]i,k |λA
≤ λA .
(2.22)
max
i

k=1

12

For the simplicity of our discussion, we also assume the system noise is diagonal Σn = σξ2 Id .
−d(i,k)
Note that λA < 1, so λA
is a large number when i and k are fart apart. So condition
(2.22) constraints the long distance interaction, measured by |[An ]i,k |, to be weak. In other
words, (2.22) models a local interaction. If we concern the unfilter covariance of the sequence
[X]i , then λA < 1 is sufficient to guarantee the covariance is localized, using Proposition 6.2
in below.
The main difficulty actually comes from the observation part. For simplicity, we require
the observations in (2.2) to be sparse in the sense that d(oi , oj ) > 2l for any i 6= j. Recall that
oi is the i-th observable component. Then for each location i ∈ {1, · · · , d}, there is at most
one location o(i) ∈ {o1 , · · · , oq } such that d(i, o(i)) ≤ l. This will significantly simplify the
analysis step and yield an explicit expression. Sparse observations are in fact quite common
in practice. Moreover, it is also possible to generalize the results here to non-sparse scenario,
by using sequential assimilation [15]. But the conditions will be much more involved.
Under the sparse observation scenario, the following function describes how does the
bn update to the one of C
bn+1 :
localized structure of C
o
n
2
(2.23)
ψλA (M, δ) = (r + δ) max λA M 1 + σo−2 M + λA σo−2 M 2 , λ2A M + σξ2 .
This function provides a way to ensure stable localized structure:

Theorem 2.5. Given a LEnKF (2.11), suppose the following holds
1) The system noise is diagonal and the observations are sparse
Σn = σξ2 Id ,

d(oi , oj ) > 2l, ∀i 6= j.

2) There is a λA < r −1 such that (2.22) holds.
3) There are constants
0 < δ∗ < min{0.25, 21 (λ−1
A − r)},

M∗ ≥

such that ψλA (M∗ , δ∗ ) ≤ M∗ with ψλA given by (2.23).

(r + 2δ∗ )σξ2
,
1 − λA

−1

4δ∗
4) Denote n∗ = 2L + ⌈ log
⌉. The sample size K exceeds
log λ−1
A


1
2
−2
−1
K > max − 2 2L log(16d n∗ δ∗ ), Γ(2rδ∗ , d) .
cδ∗ λA

(2.24)

Then the forecast ensemble covariance follows a stable localized structure (Mn , Φ, L) with
Φ(x) = λxA . In specific, the stochastic sequence Mn is dissipative every n∗ steps:
1
E0 Mn∗ ≤ M0 + (1 + 2δ∗ )M∗ .
2
The long time average condition (2.20) can be verified by
T
2n∗
1X
b0 k + M∗ ) + 2(1 + δ∗ )M∗ .
EMk ≤
(EkC
T k=1
T λLA

13

Remark 2.6. Note that
ψλA (M, 0) = max{rλA M(1 + σo−2 M)2 + rλA σo−2 M 2 , rλ2A M + rσξ2 },

With sufficiently small λA or σo−1 , ψλA (M, δ) < M can have a solution, so condition 3) holds.

2.6

Localization radius

One important and difficult issue of LEnKF implementation is how to choose the localization
radius l. The theoretical results above shed some light over this issue qualitatively. It is
worth noticing that this paper has two localization radii. l is the one used for LEnKF(2.11)
formulation, and L is used for the filter error theoretical analysis. But generally speaking L
and l should be picked so that L ≥ 4l, so we concern only of L in the following. We also
assume that Φ(x) = λxA from Theorem 2.5 for simpler discussion.
A smaller localization radius simplify the sampling task by focusing on a smaller assimilation domain, and significantly reduces the necessary sample size. This comes from
two perspectives. First, in order for the LEnKF to sample the correct localized covariance
matrix, condition 5) of Theorem 2.4 requires the sample size to grow polynomially with L,
since kΦk1 is summing over BL entries. Second, the localized covariance structure can be
very delicate at the boundary, and to maintain it one needs the random forecast covariance
to have sampling error of scale λLA . This leads to the exponential dependence of K on L, as
in condition 4) of Theorem 2.5.
On the other hand, a larger localization radius L reduces the size of the localization
inconsistency. Based on Proposition 2.3, the localization inconsistency is of order Φ(L−2l) =
L−2l
λA
, because within inequality (2.19), Bl is independent of L, and BL,l is also independent
of L if i, j are taken from {1, · · · , d}. This becomes condition 4) of Theorem 2.4, where we
need the localization radius to be large, so the inconsistency is bounded by the tolerance.

2.7

LEnKF accuracy with small noises

In practice, with frequent and accurate observations, the system noises, Σn and σo2 , are
often of scale ǫ. In this scenario, the LEnKF has its error covariance scale with ǫ in long
time, showing an accurate forecast skill. Moreover, there is no requirement that the initial
ensemble to have error of scale ǫ, meaning the LEnKF can converge to the signal Xn given
enough time.
Theorem 2.7. Suppose, the signal-observation system (2.1) satisfies the conditions of Theorem 2.5, and its LEnKF is tuned to satisfy the conditions of Theorem 2.4 except (2.20).
Then if the same LEnKF is applied to the following system
ǫ
Xn+1
= An Xnǫ + bn + ξn ,
ǫ
ǫ
Yn+1
= HXn+1
+ ζn ,

ξn+1 ∼ N (0, ǫσξ2 Id ),

ζn+1 ∼ N (0, ǫσo2Iq ),

it has small filter error covariance of scale ǫ. In particular, the ensemble covariance is of
scale ǫ in long time average
T
1X
bn k∞ ≤ 2n∗ (EkC
b0 k + ǫM∗ ) + 2(1 + δ∗ )ǫM∗ .
EkC
T n=1
T λLA

14

bn with high probability:
Moreover, the real filter covariance is dominated by C
T −1
1X
bn ◦ DL + ǫρId ))
P(ES ên ⊗ ên  r∗ (C
1−
cut
T n=0

≤

b0 k + ǫM∗ )
2δn∗ (EkC
r0
2r 1/3
+
(ρ−1 Bl2 MA2 + (ρσ
1/3 )
L
o)
T ǫ log r∗
T ǫλA log r∗

δ  −1 2 2
−1 2
2r 1/3
r 1/3 2/3
2(ρ Bl MA + (ρσ
.
)(1
+
δ
)M
+
ρ
σ
+
2
σ
+
∗
∗
1/3
ξ
ρ1/3 o
o)
log r∗

Note that ǫ appears only in terms that converge to zero with T → ∞.

Remark 2.8. We need the system to follow the conditions in Theorem 2.5 only to ensure
the stable localized structure exists. If one can find other conditions to verify that the LEnKF
follows an (Mn , Φ, L) localized structure such that Mn converges to a scale of ǫ, the conditions
in Theorem 2.5 can be replaced.

3

Numerical experiments

There is plenty of numerical evidence showing that LEnKF has good forecast skill even with
nonlinear dynamical systems. Moreover, this paper intends to understand LEnKF from a
theoretical perspective, not an empirical one. On the other hand, several new concepts and
conditions are introduced in our analysis framework. To understand their significance, we
conduct a few simple numerical experiments in this section.

3.1

Experiments setup: a stochastic turbulence model

We consider a stochastically forced dissipative advection equation on an one dimensional
periodic domain from Section 6.3 of [6]:
∂u(x, t)
∂u(x, t)
∂ 2 u(x, t)
=c
− νu(x, t) + µ
+ σx Ẇ (x, t).
∂t
∂x
∂x2
To transform it to a discrete linear system, we apply the centered difference formula with
spatial grid size h, and Euler scheme with time step ∆t. We assume W (x, t) is a white noise
in both time and space. The discretized signal-system [Xn,1 , · · · , Xn,d]T follows
√
Xn+1,i = a− Xn,i−1 + a0 Xn,i + a+ Xn,i+1 + σx ∆tWn+1,i , i = 1, . . . , d;
(3.1)
− c∆t
, a0 = 1 − 2µ∆t
− ν∆t, a+ = µ∆t
+ c∆t
.
a− = µ∆t
h2
2h
h2
h2
2h
The indices should be interpreted cyclically, that is Xn,0 = Xn,d and Xn,d+1 = Xn,1 . The
natural distance between indices is d(i, j) = min{|i − j|, ||i − j| − d|}. The system noises
Wn,i are independent samples from N (0, 1). We also initialize X0,i ∼ N (0, 1) for simplicity.
Evidently, if we formulate (3.1) in the format of (2.1), the corresponding matrix An is
constant with bandwidth l = 1. In other words it is tridiagonal. We assume one observation
is made every p components with independent Gaussian noise Bn,k ∼ N (0, 1):
Yn,k = Xn,p(k−1)+1 + σo Bn,k .
15

A simple LEnKF with domain localization radius l = 1, inflation r = 1.1 will be applied to
recover Xn . A small sample size K = 10 is taken. As comparison, we implement a standard
EnKF with the same inflation, sample size and sample noise realization. A standard Kalman
filter is also computed to indicate the optimal filter error. We are interested to see
• Does LEnKF have a close to optimal performance? Does localization play a key role?
• Is filter performance robust against dimension increase?
• Does filter performance scale with the noise strength?
• Does the LEnKF ensemble covariance localize, and is this structure stable?
• Do the a-priori conditions of Theorem 2.5 hold?
In the discussion below, we consider dimension in a wide range d = 10, 100, 1000. Yet
we will fix the grid size h in each regime. This corresponds to a sequence of domains with
increasing size, but not a fixed domain with increasing refinement. Although the latter can
also have very high dimension, localization is not a suitable tool; a proper projection to the
low effective dimension should be more effective [23]. Also it is worth noticing that there are
better ways to filter (3.1), such as Fourier domain filtering [6]. We are running LEnKF here
just to support our theoretical analysis.

3.2

Regime I: strong dissipation

We first consider a regime of (3.1) with strong uniform damping and weak advection
h = 1,

∆t = 0.1,

p = 5,

ν = 5,

c = 0.1,

µ = 0.1,

σx = σo = 1.

In this regime, the conditions of Theorem 2.5 can be verified. In particular, (2.22) can
be formulated as
−1
a− λ−1
(3.2)
A + a0 + a+ λA ≤ λA .
Direct numerical computation shows that λA = 0.5186 satisfies this relation. Furthermore,
we can verify that (δ ∗ , M∗ ) = (0.128, 0.2187) satisfy condition 3) of Theorem 2.5. Theobn follows localized structure
rem 2.5 predicts a stable stochastic sequence Mn exists so C
(Mn , Φ, 4), where Φ(x) = λx∧L
and Mn has its mean bounded by 8.8959. On the other
A
hand, Theorem 2.5 requires the sample size to be around K = 2.8 × 104 for d = 100, and
K = 7.34 × 104 for d = 106 . We will see K = 10 is sufficient for LEnKF to perform well
numerically. The overestimate is reasonable as theoretical analysis is often too conservative.
The main point of theoretical analysis is showing a logarithmic dependence of K on the
dimension.
The numerical results are presented in Figure 3.1. In subplot a) the dimension average
square forecast error
b n |2 /d
DSE := |Xn − X
of LEnKF is plotted for 100 iterations. The time mean DSE (MSE) is around 0.142 for
d = 100. This is comparable with the optimal Kalman filter MSE 0.129. Moreover, this
16

DSE

a) LEnKF

b) EnKF

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0

0

10

20

30

40

50

60

70

80

90

d=10
d=100
d=1000
optimal

100

t
c) Stability of localization structure

0.1

0

0

10

20

30

40

50

60

70

80

90

100

t
d) LEnKF MSE with small noise

6
5
−1

10
4
3

−2

10

2
1
0

−3

0

10

20

30

40

50

60

70

80

90

10

100

−1

10

Noise strength ε

t

Figure 3.1: Filter performance in stable regime I.
performance is robust for all dimensions, MSE=0.137 for d = 10 and MSE=0.143 for d =
1000, while the oscillation is stronger in d = 10 case due to averaging over a small dimension.
Since this regime is very stable, EnKF without localization also has surprisingly good
performance, as shown in subplot b). Its MSE is around 0.15, which is worse than LEnKF.
This shows that, while the conditions of Theorem 2.5 are sufficient for LEnKF to work well,
they might be too strong. It will be interesting if sharper working conditions for LEnKF
can be found. It will also be interesting if one can show such strong conditions can already
guarantee EnKF to work without localization.
Two other properties predicted by our theory are also validated. In subplot c), the
localization status Mn is plotted for all three dimensions. All three time sequences are
stable, and they are all bounded below the theoretical estimate
8.8959
√
√ from Theorem 2.5.
ǫ
ǫ
We also test LEnKF with small scale system noises σx = ǫσx , σo = ǫσo . In subplot d),
1
we plot the time mean DSE of ǫ = 1, 12 , 41 , . . . , 32
in logarithmic scales. It is clear that the
LEnKF has the correct MSE scale of ǫ as Theorem 2.7 predicted.

3.3

Regime II: strong advection

The second regime we considered has a strong advection, while the damping is weak:
h = 0.2,

∆t = 0.1,

p = 5,

ν = 0.1,

c = 2,

µ = 0.1,

σx = σo = 1.

This regime is close to unstable, since the linear system map An has spectral norm 0.99.
(3.2) does not have a solution below 1, so the conditions of Theorem 2.5 are not verifiable.
Nevertheless, we find empirically the LEnKF ensemble covariance matrices are localized. In
Figure 3.2, we demonstrate this by plotting
!
d−x
d
X
X
1
bn ]i,i+x−d |
bn ]i,i+x | +
b
|[C
|[C
Φ(x)
= E
d
i=1
i=d−x+1
17

0

10

b) EnKF

a) LEnKF
4
10

10

DSE

3
5

10
2

0

10

1

0

0

10

20

30

40

50

60

70

80

90

d=10
d=100
d=1000
optimal

100

t
c) Largest covariance component

−5

10

0

10

20

30

40

50

60

70

80

90

100

t
d) LEnKF MSE with small noise

10
8
0

10
6
4

−1

10
2
0

0

10

20

30

40

50

60

70

80

90

−1

100

0

10

10

Noise strength ε

t

Figure 3.3: Filter performance in stable regime II.
using empirical average from 1000 samples with d = 100, n = 100. The clear covariance
strength transition around x = 4 indicates that the ensemble covariance is localized. Therefore Theorem 2.4 applies and predicts that LEnKF will have a good performance.
This is indeed the case. In subplot a) of
Figure 3.3, we see that LEnKF has a forecast
skill. The MSE is around 1.63 for d = 100,
where the optimal Kalman filter MSE is
1.06. This performance does not change
much with the dimension, MSE=1.42 for
d = 10, MSE=1.72 for d = 1000. The EnKF
on the other hand is highly unstable except
Distance
for the low dimension d = 10 case. In subplot b), we see for d = 100 and 1000, the
Figure 3.2: Localization structure
DSE of EnKF grows exponentially to 1010 .
This is a phenomenon known as EnKF catastrophic filter divergence, previously studied by [6, 43]. Now this also demonstrates how
important is the localization technique. Such divergence can be resolved by introducing an
adaptive additive inflation, where the stability can be rigorously proved [42].
In this unstable regime, LEnKF retains its stability and accuracy. Since the localization
structure does not have a theoretical ground in this regime, Figure subplot c) plots only
bn . From it we see the LEnKF ensemble covariance is
the largest matrix component of C
stochastically stable for all three √
dimensions. √Like in Regime I, we also test LEnKF with
1
ǫ
small scale system noises σx = ǫσx , σoǫ = ǫσo , where ǫ = 1, 12 , 41 , . . . , 32
. Subplot d)
indicates the LEnKF has the correct MSE scaling with ǫ.
0.55

Covariance

0.5

0.45

0.4

0.35

0

18

5

10

15

20

25

30

35

40

45

50

4

Concentration of localized random matrices

In this section, we present the proof of Theorem 2.1. While part a) is more useful, it can be
established easily from part b), using a similar argument as in [21].

4.1

Entry-wise concentration

It is well known that the averages of independent Gaussian variables concentrate around
their expected values. In specific, a simplified version of theorem 1.1 from [44] is:
Theorem 4.1 (Hanson-Wright inequality). Let ξ ∼ N (0, In ) and A be an n × n matrix.
Then for any t ≥ 0




t
t2
T
T
.
,
P |ξ Aξ − Eξ Aξ| > t ≤ 2 exp −c min
kAk2HS kAk

Here c is a constant independent
of other parameters. The Hilbert-Schmidt (Frobenius) norm
P
is denoted by kAkHS = [ i,j [A]2i,j ]1/2 .
This provides us a straight forward way to control the random matrix entries [Z]i,j in
Theorem 2.1.

Lemma 4.2. Under the conditions of Theorem 2.1, let ∆ = Z − EZ. There is an absolute
constant c such that for any t ≥ 0,
P(|[∆]i,j | > σa,z t) ≤ 8 exp(−cK min{t, t2 }).
Proof. For any vector u, denote ∆u = uT [Z − EZ]u. Then by symmetry,
1
[∆]i,j = (∆ei +ej − ∆ei −ej ).
4
Recall that ei is the i-th standard basis vector. So it suffices to find a concentration bound
for ∆u with u = ei ± ej . To do that, note that uT Σz u = EuT zk zkT u, so we can decompose
∆u
∆u = K

−1

= 2K

K
X

[hu, ak + zk ihu, ak + zk i − hu, ak ihu, ak i − Ehu, zk ihu, zk i]

k=1
K
X
−1
k=1

hu, ak ihu, zk i + K

−1

K
X
k=1

(hu, zk ihu, zk i − Ehu, zk ihu, zk i)

We denote ha, bi = aT b as the inner product, and the P
two summations above as I and II in
K
T
−1
2
T
the following. Notice that hu, zk i ∼ N (0, u Σz u), K
k=1 hej , ak i = u Σa u. Moreover for
u = ei ± ej ,
uT Σz u = [Σz ]i,i + [Σz ]j,j ± 2[Σz ]i,j ≤ 2([Σz ]i,i + [Σz ]j,j ) ≤ 4σa,z .
19

(4.1)

We have the same conclusion for uT Σa u. Because hu, ak i is a deterministic scalar,
hu, ak ihu, zk i ∼ N (0, uT ak aTk u · uT Σz u)
and
I = 2K

−1

K
X
k=1

hu, ak ihu, zk i ∼ N (0, 4K −1uT Σa u · uT Σz u).

2
Because by definition of σa,z , uT Σa u · uT Σz u ≤ 16σa,z
, by the Chernoff bound for Gaussian
distributions, there is a c1 > 0 so that

P(|I| > 12 σa,z t) ≤ 2 exp(−c1 Kt).
In order to deal with II, notice that
ξ := √
So
II = K

−1

1
uT Σz u
K
X
k=1

[hu, z1 i, · · · , hu, zK i]T ∼ N (0, IK ).

(hu, zk i2 − Ehu, zk i2 ) = ξ T Aξ − Eξ T Aξ,

where A = K1 (uT Σz u)IK . Clearly, kAk ≤ 4σKa,z , and kAk2HS ≤
4.1 there is a constant c2 so that for all s ≥ 0

2
16σa,z
.
K

Therefore, by Theorem

2

s
}).
P(|II| > 21 s) ≤ 2 exp(−c2 min K{ σs2 , σa,z
a,z

−1
Let t = σa,z
s, the inequality can be written as

P(|II| > 21 σa,z t) ≤ 2 exp(−c2 K min{t, t2 }).
Because |∆u | ≤ |I| + |II|, by the union bound, if we let c = min{c1 , c2 } ,
P(|∆u | > σa,z t) ≤ P(|I| > 21 σa,z t) + P(|II| > 21 σa,z t) ≤ 4 exp(−cK min{t, t2 }).
Finally, recall the bound above holds for all u = ei ± ej , so by (4.1)
P(|[∆]i,j | > σa,z t) ≤ P(|∆ei +ej | > σa,z t) + P(|∆ei −ej | > σa,z t) ≤ 8 exp(−cK min{t, t2 }).
Entry-wise concentration now comes as a direct corollary.
Proof of Theorem 2.1 b). Let ∆ = Z − EZ. Note that kZ − EZk∞ = maxi,j=1,...,d {|[∆]i,j |},
so using the previous lemma we have our claim by the union bound
X
P(kZ − EZk∞ > σa,z t) ≤
P(|[∆]i,j | > σa,z t) ≤ 8d2 exp(−cK min{t, t2 }).
i,j

20

4.2

Summation of entry-wise deviation

One simple fact of matrix norm is that k∆k ≤ k∆k1 . This is also exploited by [21]
Lemma 4.3. Given a matrix ∆, the following holds
a) If ∆ is symmetric, then
k∆k ≤ k∆k1 = max
i

( d
X
j=1

)

|[∆]i,j | .

b) k∆k∞ ≤ k∆k always holds. If in addition ∆ has bandwidth l, then k∆k ≤ Bl k∆k∞ .
Proof. For a) part, recall that ei is the i-th standard basis vector. Notice that
±(ei eTj + ej eTi )  ei eTi + ej eTj .
Therefore
∆=

X

1X
[∆]i,j (ei eTj + ej eTi )
2 i6=j

[∆]i,i ei eTi +

i



d
X

[∆]i,i ei eTi

i=1

d

d

XX
1X
|[∆]i,j |(ei eTi + ej eTj ) =
|[∆]i,j |ei eTi  k∆k1 Id .
+
2 i6=j
i=1 j=1

For the b) part, by the definition of operator norm, and kei k = kej k = 1, we have
|ei ∆eTj | ≤ k∆k.
Taking maximum among all i and j, we have P
k∆k∞ ≤ k∆k.
Next note that k∆k = k∆∆T k1/2 ≤ maxi j |[∆∆T ]i,j |, and if ∆ is of bandwidth l, by
part a)
X
X
X
|[∆∆T ]i,j | ≤
|[∆]i,k ||[∆]j,k | ≤ Bl2 k∆k2∞ .
j

k:d(i,k)≤l j:d(j,k)≤l

Therefore k∆k ≤ Bl k∆k∞ .
Now the Theorem 2.1 a) comes as a direct corollary:
Proof of Theorem 2.1 a). Let ∆ = Z − EZ. By Lemma 4.3 a),
( d
)
X
k∆◦DL k ≤ k∆◦DL k1 = max
[DL ]i,j |[∆]i,j | ≤ kDL k1 max |[∆]i,j | = kDL k1 kZ −EZk∞ .
i

i,j

j=1

Therefore by part b) of this theorem,
P(k∆ ◦ DL k ≥ kDL k1 σa,z t) ≤ P (kZ − EZk∞ ≥ σa,z t) ≤ 8 exp(2 log d − cK min{t, t2 }).

21

5
5.1

Error analysis of LEnKF
Localization inconsistency

Lemma 5.1. Fix an L > l, if matrix A is of bandwidth l, the difference caused by commuting
localization and bilinear product with A
∆ = A[C ◦ DLcut ]AT − [ACAT ] ◦ DLcut
has nonzero entries only for indices (i, j) with |d(i, j) − L| ≤ 2l.
If in addition, matrix C follows an (M, Φ, L)-localized structure, then
|[∆]i,j | ≤ MΦ(L − 2l)kAk2∞ Bl2 ,

L − 2l ≤ d(i, j) ≤ L.

Recall that Bl is the volume constant given by (2.8).
Proof. By the matrix product rule,
X
X
[∆]i,j =
[A]i,u [C]u,v [A]j,v − 1d(i,j)≤L
[A]i,u [C]u,v [A]j,v .

(5.1)

u,v

d(u,v)≤L

If d(i, j) > L + 2l, note that [A]i,u [A]j,v 6= 0 only when d(i, u) ≤ l, d(j, v) ≤ l. But for
these terms, by the triangular inequality d(u, v) > L, and they are not included in (5.1).
Therefore (5.1)= 0.
P
If d(i, j) ≤ L, it is easy to verify that [∆]i,j = − d(u,v)>L [A]i,u [C]u,v [A]j,v . Moreover,
[A]i,u [A]j,v 6= 0 only when d(i, u) ≤ l, d(j, v) ≤ l. So if d(i, j) < L − 2l, then by triangular
inequality d(u, v) < L and [∆]i,j = 0.
Next, we assume C follows an
P(M, Φ, L)-localized structure. If L < d(i, j), then among
the nonzero terms in [∆]i,j =
d(u,v)≤L [A]i,u [C]u,v [A]j,v , d(u, v) ≥ L − 2l by triangular
inequality. This leads to
X
|[∆]i,j | ≤
kAk2∞ MΦ(L − 2l) ≤ Bl2 kAk2∞ MΦ(L − 2l).
u,v:d(u,v)≤L,d(i,u)≤l,d(j,u)≤l

Here we used that
#{(u, v) : d(u, v) ≤ L, d(v, i), d(u, i) ≤ l} ≤ #{(u, v) : d(v, i), d(u, i) ≤ l} = Bl2 .
P
If L − 2l ≤ d(i, j) ≤ L, then by [∆]i,j = − d(u,v)>L [A]i,u [C]u,v [A]j,v ,
|[∆]i,j | ≤ MΦ(L)

X

d(u,v)>L

|[A]i,u ||[A]j,v | ≤ MΦ(L)kAk2∞ Bl2 ,

where we applied the inequality
#{(u, v) : d(v, i), d(u, i) ≤ l, d(u, v) > L} ≤ #{(u, v) : d(v, i), d(u, i) ≤ l} = Bl2 .
In either case, we have the bound we claim, since Φ(L) ≤ Φ(L − 2l).
22

Proof of Proposition 2.3. Since Schur product is a linear operation, we can decompose the
localization inconsistency as
b n H)][C
bn ◦ DL ][(I − K
b n H)T AT ]
∆loc =[An (I − K
cut
n


T
T
b n H)]C
bn [(I − K
b n H) A ] ◦ DL
− [An (I − K
n
cut
2
T
T
2
T
T
bnK
b A + Σn ] ◦ DL
bnK
b A + Σn ] − [σ An K
+ [σ An K
o

n

n

o

n

n

cut

b n and Σn are of bandwidth at most l, An K
bnK
b T AT has bandwidth at most 4l
Since both K
n n
by triangular inequality. Since L ≥ 4l, so
bnK
b T AT + Σn ] = [σ 2 An K
bnK
b T AT + Σn ] ◦ DL ,
[σo2 An K
n n
o
n n
cut

In other words, ∆loc is



b n H)][C
bn ◦ DL ][(I − K
b n H)T AT ] − [An (I − K
b n H)]C
bn [(I − K
b n H)T AT ] ◦ DL ,
[An (I − K
cut
n
n
cut

b n H)k∞ . Recall that
which can be applied by Lemma 5.1. Next, we try to bound kAn (I − K
kHk = 1, kAn k ≤ MA and Lemma 4.3 b),
b n H)k∞ ≤ kAn (I − K
b n H)k ≤ MA kI − K
b n Hk ≤ MA (1 + kK
b n Hk)
kAn (I − K

b n H has bandwidth l. To see this, note that
In domain localization (2.10), K

b n H]i,j = [K
b ni H]i,j = [C
bni H T (σo2 Iq + H C
bni H T )−1 H]i,j
[K
X
bi ]i,o [(σ 2 Iq + H C
b i H T )−1 ]k,m 1j=om .
=
[C
n
o
n
k

(5.2)

m,k

bi has nonzero entries only in Ii × Ii ,
Since C
n

bni H T )−1 ]k,m = σo−2 1k=m
[(σo2 Iq + H C

if d(ok , i) > l or d(om , i) > l.

b n H]i,j = 0 if d(i, j) > l.
bni ]i,o = 0 if d(ok , i) > l. Therefore, [K
Also [C
k
b n Hk ≤ Bl kK
b n Hk∞ . Since the i-th row of K
b n H is the i-th row of
By Lemma 4.3 b), kK
i
Kn H, so by Lemma 4.3 b),
b n Hk∞ ≤ max{kK i Hk∞ } ≤ max{kK i Hk}.
kK
n
n
i

i

Moreover, by definition (2.9) and Lemma 4.3 a)

bi kk(σ 2 I + H C
b i H T )−1 k ≤ σ −2 kC
bi k ≤ σ −2 kC
b i k1 .
kKni k ≤ kC
n
o
n
o
n
o
n

bi has nonzero entries only in Ii × Ii , by Lemma 4.3,
Note that C
n
bi k1 ≤ Bl kC
bi k∞ ≤ Bl kC
bn k∞ .
kC
n
n
23

bn follows an (Mn , Φ, L) structure, kC
bn k∞ ≤ Mn . Summing up, the domain
Moreover, since C
localized Kalman gain can be bounded by
b n H)k∞ ≤ MA (1 + σ −2 Bl Mn ).
kAn (I − K
o

Then by Lemma 5.1, the localization inconsistency matrix is bounded entry-wise
|[∆]i,j | ≤ Mn MA2 (1 + σo−2 Bl Mn )2 Bl2 Φ(L − 2l),
while |[∆]i,j | = 0 if |d(i, j)−L| > 2l. So there are at most BL,l = maxi #{j, |d(i, j)−L| ≤ 2l}
nonzero entries in each row.
As a consequence
k∆loc k ≤ k∆loc k1 ≤ Mn MA2 (1 + σo−2 Bl Mn )2 Bl2 BL,l Φ(L − 2l).

5.2

Component information gain through filtering

One of the fundamental properties in Kalman filter is that the assimilation of observation
improves estimation. Mathematically, this can be represented by that the forecast covariance
matrix dominates the posterior covariance matrix. Unfortunately, with LEnKF, this natural
bn  (I − K
b n H)C
bn (I − K
b n H)T + σ 2 K
b bT
property, C
o n Kn , may no longer hold. However, we can
still show the dominance at the diagonal entries.
Proposition 5.2. The assimilation step lowers the variance at each component:
bn ]i,i ≥ [(I − K
b n H)C
bn (I − K
b n H)T + σ 2 K
b bT
[C
o n Kn ]i,i ,

i = 1, · · · , d.

bn is updated through the Kalman gain matrix
Proof. Recall that the i-th coordinate of ∆X
i
b
Kn . Therefore,
(k)

2 bi bi T
b n H)C
bn (I − K
b n H)T + σ 2 K
b bT
bi
b
bi T
[(I − K
o n Kn ]i,i = [(I − Kn H)Cn (I − Kn H) + σo Kn (Kn ) ]i,i

b ni H]i,j 6= 0 only when d(i, j) ≤ l, so
Moreover, in (5.2) we have shown that [K

b ni H)C
bn (I − K
b ni H)T + σo2 K
b ni (K
b ni )T ]i,i = [(I − K
b ni H)C
bni (I − K
b ni H)T + σo2 K
b ni (K
b ni )T ]i,i .
[(I − K

Note that the right side is the posterior Kalman covariance with the forecast covariance
bi . Therefore by
being C
n
b ni H)C
bni (I − K
b ni H)T + σo2 K
b ni (K
b ni )T = C
bni − C
bni H T (σo2 Iq + H C
bni H T )−1 H C
bni  C
bni ,
(I − K

we have

b n H)C
bn (I − K
b n H)T + σ 2 K
b bT
bi
b
[(I − K
o n Kn ]i,i ≤ [Cn ]i,i = [Cn ]i,i .
24

5.3

Sampling error

First, we have the following general integral lemma
Lemma 5.3. If Y is a nonnegative random variable that satisfies
P(Y > Mt) ≤ 8d2 exp(−cK min{t, t2 )),

c > 0, M ≥ 1.

Then for any δ ∈ (0, 1), if K ≥ Γ(Mδ −1 , d), where

x, 18
x2 log d}.
Γ(x, d) = max{9x2 , 24
c
c

We have EY ≤ δ and EY 2 ≤ 2Mδ.
Proof. Let ǫ =

δ
,
3M

and X = Y /M, we have K ≥ max{ǫ−2 , cǫ8 , cǫ22 log d}, and
P(X > t) ≤ 8d2 exp(−cK min{t, t2 )),

We will show that EX ≤ 3ǫ and EX 2 ≤ 6ǫ, which are equivalent to our
R ∞ claims. Recall the
integration by part formula for nonnegative random variables, EX = 0 P(X > x)dx,
Z ǫ
Z ∞
EX =
P(X > x)dx +
P(X > x)dx
0
ǫ
Z ∞
≤ǫ+
P(X > t)dt
ǫ
Z ∞
Z 1
2
≤ǫ+8
d exp(−cKt)dt + 8
d2 exp(−cKt2 )dt
Z1 ∞
Zǫ ∞
2
≤ǫ+8
d exp(−cKt)dt + 8
d2 exp(−cKt2 )dt.
ǫ

ǫ

Note that with our requirement on K, d2 exp(−cKǫ) ≤ 1,
Z ∞
8
8d2
exp(−cKǫ) ≤
≤ ǫ.
8d2 exp(−cKt)dt =
cK
cK
ǫ

And for t > ǫ, 8 ≤ 2ǫcKt, so
Z ∞
Z
2
8 exp(−cKt )dt ≤ ǫ
ǫ

ǫ

∞

2cKt exp(−cKt2 )dt = ǫ exp(−cKǫ2 ) ≤ ǫ.

2

As for EX , we again apply the integration by part formula
Z ∞
2
EX =
2tP(X ≥ t)dt
0
Z ∞
≤ 2ǫ +
2tP(X ≥ t)dt
ǫ
Z ∞
Z ∞
2
2
≤ 2ǫ + 8d
2t exp(−cKt)dt + 8d
2t exp(−cKt2 )dt
ǫ

= 2ǫ +
≤ 2ǫ +

2

ǫ
+ c21K 2 )
16d exp(−cKǫ)( cK
16ǫ
8
+ c216
+ cK
≤ 6ǫ.
cK
K2

We used K ≥ max{ǫ−2 , cǫ8 , cǫ22 log d} in the last line.
25

+

ǫ
8d2
cK

exp(−cKǫ2 )

bn follows (Mn , Φ, L)-localized
Corollary 5.4. Under condition 1) of Theorem 2.4, suppose C
structure. For any ǫ ∈ (0, 1), if

a) K > Γ(BL ǫ−1 , d), then the sampling error

bn+1 − rRn (C
bn )) ◦ DL k ≤ ǫ(B2 M 2 Mn + MΣ ),
En k(C
cut
l
A

b) K > Γ(rCǫ−1 , d) for any C ≥ 1, then the entry-wise sampling error
bn+1 − rRn (C
bn )k∞ ≤ ǫC −1 kRn (C
bn )k∞ .
En kC

bn+1 − rRn (C
bn )k2 ≤ ǫ2C −1 kRn (C
bn )k2 .
En kC
∞
∞

Proof. We apply Theorem 2.1 with
√
b n H)∆X
b (k) ,
ak = rAn (I − K
n

√

zk =

and DL = DLcut . Then

b n H)C
bn (I − K
b n H)T AT ,
Σa = rAn (I − K
n

b n ζ (k) +
rAn K
n

√

rξn(k) ,

b nK
b T AT + rΣn .
Σz = rσo2 An K
n n

bn ), where recall
Note that Σa  Σa + Σz and Σz  Σa + Σz = rRn (C
bn ) = An Qn AT + Σn ,
Rn (C
n

Therefore

b n H)C
bn (I − K
b n H)T + σ 2 K
b bT
Qn := (I − K
o n Kn .

1/2
1/2
bn )k∞ .
σa,z ≤ r max{[Σa ]i,i , [Σa ]i,i [Σz ]j,j } ≤ r max[Σa + Σz ]i,i = rkRn (C
i

i,j

Moreover, since Qn is positive semidefinite (PSD), so
q
kQn k∞ = max |[Qn ]i,j | ≤ max[Qn ]i,i max[Qn ]j,j = max[Qn ]i,i ≤ kQn k∞ .
i,j

i

Moreover, by Proposition 5.2,

j

i

bn ]i,i ≤ Mn .
[Qn ]i,i ≤ [C

bn ) is PSD, and by Lemma 4.3 kAn k∞ ≤ kAn k ≤ MA ,
Since Rn (C
(
)
X
bn )k∞ ≤ max[Rn (C
bn )]i,i = max [Σn ]i,i +
kRn (C
[An ]i,j [Qn ]j,k [An ]i,k
i

i

≤

Apply Theorem 2.1, since kDLcut k1 = maxi

MA2 Bl2 Mn

P

j

+ MΣ .

j:d(i,j)<L

1 = BL , we have that

bn+1 − rRn (C
bn )) ◦ DL k/kRn (C
bn )k∞ > rBL t) ≤ 8d2 exp(−cK min{t, t2 }).
Pn (k(C
cut
bn+1 − rRn (C
bn )k∞ /kRn (C
bn )k∞ > rt) ≤ 8d2 exp(−cK min{t, t2 }).
Pn (kC

Pn denotes the probability conditioned on Fn . Apply Lemma 5.3 with the both of them,
but using δ = ǫ for the first inequality and δ = ǫC −1 for the second, we have our claimed
results.
26

5.4

Error analysis

Next, we proceed to prove Theorem 2.4.
Proof of Theorem 2.4. For each time n, let rn be the smallest number such that the following
hold,
bn ◦ DL + ρId ), rn ≥ 1.
ES ên ⊗ ên  rn (C
cut
We will try to find a recursive upper bound of rn+1 in term of rn .
Step 1: tracking the filter error. Recall that the forecast error at time n + 1 is provided
by the (2.15), and its covariance conditioned on sample noise realization is
b n H)ES ên ⊗ ên (I − K
b n H)T + σ 2 K
b bT T
ES ên+1 ⊗ ên+1 = An [(I − K
o n Kn ]An + Σn
b n H)(C
bn ◦ DL )(I − K
b n H)T AT
 rn An (I − K
cut

+

b nK
b T AT
[σo2 An K
n n

n

b n H)(I − K
b n H)T AT + Σn ].
+ rn ρAn (I − K
n

By Young’s inequality (a + b)(a + b)T  2aaT + 2bbT , and that HH T = Iq ,

b n H)(I − K
b n H)T AT ≤ 2(An AT + An K
b n HH T K
b T AT )
An (I − K
n
n
n n
T T
T
b
b
≤ 2(An A + An Kn K A ).
n

Moreover, An ATn  MA2 Id 

Furthermore,

2
MA
Σ .
mΣ n

Denote DΣ = max{

n

2
2MA
, 2 },
mΣ σo2

n

then

b n H)(I − K
b n H)T AT  DΣ (Σn + σ 2 An K
bnK
b T AT ).
An (I − K
n
o
n n

b n H)(C
bn ◦ DL )(I − K
b n H)T AT + (1 + rn ρDΣ )(σ 2 An K
bnK
b T AT + Σn ).
ES ên+1 ⊗ ên+1  rn An (I − K
cut
n
o
n n
bn ) in (2.16) is
Recall that R′n (C

bn ) = An (I − K
b n H)(C
bn ◦ DL )(I − K
b n H)T AT + σ 2 An K
bnK
b T AT + Σn .
R′n (C
cut
n
o
n n

Therefore

bn ).
ES ên+1 ⊗ ên+1  max{1, rn /r, (1 + rn ρDΣ )/r} · rR′n (C

With our condition 2) on ρ,

(1 + rn ρDΣ )/r ≤

1
r

+

r−1 rn
r r

≤ max{1, rn /r},

bn ).
so ES ên+1 ⊗ ên+1  max{1, rn /r}rR′n (C

Step 2: difference between filter error covariance and its estimate.
bn+1 . Its conditional
The EnKF estimates the error covariance by the ensemble covariance C
expectation is
bn+1 = rRn (C
bn ) = r(An (I − K
b n H)C
bn (I − K
b n H)T AT + σ 2 An K
bnK
b T AT + Σn ).
En C
n
o
n n
27

(5.3)

In order to establish a control of the new filter error using localized ensemble covariance
matrix, consider the difference
bn+1 ◦ DL − rR′ (C
bn ) = (C
bn+1 ◦ DL − En C
bn+1 ◦ DL ) + r(Rn (C
bn ) ◦ DL − R′ (C
bn )).
C
cut
n
cut
cut
cut
n

The first part of (5.3) is the error caused by sampling. By Corollary 5.4, if we denote
bn+1 ◦ DL − En C
bn+1 ◦ DL k
µn+1 := kC
cut
cut

then En µn+1 ≤ (Bl2 MA2 Mn + MΣ )δ/r if K satisfies condition 5).
The second part of (5.3) is the localization inconsistency. By Proposition 2.3, we have
bn ) ◦ DL − R′ (C
bn )k ≤ Mn M 2 (1 + σ −2 Bl Mn )2 B2 BL,l Φ(L − 2l) =: νn+1 .
kRn (C
cut
n
A
o
l

Summing these two parts up,

bn )  C
bn+1 ◦ DL + r(µn+1 + νn+1 )Id .
rR′n (C
cut

Then

bn )  (1 + r (µn+1 + νn+1 ))(C
bn+1 ◦ DLcut + ρId ).
rR′n (C
ρ

bn ), so if we let rn+1 be the
Recall that in step 1, we have ES ên+1 êTn+1  max{1, rrn }rR′n (C
smallest number such that
bn+1 ◦ DLcut + ρId ),
ES ên+1 ⊗ ên+1  rn+1 (C

then

rn+1 ≥ 1,

rn+1 ≤ max{1, rrn }(1 + ρr (µn+1 + νn+1 )).

Step 3: long time stability analysis. Since r∗ ≤ r
max{0, log(rn /r)} ≤ max{0, log(rn /r∗ )} ≤ log rn − log r∗ 1rn ≥r∗ .
Taking the logarithm of (5.4), and using that log(1 + x + y 3 ) ≤ x + 2y for all x, y ≥ 0,
log rn+1 ≤ log rn − log r∗ 1rn ≥r∗ + log(1 + ρr (µn+1 + νn+1 ))
≤ log rn − log r∗ 1rn ≥r∗ + ρr µn+1 + 2( ρr νn+1 )1/3 .

Sum this inequality from n = 0, . . . , T − 1, we have
log r∗

T −1
X
n=0

1rn ≥r∗ ≤ log r0 − log rT +

T −1
X

( ρr µn+1 + 2( ρr νn+1 )1/3 ).

n=0

Because rT ≥ 1,
T −1
X
n=0

1rn ≥r∗

T −1
1 X r
log r0
( ρ µn+1 + 2( ρr νn+1 )1/3 ).
+
≤
log r∗ log r∗ n=0

28

(5.4)

Take expectation,
T −1
X
n=0

P(rn ≥ r∗ ) = E

T −1
X

1rn ≥r∗

n=0

T −1
log r0
1 X r
≤
+
( Eµn+1 + 2E( ρr νn+1 )1/3 ).
log r∗ log r∗ n=0 ρ

(5.5)

Step 4: Upper bounds for (5.5). Recall in step 2 we have that
T −1
X

r
Eµn+1
ρ

n=0

T −1
X

≤

δ
(Bl2 MA2 EMn
ρ

+ MΣ ).

n=0

Next, note the following holds because Bl ≥ 1
νn+1 = MA2 Bl2 BL,l Φ(L − 2l)Mn (1 + σo−2 Bl Mn )2 ≤ MA2 σo2 Bl3 BL,l Φ(L − 2l)(1 + σo−2 Bl Mn )3 .
With condition 4), we have
2/3

1/3

MA BL,l Bl2 Φ1/3 (L − 2l) ≤ δ,
so

1/3

2/3

1/3

Eνn+1 ≤ EMA σo2/3 BL,l Bl Φ1/3 (L − 2l)(1 + σo−2 Bl Mn ) ≤ δ(σo2/3 + σo−1/3 EMn ).

In conclusion,

1/3

2E( ρr νn+1 )1/3 ≤ 2δ ρr 1/3 (σo2/3 + σo−1/3 EMn ).
Plug these bounds to (5.5), and then use (2.20)
T −1
r0
δ
1X
P(rn ≥ r∗ ) ≤
+
(ρ−1 Bl2 MA2 +
T n=0
T log r∗ T log r∗

2r 1/3
)
(ρσo )1/3

T −1
X

EMn +

n=0

δ
1/3
(ρ−1 MΣ + 2 ρr 1/3 σo2/3 )
log r∗

r0
δ(B0 kC0 k + D0 ) −1 2 2
2r 1/3
≤
+
(ρ Bl MA + (ρσ
1/3 )
o)
T log r∗
T log r∗

δ  −1 2 2
−1
2r 1/3
r 1/3 2/3
(ρ Bl MA + (ρσ
.
)M
+
ρ
M
+
2
σ
+
0
Σ
1/3
1/3
o
ρ
o)
log r∗

For our result, simply notice that
rn ≤ r∗

6

⇔

bn + ρId ).
ES ên ⊗ ên  r∗ (C

Localized covariance for linear LEnKF systems

As discussed in the introduction, the existence of a localized covariance structure is often
assumed in practice to motivate the localization technique. Our result, Theorem 2.4, shows
that such a structure indeed can guarantee estimated performance, assuming the parameters
and sample size are properly tuned. Then it is natural to ask when does a stable localized
structure exist. This is an interesting and important question by itself, but to answer it
for general signal-observation systems with rigorous proof is beyond the scope of this paper.
Here we demonstrate how to verify a stable localized covariance for simple linear models.
29

6.1

Localized covariance propagation with weak local interactions

As discussed in Theorem 2.4, we require An to be of a short bandwidth l. In other words,
interaction in one time step exists only for components of distance l apart. When l = 1,
this type of interaction is often called nearest neighbor interaction, and it includes many
statistical physics models with proper spatial discretization.
Generally speaking, localized covariance is formed through weak local interactions. With
linear dynamics described by An , one way to enforce a weak local interaction is through
(2.22). We will show in this subsection that weak local interaction propagates a localized
bn ]i,j ∝ λd(i,j) , from diagonal entries of the covariance matrix
covariance structure of form [C
A
to entries further away from diagonal.
bn and Cn , we define the
To describe the state of localization in covariance matrices C
following quantities
n
o
n
o
cn,l = max |[C
bn ]i,j |λ−d(i,j)∧l , Mn,l = max |[Cn ]i,j |λ−d(i,j)∧l .
M
(6.1)
A
A
i,j

i,j

Then clearly, the forecast covariance matrices follow the (Mn , λxA , L) localized structure with
cn,L . The goal of this section is to show that M
cn,L is a stable stochastic sequence.
Mn = M
The following properties hold immediately because the matrices involved are PSD.

bn , define Mn,l , M
cn,l as in
Lemma 6.1. Given positive semidefinite (PSD) matrices Cn , C
cn,0 = maxi [C
bn ]i,i ,
(6.1), we have M
cn,0 ≤ M
cn,1 ≤ · · · ≤ M
cn,k ≤ M
cn,0 λ−k .
M
A

The same properties also hold for Mn,k as well.

bn ]i,j is the ensemble covariance, so for i 6= j
Proof. Recall that [C
q
bn ]i,i ||[C
bn ]j,j | ≤ max[C
bn ]i,i .
b
|[Cn ]i,j | ≤ |[C
i

Therefore

cn,0 = max |[C
bn ]i,j | = max[Cbn ]i,i .
M
i,j

i

cn,k in k is quite obvious since d(i, j) ∧ k ≤ d(i, j) ∧ (k + 1), and
The monotonicity of M
o
n
−d(i,j)∧k
b
c
b
≤ λ−k
Mn,k = max |[Cn ]i,j |λA
A max |[Cn ]i,j |.
i,j

i,j

Next, we investigate how does the forecast step change the state of localization.
Proposition 6.2. Suppose Σn = σξ2 Id and the linear dynamics admits a weak local interaction satisfying (2.22), the forecast step propagates the localization in covariance. In particbn+1 = An Cn AT + Σn , then the localization
ular, given any covariance matrix Cn , and let C
n
states described by (6.1) follows
cn+1,0 ≤ λ2 Mn,0 + σ 2 ,
M
A
ξ
30

cn+1,k ≤ max{λ2 Mn,k , M
cn+1,0 },
M
A

cn+1,k+1 ≤ max{λA Mn,k , M
cn+1,0 }.
M

Proof. Note that [Cbn+1 ]i,j = [An Cn ATn ]i,j + σξ2 1i=j . Moreover
X
|[An ]i,m [An ]j,m′ [Cn ]m,m′ |
|[An Cn ATn ]i,j | ≤
≤
≤

X

m,m′

X

m,m′

m,m′

d(i,j)∧k

Mn,k λA

−d(j,m′ )

Mn,k λA

|[An ]j,m′ |λA

−d(i,m)

|[An ]j,m′ |λA

|[An ]i,m |λA

= Mn,k λA

−d(j,m′ )

−d(i,m)

|[An ]i,m |λA

X
m

−d(i,m)

|[An ]i,m |λA

!

d(i,m)+d(j,m′ )+d(m,m′ )∧k

d(i,j)∧k

X
m

−d(j,m)

|[An ]j,m |λA

!

,

d(i,j)∧k

which by (2.22) is bounded by λ2A Mn,k λA
.
By Lemma 6.1,
bn+1 ]i,i ≤ λ2A Mn,0 + σξ2 .
cn+1,0 = max[C
M
i

Moreover,



n
o
−d(i,j)∧k
c
b
b
bn+1 ]i,i .
Mn+1,k = max max[Cn+1 ]i,j λA
, max[Cn+1 ]i,i ≤ max λ2A Mn,k , max[C
i6=j

i

i



n
o
−d(i,j)∧(k+1)
c
b
b
bn+1 ]i,i .
Mn+1,k+1 = max max[Cn+1 ]i,j λA
, max[Cn+1 ]i,i ≤ max λA Mn,k , max[C
i6=j

6.2

i

i

Preserving a localized structure with sparse observations

From now on, we require the observations to be sparse in the sense that d(oi , oj ) > 2l for any
i 6= j. Then for each location i ∈ {1, · · · , d}, there is at most one location o(i) ∈ {o1 , · · · , oq }
such that d(i, o(i)) ≤ l. If such an o(i) doesn’t exist, we set o(i) = nil, the analysis step will
not update it, and we will see the discussion for these components are trivial.
With domain localization and sparse observations, the analysis step updates the information at the i-th component using only the observation at o(i). This significantly simplifies
bi H T + σ 2 Iq )−1 , which is diagonal with entries (σ 2 + [C
bn ]o ,o )−1 in
the formulation of (H C
n
o
o
i i
Ii × Ii . As a result, the Kalman update matrix has entries

 [Cbn ]i,o(i) ,
j = o(i);
i
b n H]i,j = [K
b H]i,j = σo2 +[Cbn ]o(i),o(i)
[K
n
0,
else.
In fact, if we apply the covariance localization scheme instead of domain localization, the
Kalman gain remains the same in this setting.
In below, we investigate how does the assimilation step change the state of localization.
31

bn , define K
b n as the Kalman gain in (2.10),
Proposition 6.3. Given any covariance matrix C
and let
b n H)C
bn (I − K
b n H)T + σ 2 K
b bT
Cn = (I − K
o n Kn .

Define the state of localization using (6.1). Then
cn,0 ,
Mn,0 ≤ M

where

cn,k )
Mn,k ≤ φ(M

φ(M) = M(1 + σo−2 M)2 + σo−2 M 2 .
cn,0 = maxi |[C
bn ]i,i |, so Mn,0 ≤ M
cn,0
Proof. Based on Lemma 6.1, Mn,0 = maxi |[Cn ]i,i |, M
holds by Proposition 5.2. Next, we look at the off diagonal terms:
b
b
b
b
bn ]i,j − [Cn ]i,o(i) [Cn ]j,o(i) − [Cn ]i,o(j) [Cn ]j,o(j)
[Cn ]i,j = [C
bn ]o(i),o(i)
bn ]o(j),o(j)
σo2 + [C
σo2 + [C
bn ]i,o(i) [C
bn ]j,o(j)[C
bn ]o(i),o(j)
[C
+
bn ]o(i),o(i) )(σo2 + [C
bn ]o(j),o(j) )
(σo2 + [C
bn ]i,o(i) [C
bn ]j,o(i)
σ 2 [C
+ o
1o(i)=o(j) .
(σo2 + [Cbn ]o(i),o(i) )2

(6.2)

We have the following bounds for each term in (6.2)

bn ]i,o(i) [C
bn ]j,o(i)
[C
c2 λd(i,o(i))∧k+d(j,o(i))∧k ≤ σ −2 M
c2 λd(j,i)∧k .
≤ σo−2 M
n,k A
o
n,k A
2
b
σ + [Cn ]o(i),o(i)
o

bn ]i,o(i) [C
bn ]j,o(j)[Cbn ]o(i),o(j)
[C
bn ]o(i),o(i) )(σ 2 + [C
bn ]o(j),o(j) )
(σ 2 + [C
o

≤

o
d(i,o(i))∧k+d(j,o(j))∧k+d(o(j),o(i))∧k
−4 c3
σo Mn,k λA

c3 λ
≤ σo−4 M
n,k A

d(i,j)∧k

.

bn ]i,o(i) [C
bn ]j,o(i)
[C
c2 λd(i,o(i))∧k+d(j,o(i))∧k ≤ σ −4 M
c2 λd(i,j)∧k .
≤ σo−4 M
n,k A
o
n,k A
2
2
bn ]o(i),o(i) )
(σ + [C
o

In summary

cn,k [(1 + σ −2 M
cn,k )2 + σ −2 M
c2 ]λd(i,j)∧k = φ(M
cn,k )λd(i,j)∧k .
|[Cn ]i,j | ≤ M
o
o
n,k A
A
b
b
b
Proposition 6.4. Denote δn+1 = λ−L
A kCn+1 − rRn (Cn )k∞ /kRn (Cn )k∞ , and
o
n
2
ψλA (M, δ) = (r + δ) max λA M 1 + σo−2 M + λA σo−2 M 2 , λ2A M + σξ2 .
Then for k ≤ L − 1,

cn+1,0 ≤ (r + δn+1 )(λ2 Mn,0 + σ 2 ),
M
A
ξ
32

cn+1,k+1 ≤ ψλ (M
cn,k , δn+1 ).
M
A

Proof. Recall that
bn ) = An [(I − K
b n H)C
bn (I − K
b n H)T + σ 2 K
b bT T
Rn (C
o n Kn ]An + Σn .

Following (6.1), we define its localized status:
n
o
−d(i,j)∧l
T
2 b bT
b
b
b
Rn,l = max |[(I − Kn H)Cn (I − Kn H) + σo Kn Kn ]i,j |λA
,
i,j

n
o
bn+1,l = max |[Rn (C
bn )]i,j |λ−d(i,j)∧l .
R
A
i,j

Apply Proposition 6.3,

cn,0 ,
Rn,0 ≤ M

Then apply Proposition 6.2, we find that

2
c
bn+1,0 = kRn (C
bn )k∞ ≤ λ2 M
R
A n,0 + σξ ,

cn,k ).
Rn,k ≤ φ(M

bn+1,k+1 ≤ max{λA φ(M
cn,k ), R
bn+1,0 }.
R

Finally by Lemma 6.1,

bn )k∞ .
cn+1,0 = kC
bn+1k∞ ≤ rkRn (C
bn )k∞ + kC
bn+1 − rRn (C
bn )k∞ ≤ (r + λL δn+1 )kRn (C
M
A

2
bn )k∞ ≤ λ2 M
c
c
Since kRn (C
A n,0 + σξ , we have our bound for Mn+1,0 . Likewise,

cn+1,k+1 = max |[C
bn+1 ]i,j |λ−d(i,j)∧(k+1)
M
A
i,j

bn )]i,j |λ
≤ r max |[Rn (C
A

−d(i,j)∧(k+1)

i,j

bn+1]i,j − r[Rn (C
bn )]i,j |λ−L
+ max |[C
A
i,j

bn+1,k+1 + δn+1 kRn (C
bn )k∞
= rR
cn,k ), R
bn+1,0 } + δn+1 kRn (C
bn )k∞ ≤ ψλ (M
cn,k , δn+1 ).
≤ r max{λA φ(M
A

6.3

Stability of localized structures

Lemma 6.5. Under the conditions of Theorem 2.5, when K > Γ(rǫ−1 , d) with ǫ = min{ 2λ1A −
r δ
, }, the diagonal status defined by (6.1) satisfies:
2 2
cn+1,0 ≤ λA M
cn,0 + (r + δ)σξ2
En M

c2
c2
En M
n+1,0 ≤ λA Mn,0 +

Therefore, by Gronwall’s inequality,

2
cn,0 ≤ λn M
c
E0 M
A 0,0 + (r + δ)σξ

n
X
k=0

(r + δ)2 σξ4
1 − λA

a.s..
a.s..

c0,0 +
λkA ≤ λnA M

c2 ≤ λn M
c2
E0 M
n,0
A 0,0 +

33

(r + δ)2 σξ4
(1 − λA )2

(r + δ)σξ2
1 − λA

a.s..

a.s..

Proof. We apply Lemma 6.1, Proposition 6.3 to find that
b n H)C
bn (I − K
b n H)T + σ 2 K
b bT
c
b
k(I − K
o n Kn k∞ = Mn,0 ≤ Mn,0 = kCn k∞ ,

and by the first claim of Proposition 6.2,

2
2 b
2
bn )k∞ ≤ λ2 k(I − K
b n H)C
bn (I − K
b n H)T + σ 2 K
b bT
kRn (C
A
o n Kn k∞ + σξ ≤ λA kCn k∞ + σξ .

Also by Young’s inequality, one can show that

bn k∞ + σ 2 )2 ≤ λ3 kC
bn k2 +
bn )k2 ≤ (λ2 kC
kRn (C
∞
A
ξ
A
∞

σξ4
.
1 − λA

With ǫ = min{ 2λ1A − 2r , 2δ }, when K > Γ(rǫ−1 , d), by Corollary 5.4 b),
bn+1 − rRn (C
bn )k∞ ≤ ǫkRn (C
bn )k∞
En kC

a.s.,

bn+1 − rRn (C
bn )k2∞ ≤ 2ǫrkRn (C
bn )k2∞
En kC

a.s..

2 b
2
b
By ǫ + r ≤ λ−1
A and kRn (Cn )k∞ ≤ λA kCn k∞ + σξ ,

bn+1 k∞ ≤ En kC
bn+1 − rRn (C
bn )k∞ + rkRn (C
bn )k∞
En kC
bn )k∞ ≤ λA kC
bn k∞ + (r + δ)σ 2 .
≤ (r + ǫ)kRn (C
ξ

Likewise, because (r + 2ǫ) ≤ λ−1
A ,

bn+1k2 ≤ En kC
bn+1 − rRn (C
bn )k2 + r 2 kRn (C
bn )k2
En kC
∞
∞
∞
b
b
b
+ 2rkRn (Cn )k∞ En kCn+1 − rRn (Cn )k∞
bn )k2∞
≤ (2ǫr + r 2 + 2ǫr)kRn (C
bn )k2
≤ (r + 2ǫ)2 kRn (C
∞

(r + δ)2 σξ4
2
b
≤ λA k C n k ∞ +
.
1 − λA

Lemma 6.6. Suppose the following holds
log 4δ∗−1
,
n∗ ≥ 2L +
log λ−1
A

1
1
δ∗ ≤ , δ∗ ≤ (λ−1
− r),
4
2 A

and the sample size satisfies (2.24). Then
cn∗ ,L ≤ 1 M
c0,L + (1 + 2δ∗ )M∗
E0 M
2
34

a.s..

c0,L >
Proof. Case 1: if M

4(r+δ∗ )σξ2
λL
A (1−λA )

. By Lemma 6.1

ck,0 ≤ M
ck,L ≤ λ−L M
ck,0.
M
A

Then by Lemma 6.5

cn,L ≤ E0 λ−L M
cn,0 ≤ λn−L M
c0,0 +
E0 M
A
A

(r + δ∗ )σξ2
(r + δ∗ )σξ2
n−L c
≤
λ
M
+
0,L
A
λLA (1 − λA )
λLA (1 − λA )

a.s..

n∗ −L
By our choice of n∗ , λA
≤ 41 , so we have our claim, since

c0,L ≤
Case 2: if M

(r + δ∗ )σξ2
1c
c0,L +
cn∗ ,L ≤ 1 M
≤ M
E0 M
0,L
L
4
λA (1 − λA )
2

4(r+δ∗ )σξ2

λL
A (1−λA )

a.s..

. Consider the event

U = {δk ≤ δ∗ , ∀k ≤ n∗ }.

Denote its complementary set as U c . Then the expectation can be decomposed as
q
p
cn∗ ,L 1U + P0 (U c ) E0 M
cn∗ ,L ≤ E0 M
cn∗ ,L 1U + E0 M
cn∗ ,L 1U c ≤ E0 M
c2 ,
E0 M
n∗ ,L

where we applied the Cauchy inequality for the U c part, and P0 is the probability conditioned
on F0 . We will find a bound for each of the two parts.
If U holds, then δn+1 ≤ δ∗ for n ≤ n∗ − 1. By Proposition 6.4,
2
2
cn+1,0 ≤ (r + δ∗ )(λ2 M
c
c
M
A n,0 + σξ ) ≤ λA Mn,0 + (r + δ∗ )σξ .

Then by the Gronwall’s inequality, under U,

c0,0 ≤ M
c0,L ≤
Because M

c
cn,0 ≤ λn M
M
A 0,0 +

4(r+δ∗ )σξ2

λL
A (1−λA )

c0,0 ≤
λnA0 M

(r + δ∗ )σξ2
.
1 − λA

, so after n0 = n∗ − L ≥ L + ⌈− log(4rδ∗−1 + 4)/ log λA ⌉

δ∗ σξ2
,
1 − λA

cn0 ,0 ≤
so M

(r + 2δ∗ )σξ2
≤ M∗ .
1 − λA

In the next 1 ≤ k ≤ L steps, since δn ≤ δ∗ when U holds, because ψλA is increasing, by
Proposition 6.4
cn0 +k,k ≤ ψλ (M
cn0 +k−1,k−1, δ∗ ),
M
A

cn0 +L,L ≤ M∗ . Therefore by n∗ = n0 + L,
we can derive that M
cn∗ ,L 1U ≤ M∗ ,
E0 M
35

a.s..

In order to conclude our claim, it suffices to show that
c2 ≤ δ 2 M 2 ,
P0 (U c )E0 M
n,L
∗
∗

c0,L ≤
Apply Lemma 6.5 with δ = δ∗ , recall that M

a.s..

4(r+δ∗ )σξ2

λL
A (1−λA )

(6.3)

n∗ −2L
and 16λA
≤ 1,

(r + δ∗ )2 σξ4
(1 − λA )2
16(r + δ∗ )2 σξ4 (r + δ∗ )2 σξ4
n∗
≤ 2M∗2 .
+
≤ λA 2L
2
2
(1 − λA )
λA (1 − λA )

2
cn2 ,L ≤λ−2L E0 M
c0,L
cn2 ,0 ≤ λn∗ −2L M
E0 M
+
A
A
∗
∗

Moreover, by Theorem 2.1 b)
2
P(δn+1 > δ∗ |Fn ) ≤ 8d2 exp(−cKλ2L
A δ∗ ) ≤

δ∗2
.
2n∗

where the final bound comes with the sample K satisfying (2.24). Therefore, by the law of
iterated expectation,
c

P0 (U ) ≤

n∗
X

P0 (δk > δ∗ ) =

k=1

nX
∗ −1
k=0

1
E0 P(δk+1 > δ∗ |Fk ) ≤ δ∗2 ,
2

and (6.3) comes as a result.
cn,L . So
Proof of Theorem 2.5. Recall that Mn = M

1
E0 Mn∗ ≤ M0 + (1 + δ∗ )M∗
2

has been proved by Lemma 6.6. This leads to the following using Gronwall’s inequality,
E0 Mjn∗ ≤

1
M0 + 2(1 + δ∗ )M∗ .
2j

Next, for k = 1, · · · , n∗ − 1, apply Lemma 6.5 with δ = δ∗
ck,L ≤ λ−L EM
ck,0 ≤ λ−L EM
c0,0 +
EMk = EM
A
A

(r + δ∗ )σξ2
b
≤ λ−L
A (EkC0 k + M∗ ),
λLA (1 − λA )

c0,0 = kC
b0 k∞ ≤ kC
b0 k by Lemma 6.1. Then if k + mn∗ ≤ T ,
because M
m
X
j=0

EMk+jn∗ =

m
X
j=0

EEk Mk+jn∗

m
X
1
≤
EMk + 2(1 + δ∗ )M∗
2j
j=0

ck,L + 2(m + 1)(1 + δ∗ )M∗
≤ 2EM
b0 k + M∗ ) + 2(m + 1)(1 + δ∗ )M∗ .
≤ 2λ−L (EkC
A

Summation of the inequality above with k = 0, · · · , n∗ − 1, we obtain our final claim.
36

6.4

Small noise scaling

Proof of Theorem 2.7. It suffices to verify the conditions of Theorems 2.4 and 2.5 under the
small noise scaling.
First, we check Theorem 2.5. Condition 1) is invariant except that Σn = ǫσξ2 Id . Condition
2) concerns only of An , so it and λA are also invariant under small noise scaling. For condition
3), if it holds without small noise scaling, that is
n
o
2
−2
−2
2
2
2
(r + δ∗ ) max λA M∗ 1 + σo M∗ + λA σo M∗ , λA M∗ + σξ ≤ M∗ .
This leads to

o
n
2
(r + δ∗ ) max λA (ǫM∗ ) 1 + (ǫσo2 )−1 (ǫM∗ ) + λA (ǫσo2 )−1 (ǫM∗ )2 , λ2A ǫM∗ + ǫσξ2 ≤ ǫM∗ .

Moreover, condition 3) requires that

(r + δ∗ )σξ2
M∗ ≥
1 − λA

(r + δ∗ )ǫσξ2
ǫM∗ ≥
.
1 − λA

⇒

Therefore, with small scaling, condition 3) holds with the same δ∗ , while M∗ is replaced by
ǫM∗ . Condition 4) is invariant under the small noise scaling, since δ∗ and λA are invariant.
As a consequence, Theorem 2.5 implies the following:
T
2n∗
1X
b0 k + ǫM∗ ) + 2(1 + δ∗ )ǫM∗ .
EMk ≤
(EkC
T k=1
T λLA

(6.4)

ck,L ≥ kC
bk k∞ by Lemma 6.1.
This yields the first claimed result, since Mk = M
Next we check the conditions of Theorem 2.4. For condition 1), mΣ and MΣ need to
be replaced by ǫσξ2 since we assume Σn = ǫσξ2 Id . Condition 2) still holds with (r0 , ρ) →
(ǫ−1 r0 , ǫρ) since
b0 + ρId )
Eê0 ⊗ ê0  r0 (C

⇒

b0 + ǫρId ).
Eê0 ⊗ ê0  (ǫ−1 r0 )(C

Condition 3) is guaranteed by (6.4) above, with M0 = 2(1 + δ∗ )ǫM∗ . Condition 4) and
condition 5) are both invariant, as it concerns only geometry quantities. Finally it suffices
to plug in all the estimates for the result, and find
T −1
1X
bn ◦ DLcut + ǫρId ))
1−
P(ES ên ⊗ ên  r∗ (C
T n=0

≤

b0 k + ǫM∗ )
2δn∗ (EkC
r0
2r 1/3
+
(ρ−1 Bl2 MA2 + (ρσ
1/3 )
L
o)
T ǫ log r∗
T ǫλA log r∗

δ  −1 2 2
−1 2
2r 1/3
r 1/3 2/3
+
2(ρ Bl MA + (ρσ
.
)(1
+
δ
)M
+
ρ
σ
+
2
σ
∗
∗
1/3
ξ
ρ1/3 o
o)
log r∗

Note that in above some ǫ terms are upper-bounded by 1, so the inequality has a simpler
form.
37

7

Conclusion and discussion

Ensemble Kalman filter (EnKF) is a popular tool for high dimensional data assimilation
problems. Domain localization is an important EnKF technique that exploits the natural
localized covariance structure, and simplifies the associated sampling task. We rigorously
investigate the performance of localized EnKF (LEnKF) for linear systems. We show in
Theorem 2.4 that in order for the filter error covariance to be dominated by the ensemble
covariance, 1) the sample size K needs to exceed a constant that depends on the localization
radius and the logarithmic of the state dimension, 2) the forecast covariance has a stable
localized structure. Condition 2) is necessary for an intrinsic localization inconsistency to be
bounded. This condition is usually assumed in LEnKF operations, but it can also be verified
for systems with weak local interaction and sparse observation by Theorem 2.5.
While the results here provide the first successive explanation of LEnKF performance
with almost dimension independent sample size, there are several issues that require further
study. In below we discuss a few of them.
1. There are several ways to apply the localization technique in EnKF. We discuss here
only the domain localization with standard EnKF procedures. In principle, our results
can be generalized to the covariance localization/tempering technique, and also the
popular ensemble square root implementation. But such generalization will not be
trivial, as the Kalman gain will not be of a small bandwidth, and localization techniques
will have unclear impact on the square root SVD operation.
2. This paper studies the sampling effect of LEnKF and shows the sampling error is controllable. Yet LEnKF without sampling error, in other words, LEnKF in the large
ensemble limit, is not well studied mathematically. The effect of the localization techniques on the classical Kalman filter controllability and observability condition is not
known. This may lead to practical guidelines in the choice of localization radius.
3. Theorem 2.5 provides the first proof that LEnKF covariance has a stable localized
structure. But the conditions we impose here are quite strong, while localized structure
is taken for granted in practice. How to show it in general nonlinear settings is a very
interesting question.

Acknowledgement
This research is supported by the NUS grant R-146-000-226-133, where X.T.T. is the principal investigator. The author thanks Andrew J. Majda, Lars Nerger and Ramon van Handel
for their discussion on various parts of this paper.

References
[1] C. Snyder, T. Bengtsson, and P. J. Bickel. Obstacles to high-dimensional particle
filtering. Mon. Wea. Rev., 136(12):4629–4640, 2008.

38

[2] P. J. van Leeuwen. Particle filtering in geophysical systems. Mon. Wea. Rev., 137:4089–
4114, 2009.
[3] J. L. Anderson. An ensemble adjustment Kalman filter for data assimilation. Mon.
Weather Rev., 129(12):2884–2903, 2001.
[4] T. M. Hamill, C. Whitaker, and C. Snyder. Distance-dependent filtering of background
error covariance estimates in an ensemble Kalman filter. Mon. Weather Rev., 129:2776–
2790, 2001.
[5] G. Evensen. The ensemble Kalman filter: Theoretical formulation and practical implementation. Ocean dynamics, 53(4):343–367, 2003.
[6] A. J. Majda and J. Harlim. Filtering complex turbulent systems. Cambridge University
Press, Cambridge, UK, 2012.
[7] E. Kalnay. Atmospheric modeling, data assimilation, and predictability. Cambridge
university press, 2003.
[8] P. L. Houtekamer and H. L. Mitchell. Data assimilation using an ensemble kalman filter
technique. Mon. Wea. Rev., 126(3):796–811, 1998.
[9] J. S. Whitaker and T. M. Hamill. Ensemble data assimilation without perturbed observations. Mon. Wea. Rev., 130(7):1913–1924, 2002.
[10] T. Miyoshi and S. Yamane. Local ensemble transform Kalman filtering with an AGCM
at a T159/L48 resolution. Mon. Wea. Rev., 135(11):3841–3861, 2007.
[11] B. R. Hunt, E. J. Kostelich, and I. Szunyogh. Efficient data assimilation for spatiotemporal chaos: a local ensemble transform Kalman filter. Physica D, 230(1):112–126,
2007.
[12] K. Bergemann and S. Reich. A localization technique for ensemble Kalman filters.
Quart. J. Roy. Meteor. Soc., 136(648):701–707, 2010.
[13] T. Janjić, L. Nerger, A. Albertlla, J. Schröter, and S. Skachoko. On domain localization
in ensemble-based Kalman filter algorithms. Mon. Wea. Rev., 139(7):2046–2060, 2011.
[14] L. Nerger, T. Janjić, J. Schröter, and W. Hiller. A regulated localization scheme for
ensemble-based Kalman filters. Quart. J. Roy. Meteor. Soc., 138:802–812, 2012.
[15] L. Nerger. On serial observation processing in localized ensemble Kalman filters. Mon.
Wea. Rev., 143(5):1554–1567, 2015.
[16] H. R. Künsch and Sylvian Robert. Localizing the ensemble Kalman particle filter. Tellus
A: Dynamic Meteorology and Oceanography, 69(1):1282016, 2017.
[17] M. D. L. Chevrotiere and J. Harlim. A data-driven method for improving the correlation
estimation in serial ensemble Kalman filters. Mon. Wea. Rev., 145(3):985–1001, 2017.
39

[18] P. Rebeschini and R. Van Handel. Can local particle filters beat the curse of dimensionality? Ann. Appl. Probab., 25:2809–2866, 2015.
[19] A. J. Majda and Y. Lee. State estimation and prediction using clustered particle filters.
Proc. Natl. Acad. Sci., 113(51):14609–14614, 2016.
[20] J. Poterjoy. A localized particle filter for high-dimensional nonlinear systems. Mon.
Wea. Rev., 144(1):59–76, 2016.
[21] P. J. Bickel and E. Levina. Regularized estimation of large covariance matrices. The
Annals of Statistics, 36(1):199–227, 2008.
[22] D. Sanz-Alonso and A. M. Stuart. Long time asymptotics of the filtering distribution
for partially observed chaotic dynamical systems. SIAM/ASA J. Uncertainty Quantification, 3:1200–1220, 2015.
[23] A. J. Majda and X. T. Tong. Performance of ensemble Kalman filters in large dimensions. Accepted by Comm. Pure Appl. Math. arXiv: 1606.09321, 2017.
[24] D. T. B. Kelly and A. M. Stuart. Ergodicity and accuracy of optimal particle filters for
Bayesian data assimilation. arXiv:1611.08761, 2017.
[25] J. S. Whitaker, T. M. Hamill, X. Wei, Y. Song, and Z. Toth. Ensemble data assimilation
with the NCEP global forecast system. Mon. Wea. Rev., 136(2):463–482, 2008.
[26] A. N. Bishop, P. Del Moral, and S. D. Pathiraja. Perturbations and projections of
Kalman-Bucy semigroups motivated by methods in data assimilation. arXiv:1701.05978,
2017.
[27] A. J. Majda and X. T. Tong. Rigorous accuracy and robustness analysis for two-scale
reduced random Kalman filters in high dimensions. arXiv: 1606.09087, 2016.
[28] R. S. Liptser and A. N. Shiryaev. Statistics of random processes. I, II,, volume 5 of
Applications of Mathematics. Springer-Verlag, 2001.
[29] J. Mandel. Efficient implementation of the ensemble Kalman filter. Technical Report
UCDHSC/CCM Report No. 231, University of Colorado at Denver and Health Sciences
Center, 2006.
[30] G. Burgers, P. J. van Leeuwen, and G. Evensen. Analysis scheme in the ensemble
Kalman filter. Mon. Wea. Rev., 126(6):1719–1724, 1998.
[31] G. Gaspari and S. E. Cohn. Construction of correlation functions in two and three
dimensions. Quarterly journal of the Royal Meteorological Society, 125(554):723–757,
1999.
[32] Lorenc. A. C. The potential of the ensemble Kalman filter for NWP-a comparison with
4D-Var. Quart. J. Roy. Meteor. Soc., 129(595):3183–3203, 2003.

40

[33] P. J. Bickel and M. Lindner. Approximating the inverse of banded matrices by banded
matrices with applications to probability and statistics. Theory of Probaility & its
Applications, 56(1):1–20, 2012.
[34] R. Furrer and T. Bengtsson. Estimation of hign-dimensional prior and posterior covariance matrices in Kalman filter variants. Journal of Multivariate Analysis, 98:227–255,
2007.
[35] Li. H., E. Kalnay, and T. Miyoshi. Simultaneous estimation of covariance inflation
and observation errors within an ensemble kalman filter. Quart. J. Roy. Meteor. Soc.,
135:523–533, 2009.
[36] R. Vershynin. Introduction to the non-asymptotic analysis of random matrices. In
Y Eldar and G Kutyniok, editors, Compressed Sensing, Theory and Applications, pages
210–268. Cambridge University Press, 2011.
[37] J. Mandel, L. Cobb, and J. D. Beezley. On the convergence of the ensemble Kalman
filter. Applications of Mathematics, 56(6):533–541, 2011.
[38] E. Kwiatkowski and J. Mandel. Convergence of the square root ensemble Kalman filter
in the large ensemble limit. SIAM/ASA J. Uncertainty Quantification, 3(1):1–17, 2015.
[39] K. J. Law, H. Tembine, and R. Tempone. Deterministic mean-field ensemble Kalman
filtering. SIAM J. Scientific Computing, 38(3):A1251–A1279, 2016.
[40] D. T. B. Kelly, K. J. Law, and A. M. Stuart. Well-posedness and accuracy of the
ensemble Kalman filter in discrete and continuous time. Nonlinearity, 27:2579–2603,
2014.
[41] X. T. Tong, A. J. Majda, and D. Kelly. Nonlinear stability and ergodicity of ensemble
based Kalman filters. Nonlinearity, 29:657–691, 2016.
[42] X. T. Tong, A. J. Majda, and D. Kelly. Nonlinear stability of the ensemble Kalman
filter with adaptive covariance inflation. Commun. Math. Sci., 14(5):1283–1313, 2016.
[43] D. Kelly, A. J. Majda, and X. T. Tong. Concrete ensemble Kalman filters with rigorous
catastrophic filter divergence. Proc. Natl. Acad. Sci., 112(34):10589–10594, 2016.
[44] M. Rudelson and R. Vershynin. Hanson-Wright inequality and sub-gaussian concentration. Electron. Commun. Probab., 18(82):1–9, 2013.

41

