arXiv:1702.08061v2 [stat.ME] 28 Feb 2017

The Ensemble Kalman Filter:
A Signal Processing Perspective
Michael Roth, Gustaf Hendeby, Carsten Fritsche,
and Fredrik Gustafsson
Dept. of Electrical Engineering,
Linköping University,
Linköping, Sweden
michael.roth@liu.se
March 1, 2017
Abstract
The ensemble Kalman filter (EnKF) is a Monte Carlo based implementation of the Kalman filter (KF) for extremely high-dimensional, possibly nonlinear and non-Gaussian state estimation problems. Its ability
to handle state dimensions in the order of millions has made the EnKF
a popular algorithm in different geoscientific disciplines. Despite a similarly vital need for scalable algorithms in signal processing, e.g., to make
sense of the ever increasing amount of sensor data, the EnKF is hardly
discussed in our field.
This self-contained review paper is aimed at signal processing researchers and provides all the knowledge to get started with the EnKF.
The algorithm is derived in a KF framework, without the often encountered geoscientific terminology. Algorithmic challenges and required extensions of the EnKF are provided, as well as relations to sigma-point
KF and particle filters. The relevant EnKF literature is summarized in
an extensive survey and unique simulation examples, including popular
benchmark problems, complement the theory with practical insights. The
signal processing perspective highlights new directions of research and facilitates the exchange of potentially beneficial ideas, both for the EnKF
and high-dimensional nonlinear and non-Gaussian filtering in general.

1

Introduction

Numerical weather prediction [1] is an extremely high-dimensional geoscientific
state estimation problem. The state x comprises physical quantities (temperature, wind speed, air pressure, etc.) at many spatially distributed grid points,
which often yields a state dimension n in the order of millions. Consequently,
the Kalman filter (KF) [2, 3] or its nonlinear extensions [4, 5] that require the
1

storage and processing of n × n covariance matrices cannot be applied directly.
It is well-known that the application of particle filters [6, 7] is not feasible either.
In contrast, the ensemble Kalman filter [8, 9] (EnKF) was specifically developed
as algorithm for high-dimensional n. The EnKF
• is a random-sampling implementation of the KF;
• reduces the computational complexity of the KF by propagating an ensemble of N < n state realizations;
• can be applied to nonlinear state-space models without the need to compute Jacobian matrices;
• can be applied to continuous-time as well as discrete-time state transition
functions;
• can be applied to non-Gaussian noise densities;
• is simple to implement;
• does not converge to the Bayesian filtering solution for N → ∞ in general;
• often requires extra measures to work in practice.
Also in the field of stochastic signal processing (SP) and Bayesian state estimation, high-dimensional problems become more and more relevant. Examples
include SLAM [10] where x contains an increasing number of landmark positions, or extended target tracking [11, 12] where x can contain many parameters
to describe the shape of the target. Furthermore, scalable SP algorithms are
required to make sense of the ever increasing amount of data from sensors in
everyday devices. EnKF approaches hardly appear in the relevant SP journals,
though. In contrast, vivid theoretical development is documented in geoscientific journals under the umbrella term data assimilation (DA) [1]. Hence, a
relevant SP problem is being addressed with only little participation from the
SP community. Conversely, much of the DA literature makes little reference to
relevant SP contributions. It is our intention to bridge this interesting gap.
There is further overlap that motivates for a closer investigation of the EnKF.
First, the basic EnKF [9] can be applied to nonlinear and non-Gaussian statespace models because it is entirely sampling based. In fact, the state evolution in
geoscientific applications is typically governed by large nonlinear black box prediction models derived from partial differential equations. Furthermore, satellite
measurements in weather applications are nonlinearly related to the states [1].
Hence, the EnKF has long been investigated as nonlinear filter. Second, the
EnKF literature contains so called localization methods [13, 14] to systematically approach high-dimensional problems by only acting on a part of the state
vector in each measurement update. These ideas can be directly transferred to
sigma point filters [5]. Third, the EnKF offers several interesting opportunities
to apply SP techniques, e.g., via the application of bootstrap or regularization
methods in the EnKF gain computation.
2

The contributions of this paper aim at making the EnKF more accessible
to SP researchers. We provide a concise derivation of the EnKF based on the
KF. A literature review highlights important EnKF papers with their respective contributions, and facilitates an easier access to the extensive and rapidly
developing DA literature on the EnKF. Moreover, we put the EnKF in context
with popular SP algorithms such as sigma point filters [4, 5] and the particle
filter [6, 7]. Our presentation forms a solid basis for further developments and
the transfer of beneficial ideas and techniques between the fields of SP and DA.
The structure of the paper is as follows. After an extensive literature review
in Sec. 2, the EnKF is developed from the KF in Sec. 3. Algorithmic properties
and challenges of the EnKF and the available approaches to face them are
discussed in Sec. 4 and 5, respectively. Relations to other filtering algorithms
are discussed in Sec. 6. The theoretical considerations are followed by numerical
simulations in Sec. 7 and some concluding remarks in Sec. 8.

2

Filtering and EnKF literature

The following literature review provides important landmarks for the EnKF
novice.
State-space models and the filtering problem have been investigated since
the 1960s. Early results include the Kalman filter (KF) [2] as algorithm for
linear systems, and the Bayesian filtering equations [15] as theoretical solution
for nonlinear and non-Gaussian systems. Because the latter approach cannot be
implemented in general, approximate filtering algorithms are required. With a
leap in computing capacity, the 1990s saw major developments. The samplingbased sigma point Kalman filters [4, 5] started to appear. Furthermore, particle
filters [6, 7] were developed to approximately implement the Bayesian filtering
equations through sequential importance sampling.
The first EnKF [8] was proposed in a geoscientific journal in 1994 and introduced the idea of propagating ensembles to mimic the KF. A systematic error
that resulted in an underestimated uncertainty was later corrected by processing “perturbed measurements”. This randomization is well motivated in [9] but
also used in [13].
Interestingly, [8] remains the most cited EnKF paper1 , followed by the
overview article [16] and the monograph [17] by the same author. Other insightful overviews from a geoscientific perspective are [18, 19]. Many practical
aspects of operational EnKF for weather prediction and re-analysis are described
in [19–21]. Whereas the aforementioned papers were mostly published in geoscientific outlets, a special issue of the IEEE Control Systems Magazine appeared
with review articles [22–24] and an EnKF case study [25]. Still, the above material was written by EnKF researchers with a geoscientific focus and in the
application specific terminology. Furthermore, references to the recent SP literature and other nonlinear KF variants [5] are scarce.
1 With

over 3000 citations between 1994 and 2016.

3

Some attention has been devoted to the EnKF also beyond the geosciences.
Convergence properties for N → ∞ have been established using different theoretical analyses of the EnKF [26–28]. Statistical perspectives are provided in
the thesis [29] and the review [30]. A recommended SP view that connects
the EnKF with Bayesian filtering and particle methods, including convergence
results for nonlinear systems, is [31]. Examples of the EnKF as tool for tomographic imaging and target tracking are described in [32] and [33], respectively.
Brief introductory papers that connect the EnKF with more established SP
algorithms include [34] and [35]. The latter also served as basis for this article.
The majority of EnKF advances is still documented in geoscientific publications. Notable contributions include deterministic EnKF that avoid the randomization of [9] and propagate an ensemble of deviations from the ensemble
mean [16, 36–38]. Their common basis as square root EnKF and the relation
to square root KF [3] is discussed in [39]. The computational advantages in
high-dimensional EnKF with small ensembles (N  n) come at the price of
adverse effects, including the risk of filter divergence. The often encountered
underestimation of uncertainty can be counteracted with ensemble inflation [40].
A scheme with two EnKF in parallel that provide each other with gain matrices
to reduce unwanted “inbreeding” has been suggested in [13]. The benefit of such
a double EnKF is, however, debated [38, 41]. The low-rank approximation of
covariance matrices can yield spurious correlations between supposedly uncorrelated state components and measurements. Localization techniques such as
local measurement updates [13, 16, 42] or covariance tapering [14, 43] let the
measurement only affect a part of the state vector. In other words, localization
effectively reduces the dimension of each measurement update. Inflation and
localization are essential components of operational EnKF [19]. A list of further
advances in the geoscientific literature is provided in the appendix of [17].
An interesting development for SP researchers is the reconsideration of particle filters (PF) for high-dimensional geoscientific problems, with seemingly little
reference to SP literature. An early example is [44]. The well-known challenges,
mostly related to the problem of importance sampling in high dimensions, are
reviewed in [45, 46]. Several recent approaches [47–49] were successfully tested
on a popular EnKF benchmark problem [50] that is also investigated in the
simulation examples of this paper. Combinations of the EnKF with the deterministic sampling of sigma point filters [5] are given in [51] and [52]. However,
the benefit of the unscented transformation [5, 53] in [52] is debated in [54]. Ideas
to combine the EnKF with Gaussian mixture approaches are given in [55–57].

3

A Signal Processing Introduction to the Ensemble Kalman Filter

The underlying framework of our filter presentation are discrete-time state-space
models [3, 15]. The Kalman filter and many EnKF variants are built upon the

4

linear model
xk+1 = F xk + Gvk ,

(1a)

yk = Hxk + ek ,

(1b)

with the n-dimensional state x and the m-dimensional measurement y. The
initial state x0 , the process noise vk , and the measurement noise ek are described by E(x0 ) = x̂0 , E(vk ) = 0, E(ek ) = 0, cov(x0 ) = P0 , cov(vk ) = Q, and
cov(ek ) = R. In the Gaussian case, these moments completely characterize the
distributions of x0 , vk , and ek .
Nonlinear relations in the state evolution and measurement equations can
be described by a more general model
xk+1 = f (xk , vk ),

(2a)

yk = h(xk , ek ).

(2b)

More general noise and initial state distributions can, for example, be characterized by probability density functions p(x0 ), p(vk ), and p(ek ).
Both (1) and (2) can be time-varying but the time indices for functions and
matrices are omitted for convenience.

3.1

A brief Kalman filter review

The KF is an optimal linear filter [3] for (1) that propagates state estimates
x̂k|k and covariance matrices Pk|k .
The KF time update or prediction is given by
x̂k+1|k = F x̂k|k ,

(3a)
T

T

Pk+1|k = F Pk|k F + GQG .

(3b)

The above parameters can be used to predict the output of (1) and its uncertainty via
ŷk|k−1 = H x̂k|k−1 ,

(4a)
T

Sk = HPk|k−1 H + R.

(4b)

The measurement update adjusts the prediction results according to
x̂k|k = x̂k|k−1 + Kk (yk − ŷk|k−1 )
= (I − Kk H)x̂k|k−1 + Kk yk ,
Pk|k = (I − Kk H)Pk|k−1 (I − Kk H)T + Kk RKkT ,

(5a)
(5b)
(5c)

with a gain matrix Kk . Here, (5b) resembles a deterministic observer and only
requires all eigenvalues of (I − Kk H) inside the unit circle to obtain a stable
filter. The optimal Kk in the minimum variance sense is given by
Kk = Pk|k−1 H T Sk−1 = Mk Sk−1 ,
5

(6)

where Mk is the cross-covariance between the state and output predictions.
Alternatives to the covariance update (5c) exist, but the shown Joseph form [3]
will simplify the derivation of the EnKF. Furthermore, it is valid for all gain
matrices Kk beyond (6) and numerically well-behaved. It should be noted that
it is numerically advisable to obtain Kk by solving Kk Sk = Mk rather than
explicitly computing Sk−1 [58].

3.2

The ensemble idea

The central idea of the EnKF is to propagate an ensemble of N < n (often
(i)
N  n) state realizations {xk }N
i=1 instead of the n-dimensional estimate x̂k|k
and the n × n covariance Pk|k of the KF. The ensemble is processed such that
x̄k|k =
P̄k|k =

XN

(i)
x ≈ x̂k|k ,
i=1 k
XN

(i)
1
xk − x̄k|k
N −1
i=1
1
N

(7a)
(i)

xk − x̄k|k

T

≈ Pk|k .

(7b)

Reduced computational complexity is achieved because the explicit computation
of P̄k|k is avoided in the EnKF recursion. Of course, this reduction comes at
the price of a low rank approximation in (7b) that entails some negative effects
and requires extra measures.
For our development it is convenient to treat the ensemble as an n × N
(i)
matrix Xk|k with columns xk . This allows for the compact notation of the
ensemble mean and covariance
x̄k|k =
P̄k|k =

1
N Xk|k 1,
1 e
eT
N −1 Xk|k Xk|k ,

(8a)
(8b)

where 1 = [1, . . . , 1]T is an N -dimensional vector and
ek|k = Xk|k − x̄k|k 1T = Xk|k (IN −
X

T
1
N 11 )

(9)

is an ensemble of deviations from x̄k|k , sometimes called ensemble anomalies [17].
The matrix multiplication in (9) provides a compact way to write the anomalies,
but is not the most efficient way to compute them.

3.3

The EnKF time update

The EnKF time update is referred to as forecast in the geoscientific literature.
In analogy to (3), a prediction ensemble Xk+1|k is computed that carries the
information in x̂k+1|k and Pk+1|k . An ensemble of N independent process noise
(i)
realizations {vk }N
i=1 with zero mean and covariance Q, stored as matrix Vk , is
used in
Xk+1|k = F Xk|k + GVk .
(10)
An extension to nonlinear state transitions (2a) is given by
Xk+1|k = f (Xk|k , Vk ),
6

(11)

where we generalized f to act on the columns of its input matrices. Apparently, the EnKF time update amounts to a one-step-ahead simulation of Xk|k .
Consequently, also continuous-time dynamics can be considered by, for example, numerically solving partial differential equations to obtain Xk+1|k . Also
non-Gaussian initial state and process noise distributions with arbitrary densities p(x0 ) and p(vk ) can be employed as long as they allow sampling. Perhaps
because of this flexibility, the time update is often omitted in the EnKF literature [9, 13].

3.4

The EnKF measurement update

The EnKF measurement update is referred to as analysis in the geoscientific
literature. A prediction or forecast ensemble Xk|k−1 is processed to obtain the
filtering ensemble Xk|k that encodes the KF mean and covariance. We assume
that a gain K̄k = Kk is given and postpone its computation to the next section.
With K̄k available, the KF update (5b) can be applied to each ensemble
member as follows [8]
Xk|k = (I − K̄k H)Xk|k−1 + K̄k yk 1T .

(12)

The resulting ensemble average (8a) approximates the correct x̂k|k . However,
with yk 1T known, the sample covariance (8b) gives only the first term of (5c)
and therefore fails to resemble Pk|k . A solution [9] is to account for the missing
term K̄k RK̄kT by adding artificial zero-mean measurement noise realizations
(i)
{ek }N
i=1 with covariance R, stored as matrix Ek , according to
Xk|k = (I − K̄k H)Xk|k−1 + K̄k yk 1T − K̄k Ek .

(13)

Then, Xk|k correctly resembles x̂k|k and Pk|k . The model (1) is implicit in (13)
because the matrix H appears. If we, in analogy to (4), define a predicted
output ensemble
Yk|k−1 = HXk|k−1 + Ek
(14)
that encodes ŷk|k−1 and Sk , we can reformulate (13) to an update that resembles
(5a):
Xk|k = Xk|k−1 + K̄k (yk 1T − Yk|k−1 ).
(15)
In contrast to (13), the update (15) is entirely sampling based. As a consequence,
we can extend the algorithm to nonlinear measurement models (2b) by replacing
(14) with
Yk|k−1 = h(Xk|k−1 , Ek ),
(16)
where we generalized h to accept matrix inputs similar to (11).
In the EnKF literature, the prevailing view of inserting artificial noise is that
perturbed measurements yk 1T − Ek are processed. This might appear unusual
from an SP perspective since it suggests that information is distorted before
processing. The introduction of output ensembles Yk|k−1 , in contrast, yields a
direct connection to (4) and highlights the similarities between (15) and (5a).
7

An interesting point [57] is that the measurement enters linearly in (13)
and (15) and merely shifts the ensemble locations. This highlights the EnKF
roots in the linear KF in which Pk|k also remains unchanged by yk .

3.5

The EnKF gain

The optimal gain (6) in the KF is computed from the covariance matrices of
the predicted state and output. In the EnKF, the required Mk and Sk are not
available but must be approximated from the prediction ensembles (10) or (11),
and (14) or (16).
A straightforward way to compute the EnKF gain K̄k is to first compute
the deviations or anomalies
ek|k−1 = Xk|k−1 (IN − 1 11T ),
X
N
Yek|k−1 = Yk|k−1 (IN − N1 11T ),

(17a)
(17b)

and second the sample covariance matrices
M̄k =

1 e
eT
N −1 Xk|k−1 Yk|k−1 ,

(17c)

S̄k =

1 e
eT
N −1 Yk|k−1 Yk|k−1 .

(17d)

The computations (17) are entirely sampling-based, which is useful for the nonlinear case but introduces extra sampling errors. An obvious improvement for
additive measurement noise ek with covariance R is given in Sec. 5.2, together
with the square root EnKF that avoid the insertion of Ek altogether.
Similar to the KF, the gain K̄k should be obtained from the solution of a
linear system of equations
T
ek|k−1 Ye T
K̄k Yek|k−1 Yek|k−1
=X
k|k−1 .

4

(18)

Some Properties and Challenges of the EnKF

After a brief review of convergence results and the computational complexity
of the EnKF, we discuss adverse effects that can occur in EnKF with finite
ensemble size N .

4.1

Asymptotic convergence results

In linear systems the EnKF mean and covariance (7) converge to the KF results (5) as N → ∞. This result has been established from different theoretical
perspectives [26–28, 31].
For nonlinear systems the convergence is not as tangible. An investigation of
the EnKF as particle system is given in [31], with the outcome that the EnKF
does not give the Bayesian filtering solution except for the linear Gaussian case.
An illustration of this property is given in the example of Sec. 7.2.

8

4.2

Computational complexity

For the complexity analysis we assume that we are only interested in the filtering
results and that n > N > m, that is, the number of measurements is less than
the ensemble size and state dimension.
The KF propagates the n-dimensional mean vector x̂k|k and the n×n covariance matrix Pk|k with n(n + 1)/2 unique entries. These storage requirements of
O(n2 /2) dominate for large n > m. The EnKF requires the storage of only nN
values. The space required to store the Kalman gain and other intermediate
results is similar for the KF and EnKF. A reduction via sequential processing
of measurements, as explained in Sec. 5.1, is possible for both.
For large n the computational bottleneck of the KF is the covariance time
update (3b). Without considering any potential structure in F , slightly less
than O(n3 ) floating point operations (flops) are required. Contemporary matrix
multiplication routines [58] achieve a reduction to roughly O(n2.4 ). The EnKF
time update requires the propagation of N realizations. If each propagation
costs O(n2 ) flops, then time update is achieved in O(n2 N ) flops.
The computation of the KF gain requires O(n2 m) flops for the computation
of Mk and Sk . The solution of (6) for Kk amounts to O(m3 ). The actual
measurement update (5) adds further O(n2 m) flops. For large n, the total cost
is O(n2 m). In contrast, the EnKF parameters M̄k and S̄k can be computed
in O(nmN ) flops which, again, dominates the total cost of the measurement
update for large n. So, the EnKF reduces the flop count by a factor N
n.

4.3

Sampling and coupling effects for finite ensemble size

A serious issue in the EnKF is a commonly noted tendency to underestimate
the state uncertainty when using N < n ensemble members [13, 18, 19]. In
other words, the EnKF becomes over-confident and is likely to diverge [3] for
too small N . A number of causes and related effects can be noted.
First, an ensemble Xk|k−1 with too few members might not cover the relevant regions of the state-space well enough after the time update (10). The
underestimated spread persists in the measurement update (13) or (15) and also
Xk|k shows too little spread.
Second, the ensemble can only transport limited information and provide a
sampling covariance P̄k|k , (7b) or (8b), of at most rank N − 1. Consequently,
identically zero entries of Pk|k are difficult to reproduce and unwanted spurious
correlations show up in P̄k|k . An example would be an unreasonably large
correlation between the temperature at two distant locations on the globe. Of
course, these correlations also affect M̄k and S̄k , and thus the EnKF gain K̄k
in (18). As a result, state components that are actually uncorrelated to yk are
erroneously updated in (13) or (15). Again, this leads to a reduction in ensemble
spread.
Third, the ensemble members are nonlinearly coupled because the gain (18)
is computed from the ensemble. This “inbreeding” [13] increases with each
measurement update. An interesting side effect is that the ensemble is not

9

independent and Gaussian, even for linear Gaussian problems. To illustrate
this, we combine (18) and (15) to obtain
−1
e
eT
ek|k−1 Ye T
Xk|k = Xk|k−1 + (X
(yk 1T − Yk|k−1 )
k|k−1 )(Yk|k−1 Yk|k−1 )

(19)

and consider a linear model (1) with n = 1, H = 1, and a zero-mean Xk|k−1 .
Then, one member of Xk|k is given by
(j)
2
j=1 (xk|k−1 )
(y
PN
(j)
(j) 2 k
j=1 (xk|k−1 + ek )

PN

(i)
xk|k

=

(i)
xk|k−1

+

(i)

(i)

− xk|k−1 − ek|k−1 ),

(20)
(i)

which clearly shows the nonlinear dependencies that impede Gaussianity of xk|k .
Although similar conclusions hold for the general case, concise effects on the
ensemble spread are difficult to analyze. Some special cases (n = 1 and n = m,
H = I, R ∝ I) are investigated in [26] and shown to produce an underestimated
P̄k|k .
Finally, the random sampling in the measurement update by inserting measurement noise in (14) or (16) adds to the EnKF error budget. The inherent
sampling errors can be reduced by using the square root EnKF of Sec. 5.2.
Experiments suggest that there is a threshold for N above which the EnKF
works. A good example is given in [42]. Sec. 5 discusses methods such as
inflation and localization that can reduce this minimum N .

5

Important Extensions to the EnKF

The previous section highlighted some of the challenges of the EnKF. Here,
we summarize the important extensions that are often essential to achieve a
working EnKF with only few ensemble members.

5.1

Sequential updates

For the KF it is algebraically equivalent to carry out m measurement updates (5) with the scalar components of yk instead of a batch update with
the m-dimensional yk , if the measurement noise covariance R is diagonal [3].
Although often treated as a side note only, this technique is very useful. It
yields a more flexible algorithm with regard to the availability of measurements
at each time step k and reduces the computational complexity. After all, the
Kalman gain (6) merely requires a scalar division for each component of yk . An
extension to block-diagonal R is imminent.
Motivated by the large number of measurements in geoscientific problems,
sequential updates have also been suggested for the EnKF [14]. Because of the
randomness inherent to the EnKF, there is no algebraic equivalence between
sequential and batch updates. Hence, the order in which measurements are
processed has an effect on the filtering results.

10

Furthermore, an unusual alternative interpretation of sequential updates can
be found in the EnKF literature. Namely, measurement updates are carried out
“grid point by grid point” [13, 16, 42], that is, an iteration is carried out over
state rather than measurement components. We will return to this aspect in
Sec. 5.4.

5.2

Model knowledge in the EnKF and square-root filters

The sampling based derivation of the EnKF in equations (10) through (18)
facilitates a compact presentation. However, the randomization through Ek
in (14) or (16) adds Monte Carlo sampling errors to the EnKF budget. This
section discusses how these errors can be reduced for linear systems (1). Similar
results for nonlinear systems with additive noise follow easily. The interpretation
of ensembles as (rectangular) matrix square roots is a common theme in the
ek|k can be seen as an n × N
following approaches. In (8b), for instance, √N1−1 X
square root of P̄k|k .
A first thing to note is that the cross covariance Mk in the KF and its
ensemble equivalent M̄k should not be influenced by additive measurement noise
ek . Therefore, it is reasonable to replace Yek|k−1 with
ek|k−1
Zek|k−1 = H X

(21a)

so as to reduce the Monte Carlo variance of (17) using
M̄k =
S̄k =

1 e
eT
N −1 Xk|k−1 Zk|k−1 ,

(21b)

1 e
eT
N −1 Zk|k−1 Zk|k−1

(21c)

+ R.

The Kalman gain K̄k is then computed as in the KF (6). Alternatively, a matrix
1
1
T
square-root R 2 with R 2 R 2 = R can be used to factorize
#
"
h
i √1 Z
eT
1
1
k|k−1
e
N −1
S̄k = √N −1 Zk|k−1 R 2
.
(22)
T
R2
A QR decomposition [58] of the right matrix then yields a triangular m × m
square root of S̄k , and the computation of K̄k simplifies to forward and backward
substitution. Such ideas have their origin in sigma point KF variants [59].
The KF permits offline computation of the covariance matrices Pk|k for all
k because they do not depend on the measurements. In an EnKF for a linear
system (1) we can mimic this behavior by propagating zero-mean ensembles
ek|k that only carry the information of Pk|k . This is the central idea of different
X
square root EnKF [39] which were suggested in [36–38]. The name stems from
1
2
a relation
toT square root KF [3] which propagate n × n matrix square roots Pk|k
1
2
2
with Pk|k Pk|k = Pk|k . Most importantly, the artificial measurement noise and
the inherent sampling error can be avoided.

11

The following1 derivation [39] rewrites an alternative expression for (5c) using
2
ek|k−1 :
a square root Pk|k−1
and its ensemble approximation N 1−1 X
Pk|k = (I − Kk H)Pk|k−1
1

T

T

1

2
2
2
2
= Pk|k−1
(I − Pk|k−1
H T Sk−1 HPk|k−1
)Pk|k−1

≈

1 e
N −1 Xk|k−1 (I

−

(23a)

−1 e
1 eT
eT
N −1 Zk|k−1 S̄k Zk|k−1 )Xk|k−1 ,

where (21a) was used. The next step is to factorize
(I −

−1 e
1 eT
N −1 Zk|k−1 S̄k Zk|k−1 )

1

T

= Πk2 Πk2 ,

(23b)

which requires the left hand side to be positive definite. This property is easily
established for the positive definite S̄k of (21c) after realizing that the left hand
side of (23b) is a Schur complement [58] of a positive definite matrix.
1

Finally, the N × N matrix Πk2 can be used to create a deviation ensemble
1

ek|k = X
ek|k−1 Π 2
X
k

(24)

that correctly encodes Pk|k without using any random perturbations. Other
variants update the deviation ensemble via a multiplication from the left [36],
1

which is more costly for large n. Some more conditions on Πk2 must be met for
ek|k to remain zero-mean [60, 61].
X
The actual filtering is achieved by updating a single estimate according to
x̄k|k = (I − K̄k H)F x̄k−1|k−1 + K̄k yk ,

(25)

where K̄k is computed from the deviation ensembles.
There are indications that in nonlinear and non-Gaussian systems the sampling based EnKF variants should be preferable over their square root counterparts: A low-dimensional example is studied in [62]; the impression is confirmed
for a high-dimensional problem in [63].

5.3

Ensemble inflation

Ensemble or covariance inflation is a measure to counteract the tendency of
the EnKF to underestimate the state uncertainty for small N , and an important ingredient in operational EnKF [18]. The spread of the prediction ensemble Xk|k−1 is increased according to
ek|k−1
Xk|k−1 = x̄k|k−1 1T + cX

(26)

with a factor c > 1. In the EnKF context, this heuristic has been proposed
in [40]. Related concepts are dithering in the PF [7] and the “fudge factor”
to increase Pk|k−1 in the KF [64]. Extensions to adaptive inflation, where c is
adjusted online, are discussed in [23].
12

5.4

Localization

Localization is a technique to address the issue of spurious correlations in the
EnKF, and a crucial feature of operational EnKF [18, 19]. The underlying idea
applies equally well to the EnKF and the KF, and can be used to systematically
update only a part of the state vector with each measurement.
In order to explain the concept, we regard the KF measurement update for
a linear system (1) with a low-dimensional2 measurement yk . Let x = xk|k−1
and P = Pk|k−1 for notational convenience. It is possible to permute the state
components such that
 


P1 P12
0
x1


T
T
P2 P23
x = x2  , H = H1 0 0 , P = P12
.
(27)
T
x3
0 P23 P3
Only the part x1 appears in the measurement equation (1b) yk = H1 x1 + ek .
While x2 is correlated to x1 , there is zero correlation between x1 and x3 . As a
consequence, many submatrices of P vanish in the computation of

T
P H T = H1 P1 H1 P12 0 ,
(28a)
HP H T = H1 P1 H1T ,
and do not contribute to the Kalman gain (6)


P1 H1T
T
H1T  (H1 P1 H1T + R)−1 .
Kk = P12
0

(28b)

(28c)

A KF measurement update (5) with the above Kk does not affect the x3 estimate
or covariance. Hence, there is a lower-dimensional measurement update that
only alters the statistics of x1 and x2 .
Localization in the EnKF enforces the above structure using two prevailing
techniques, local updates [13, 16, 42] and covariance tapering [14, 43]. Both
rely on prior knowledge of the covariance structure. For example, the state
components are often connected to geographic locations in geoscientific applications. From the underlying physics it is reasonable to assume zero correlation
between distant states. Unfortunately, this viewpoint is not transferable to
high-dimensional problems in general.
Local updates were introduced for the sampling based EnKF in [13] and for
different square root EnKF in [16, 42]. Nonlinear measurement functions (2b)
are linearized in the latter two. All of the above references update the state
vector “grid point by grid point”, which appears unusual from a KF perspective [3]. In an iteration, local state vectors of small dimension (< N ) are chosen
and updated with a subset of supposedly relevant measurements. These “full
rank” updates avoid some of the problems associated with small N and large n.
However, discontinuities between state components are introduced [65]. Some
2 We

assume that the components can be processed sequentially.

13

heuristics to combine the local ensembles and further implementation details
can be found in [42, 66].
Under the assumption of the structure in (27), a local analysis would amount
to an EnKF update of the x1 - and x2 -components only, to avoid errors in x3 .
Covariance tapering was introduced in [13]. It contradicts the EnKF idea in
the sense that the ensemble covariance P̄k|k−1 of Xk|k−1 is processed. However,
it will become clear that not all entries of P̄k|k−1 must be computed. Prior
knowledge of a covariance structure as in (27) is used to create an n × n matrix
ρ with entries in [0, 1], and a tapered covariance (ρ ◦ P̄k|k−1 ) is computed. Here,
◦ denotes the element-wise Hadamard or Schur product [58]. A typical ρ has
ones on the diagonal and decays smoothly to zero for unwanted off-diagonal
elements. The standard choice uses a compactly supported correlation function
from [67] and is discussed in [14, 43, 65]. Subsequently, the Kalman gain is
computed as in the KF (6) using
M̄k = (ρ ◦ P̄k|k−1 )H T ,

(29a)
T

S̄k = H(ρ ◦ P̄k|k−1 )H + R,

(29b)

where we assumed a linear measurement relation (1b).
There are some technicalities associated with the tapering operation. Only
positive semi-definite ρ guarantee that (ρ ◦ P̄k|k−1 ) is a valid covariance [26].
Full rank ρ yield an increased rank in (ρ ◦ P̄k|k−1 ) [14]. However, low rank ρ do
not necessarily decrease the rank of (ρ ◦ P̄k|k−1 ). A closely related problem to
finding valid (positive semi-definite or definite) ρ is the creation of covariance
functions and kernels in Gaussian processes [68]. Here, a methodology to create
more complicated kernels from simpler ones could be used to create ρ.
Unfortunately, the Hadamard product cannot be formulated as an operation
on the ensembles in general. Still, the computational requirements can be limited by only working with the non-zero elements of (ρ ◦ P̄k|k−1 ). Furthermore,
it is common to avoid the computation of P̄k|k−1 using
M̄k = ρM ◦ M̄k ,

(30)

instead of (29a) and to skip the tapering in Sk altogether [43]. After all, for
low-dimensional yk (small m) M̄k has the strongest influence on the gain K̄k .
Also the matrix ρM is constructed from prior knowledge about the correlation.
In the geoscientific context, where the state components and measurements are
associated with geographic locations, this is easy. In general, however, it might
not be possible to devise an appropriate ρM . Other variants [14, 26, 65] with
tapering for S̄k exist and have in common that they are only identical to (29)
for H = I.
Some relations between local updates and covariance tapering are discussed
in [65]. For the structure in (27) we can suggest a rank-1 taper ρ that establishes
a correspondence between the two concepts. Let r1 and r2 be vectors of the same
dimensions as x1 and x2 , respectively, that contain all ones. Let r3 be a zero
vector of the same dimension as x3 and rT = [r1T , r2T , r3T ]. Then ρ = rrT removes
14

all entries from P̄k|k−1 that would disappear in (28) anyhow. Furthermore, the
Hadamard product for the rank-1 ρ can be written as an operation on the
ek|k−1 using
ensemble X
(rrT ) ◦ P̄k|k−1 = diag(r)P̄k|k−1 diag(r)


T
ek|k−1 diag(r)X
ek|k−1 .
= N 1−1 diag(r)X

(31)

The multiplication with diag(r) merely removes the rows corresponding to x3 ,
which establishes an equivalence between local updates and covariance tapering.
By picking a smoothly decaying r we can furthermore avoid the discontinuities
associated with local updates.

5.5

The EnKF gain and least squares

A parallel to least squares problems can be disclosed by closer inspection of
the equation (18) that is used to compute the EnKF gain K̄k . Perhaps more
apparent in the transpose of (18), in
T
eT
Yek|k−1 Yek|k−1
K̄kT = Yek|k−1 X
k|k−1 ,

(32a)

appear the normal equations of the least squares problems
T
eT
Yek|k−1
K̄kT = X
k|k−1

(32b)

ek|k−1 . Hence, the EnKF iteration
that are to be solved for each row of K̄k and X
can be carried out without explicitly computing any sample covariance matrices
if instead efficient solutions to the problem (32b) are employed. Furthermore,
the problem (32b) could be modified using regularization [69] to enforce sparsity
in K̄k . This would be an alternative approach to the localization methods
discussed earlier. Related ideas to improve the Kalman gain using bootstrap
methods [69] for computing M̄k and S̄k in (17) are discussed in [70, 71].

6

Relations to other algorithms

The EnKF for nonlinear systems (2) differs from other sampling based nonlinear
filters such as sigma point KF [5] or particle filters (PF) [7]. One reason for this
is that the EnKF approximates the KF algorithm (with the side effect that it
can be applied to (2)) rather than trying to solve the nonlinear filtering problem
directly.
The biggest difference between the EnKF and sigma point filters [5] such
as the unscented KF [4, 53] or divided difference KF [59] is the measurement
update. Whereas the EnKF updates its ensembles, the latter carry out the
KF measurement update (5) using approximately computed mean values and
covariance matrices. That is, the samples or sigma points are condensed into a
filtering estimate x̂k|k and its covariance Pk|k , which entails a loss of information
15

and can be seen as an inherent Gaussian assumption on the filtering density
p(xk |y1:k ). In contrast, the EnKF can preserve more information and deviations
from Gaussianity in the ensemble. Similarities appear in the gain computations
of the EnKF and sigma point KF. In both, the Kalman gain appears as a
function of the sampling covariance matrices, although with the deterministic
sigma points and weights in the latter. With their origin in the KF, both sigma
point filters and the EnKF can be expected to share difficulties with multimodal
posterior distributions.
Similar to the EnKF, the PF propagates N state realizations that are called
particles. For the bootstrap particle filter [6], the prediction step corresponds
to the EnKF time update (11). Apart from that, however, the differences dominate. First, the PF is designed as an approximate solution of the Bayesian
filtering equations [15] using sequential importance sampling [7]. For N → ∞,
the PF solution recovers the true filtering density. Second, the samples in basic
PF variants are generated from a proposal distribution only once every time instance and then left untouched. The measurement update amounts to updating
the particle weights, which leads to a degeneracy problem for large n. In the
EnKF, in contrast, the ensemble members are influenced by the time and the
measurement update. Third, the PF relies on a crucial resampling step that is
not present in the EnKF.
In summary, the EnKF appears as a distinct algorithm besides sigma point
KF and PF. Its properties and potential for nonlinear problems remain to be
fully investigated. Existing results that the EnKF does not converge to the
Bayesian filtering recursion [31] remain to be interpreted in a constructive manner.

7

Instructive Simulation Examples

Four examples are discussed in greater detail, among them one popular benchmark problem of the SP and DA literature each.

7.1

A scalar linear Gaussian model

The first example illustrates the tendency of the EnKF to underestimate the
state uncertainty. A related example is studied in [38]. We compare the EnKF
variance P̄k|k to the Pk|k of the KF via Monte Carlo simulations on the simple
scalar state-space model
xk+1 = xk + vk ,

(33a)

yk = xk + ek .

(33b)

16

Density of the EnKF variances

80
EnKF P

60

k|k

pdf

true P k|k
EnKF P

40

EnKF P

k|k
k|k

mean
median

20
0
0

0.005

0.01

0.015
0.02
0.025
P and EnKF variances

0.03

0.035

0.04

k|k

Figure 1: Distribution of EnKF variances P̄k|k with k = 10 and N = 5 ensemble
members for 10000 runs on the same trajectory. Also shown is the mean and
median of all outcomes and the desired KF variance Pk|k .
The initial state x0 , the process noise vk , and the measurement noise ek are
specified by the probability density functions
p(x0 ) = N (x0 ; 0, 0.1),

(33c)

p(vk ) = N (vk ; 0, 0.1),

(33d)

p(ek ) = N (ek ; 0, 0.01).

(33e)

A trajectory of (33) is simulated and a KF is used to compute the optimal
variances Pk|k . Because the model is time-invariant, the Pk|k quickly converge
to a constant value. For k > 3 Pk|k = 0.0092 is obtained.
Next, 10000 Monte Carlo experiments with a sampling based EnKF with
N = 5 are performed. The distribution of obtained P̄k|k for k = 10 is illustrated
in Fig. 1. The vertical lines indicate the Pk|k of the KF and the median and
mean of the P̄k|k outcomes. The average P̄k|k over the Monte Carlo realizations
is close to the desired Pk|k . However, there is a large spread among the P̄k|k and
the distribution is skewed toward zero with its median below Pk|k . Although
N > n, there is a tendency to underestimate Pk|k .
In order to clarify the reason for this behavior and whether it has to do with
the coupling between the EnKF K̄k and the ensemble members, we repeat the
experiment with an EnKF that uses the gain of the stationary KF for all k. The
resulting outcomes are illustrated in Fig. 2. Now, the average P̄k|k is correct.
However, the median shows that there is still more probability mass below Pk|k .
The tendency to underestimate Pk|k and the remaining spread must be due to
random sampling errors. For larger N the effect vanishes, and the median and
mean of P̄k|k appear similar for N ≥ 10.

7.2

The particle filter benchmark

In the second example we show that the EnKF does not converge to the Bayesian
filtering solution in nonlinear systems as N → ∞ [31]. A well-known benchmark

17

Density of the EnKF variances

80
EnKF P

60

k|k

pdf

true P k|k
EnKF P

40

EnKF P

k|k
k|k

mean
median

20
0
0

0.005

0.01

0.015
0.02
0.025
P and EnKF variances

0.03

0.035

0.04

k|k

Figure 2: Distribution of EnKF variances P̄k|k , but computed with the correct
Kalman gain. Otherwise, similar to Fig. 1.
problem from the PF literature [6] is used. The model is specified by
xk+1 =
yk =

xk
xk
+ 25
+ 8 cos(1.2(k + 1)) + vk ,
2
1 + x2k

(34a)

x2k
+ ek ,
20

(34b)

with independent vk ∼ N (0, 10), ek ∼ N (0, 1), and x0 ∼ N (0, 1). Because
the model is scalar, the Bayesian filtering densities p(xk | y1:k ) can be computed
numerically using point mass filters (PMF) [72]. A sampling based EnKF with
N = 500 is tested and kernel density estimates are used to obtain approximations of p(xk | y1:k ) from the ensembles. For comparison, we include a closely
related sigma point KF variant that uses Monte Carlo integration with N = 500
samples [5]. The only difference to the EnKF is that this Monte Carlo KF
(MCKF) carries out the KF measurement update (5) to propagate a mean and
a variance. We illustrate the results as Gaussian densities.
Fig. 3 shows the prediction results for k = 150. The PMF reference solution
is bimodal with one mode close to the true state. The reason for this lies in
the squared xk in (34b). The EnKF prediction resembles the PMF well except
for the random variations in the kernel density estimate. The MCKF cannot
represent the multimodality but the Gaussian bell covers the relevant regions.
The filtering results for k = 150 are shown in Fig. 4. The PMF reference
solution has much narrower peaks after including yk . The EnKF provides a
skewed density that does not resemble p(xk | y1:k ) even though the EnKF prediction approximated p(xk | y1:k−1 ) well. This is the main take-away result and
confirms [31]. Again, the MCKF exhibits a large variance.
Further filtering results for the PMF and EnKF are shown in Fig. 5. It
can be seen that the EnKF solutions sometimes resemble the PMF very well,
but not always. Similar statements can be made for the prediction results.
Dots in Fig. 5 illustrate the mean values as state estimates. Especially for the
PMF, it can be seen that the mean (though optimal in a minimum variance
sense [3]) is debatable for multimodal densities. Often, all estimates are quite

18

k=150

0.08

PMF
EnKF
MCKF
True

p(x)

0.06
0.04
0.02
0
-30

-20

-10

0

10

20

30

x

Figure 3: Prediction densities p(xk | y1:k−1 ) by the PMF, EnKF, and MCKF for
k = 150. The true state is illustrated as green dot. The PMF serves as reference
solution.

k=150

0.4

PMF
EnKF
MCKF
True

p(x)

0.3
0.2
0.1
0
-30

-20

-10

0

10

20

30

x

Figure 4: Filtering densities p(xk | y1:k ) by PMF, EnKF, and MCKF for k = 150.
Otherwise similar to Fig. 3.

19

PMF
EnKF
True state
PMF mean
EnKF mean
MCKF mean

Filtering densities at time k

125

124

123

122

121

120

-30

-20

-10

0

10

20

30

x

Figure 5: Consecutive filtering densities p(xk | y1:k ) by PMF, EnKF, and MCKF
for k = 120, . . . , 125. Also illustrated are the mean values of the respective
densities and the true state.
close. Fig. 6 provides the estimation error densities obtained from 100 Monte
Carlo experiments with 151 time steps each. The PMF mean estimates exhibit
a larger peak around 0. The estimation errors for the EnKF and MCKF appear
similar. This is surprising because the latter employs a Gaussian approximation
at each time step. Both error densities have heavier tails than the PMF density.
All estimation errors appear unbiased.

7.3

Batch smoothing using the EnKF

We here show how to use the EnKF as smoothing algorithm by estimating
batches of states. This allows us to compare its performance for N < n in
problems of arbitrary dimension.
First, we formulate an “augmented state” that comprises an entire trajectory
of L + 1 steps,

T
ξ = xT0 xT1 . . . xTL ,
(35)
with dimension n = (L + 1)nx . Second, we note that the measurements yk ,
k = 1, . . . , L, have uncorrelated measurement noise and known relations to the
components of ξ. For linear systems (1), the predicted mean and covariance of ξ
can be easily derived, and smoothed estimates of all xk , k = 0, . . . , L, can be
obtained by sequentially processing all yk in KF measurement updates for ξ.

20

Error pdf

0.3
PMF
EnKF
MCKF

0.2

0.1

0
-15

-10

-5

0

5

10

15

Estimation error

Figure 6: Density of the estimation errors obtained from 100 Monte Carlo runs
with 151 time steps each.
Also other smoothing variants and the Rauch-Tung-Striebel (RTS) algorithm
can be derived from state augmentation approaches [3]. Due to its sequential
nature, however, the RTS smoother does not provide joint covariance matrices of
xk and xk+i for i 6= 0. Except for this and the higher computational complexity
of working with ξ, the batch and RTS smoothers are equivalent for (1).
An EnKF approach to batch smoothing mimics the above. A prediction
ensemble for ξ is obtained by simulating N trajectories for random process
noise and initial state realizations. This can also be carried out for nonlinear
models (2). Then, sequential EnKF measurement updates are performed for
all yk .
For our experiments we use a tracking problem with a constant velocity
model [64] and position measurements. The low-dimensional state is given by

x= x

y

ẋ ẏ

T

(36a)

and comprises the Cartesian position [m] and velocity [m/s] of an object. The
parameters of (1) are given by


 T2 


I2 TI2
I
F =
,
G= 2 2 ,
H = I2 0 ,
(36b)
0
I2
TI2
with T = 1 s. The initial state x0 is Gaussian distributed with

x̂0 = 0

0

15

−10

T

,

P0 = diag(502 , 502 , 202 , 202 ),

and the process and measurement noise covariances are


2000 1000
Q = diag(10, 50), R =
.
1000 1980

(36c)

(36d)

With nx = 4 and L = 49 we obtain n = 200 as dimension of ξ. The RTS
solution is compared to EnKF of ensemble size N = {10, 20, 50}. Monte Carlo
errors are reduced using (21) in the gain computations.
21

Position in km

0.2
0
-0.2
-0.4
-1.5

-1

-0.5

0

0.5

Position in km

Figure 7: Illustration of a representative trajectory (black), the RTS smoothing
solution (cyan), and an initial ensemble (N = 50, orange). Red circles depict
the measurements. Most ensemble trajectories go beyond the plot area.

Position in km

0.2
0
-0.2
-0.4
-1.5

-1

-0.5

0

0.5

Position in km

Figure 8: The ensemble of Fig. 7 after a measurement update with yL only.
Some ensemble trajectories leave and re-enter the plot area.
A realization of a true trajectory and its measurements is provided in Fig. 7
together with the RTS estimate and an ensemble of N = 50 trajectories. The
latter are the initial ensemble of an EnKF. The ensemble is well gathered around
the initial position but fans out wildly. Fig. 8 shows the ensemble after an update
with yL only. The measurement at the end of the trajectory provides an anchor
point and quickly reduces the spread of the ensemble. Fig. 9 shows the result
after processing all measurements in sequential order from first to last. The true
trajectory and the RTS estimate are mostly covered well by the ensemble. The
EnKF with N = 50 appears consistent in this respect. Position errors for the
RTS and the EnKF are provided in Fig. 10. The EnKF performs slightly worse
than the RTS but still gives good results for N = 50, without extra inflation or
localization.
The next experiment explores the EnKF for N = 10. Fig. 11 shows the ensemble after processing all measurements. The ensemble is compactly gathered
but does not cover the true trajectory well. The EnKF is overconfident.
A last experiment explores how well an EnKF with N = 20 captures the
uncertainty of the state estimate. Furthermore, we discuss effects of the order

22

Position in km

0.2
0
-0.2
-0.4
-1.5

-1

-0.5

0

0.5

Position in km

Figure 9: The ensemble of Fig. 7 after updating with all measurements in the
order y1 , . . . yL . The RTS solution is covered well.

Position error

60

RTS
EnKF

40

20

0

5

10

15

20

25

30

35

40

45

Time step k

Figure 10: Position errors of the RTS (cyan) and the EnKF (N = 50, orange)
after updating with all measurements in the order y1 , . . . yL .

Position in km

0.2
0
-0.2
-0.4
-1.5

-1

-0.5

0

0.5

Position in km

Figure 11: An ensemble with N = 10 after updating with all measurements in
the order y1 , . . . yL . The smaller ensemble is more condensed and does not cover
the RTS solution well.

23

500

0

400

time step k

10
300

20

200
100

30

0

40

-100
-200

0

10

20

30

40

time step k

Figure 12: Exact position covariance matrix cov(xi , xj ) after including all measurements.
500

0

400

time step k

10
300

20

200
100

30

0

40

-100
-200

0

10

20

30

40

time step k

Figure 13: EnKF (20 members) position covariance matrix cov(xi , xj ) after
including all measurements in the order y1 , . . . yL .
in which the measurements are processed. Specifically, we compare the ensemble
covariance of the positions xk to the exact cov(xk , xi ), i, k = 0, . . . , L, obtained
by KF updates for the augmented state ξ.
The exact covariance after processing all measurements is illustrated in
Fig. 12. Row k in the matrix defines the covariance function between xk and the
remaining x positions. The banded structure indicates that subsequent positions
are more related than, say, x0 and xL .
Fig. 13 shows the corresponding EnKF covariance after processing the measurements from y1 to yL . The off-diagonal elements do not decay uniformly
as in Fig. 12 and spurious positive and negative correlations appear. Furthermore, the correct temporal order of measurements entails an unwanted structure. Later xk are rated more uncertain according to the lighter areas in the
lower right corner of Fig. 13.
A covariance after processing the measurements in random order is shown
in Fig. 14. The spurious correlations persist but the diagonal elements appear
more homogeneous.
From the above experiments we conclude that the EnKF can provide good
estimates for ensembles with N < n. However, there is a minimum N required

24

500

0

400

time step k

10
300

20

200
100

30

0

40

-100
-200

0

10

20

30

40

time step k

Figure 14: EnKF (20 members) position covariance matrix cov(xi , xj ) after
including all measurements in random order.
to obtain consistent results without further measures such as localization or
inflation. We have shown adverse effects such as ensembles with too little spread
and spurious correlations.

7.4

The 40-dimensional Lorenz model

Our final example is a benchmark problem from the EnKF literature. We
investigate the 40-dimensional Lorenz-96 model3 from [50] that is used in, e.g.,
[36, 38, 42, 47, 49, 60, 66]. The state x mimics an atmospheric quantity at
equally spaced locations along a circle. Its evolution is specified by the nonlinear
differential equation


ẋ(j) = x(j + 1) − x(j − 2) x(j − 1) − x(j) + F(j),
(37)
where j = 1, . . . , 40 indexes the components of x, with the convention that
x(0) = x(40) etc. Instead of the commonly used forcing term F(j) = 8, we
assume time-dependent Fk (j) ∼ N (8, 1) that are constant for time intervals
T = 0.05 only and act as process noise. A Runge-Kutta method (RK4) is used
to discretize (37) to obtain the nonlinear state difference equation (2a) with
xk = xk and vk = Fk . The step size T corresponds to about six hours if x were
an atmospheric quantity on a latitude circle of the earth [50]. Although the
model (37) is said to be chaotic, the effects are only mild for short integration
times T. In our experiments all n = 40 states are measured with additive
Gaussian noise ek ∼ N (0, I). The initial state is Gaussian with x0 ∼ N (0, P0 ),
where P0 is drawn from a Wishart distribution with seed matrix In and n
degrees of freedom. Fig. 15 illustrates how the state evolves over several time
steps. There is a tendency for peaks to move “westwards” as k increases.
We note that there are also alternative approaches for estimating x, for
example, by first linearizing and then discretizing (37). However, we adopt the
RK4 discretization of the EnKF literature that yields a state transition that is
3 Also

known as the Lorenz-96, L95, L96, or L40 model.

25

State at time step k

60
55
50
45
40
35
30
5

10

15

20

25

30

35

40

State component

Figure 15: State evolution for the Lorenz model. Each horizontal line “carries”
a 40-dimensional state vector.

State component

0.4
0.3

10

0.2

20

0.1
0

30

-0.1

40
10

20

30

40

State component

Figure 16: Prediction covariance P̄k|k−1 for k = 30 obtained from an EnKF
with N = 1000. The banded structure justifies the use of localization.
easy to evaluate but difficult to linearize. Because of this, the EKF [3] cannot
be applied easily and we obtain a challenging benchmark problem.
We use sampling based EnKF to estimate long state sequences of L = 104
time steps. Following [38, 42], the performance is assessed by the error
r
1
(x̂k|k − xk )T (x̂k|k − xk ),
(38)
εk =
n
where x̂k|k is the ensemble mean. We use the average εk for k = 100, . . . , L,
denoted by ε̄, as quantitative performance measure for different EnKF. Useful
EnKF must yield ε̄ < 1, which is the error when simply taking x̂k|k = yk .
First, we compute a reference solution using an EnKF with N = 1000.
Without any localization or inflation ε̄ = 0.29 is achieved. Fig. 16 shows the
sample covariance P̄k|k−1 of a prediction ensemble Xk|k−1 , our best guess of
the true covariance. The banded structure reveals that the problem is suitable
for localization. Hence, we construct a matrix ρ for covariance tapering from a
compactly supported correlation function [67] that is also used in [14, 26, 38, 43]
26

State component

1
0.8

10

0.6

20
0.4

30

0.2

40

0

10

20

30

40

State component

Figure 17: The employed tapering matrix ρ.
and appears to be the standard choice. The chosen ρ is a Toeplitz matrix because
the components of xk are at equidistant locations, and shown in Fig. 17.
Next, EnKF with different ensemble sizes N , ensemble inflation factors c,
with or without tapering, are compared. The obtained errors ε̄ are summarized
in Table 1. For N = n = 40, we obtain a worse ε̄ than for N = 1000. While
Table 1: Averaged errors ε̄ for
N
c
ρ
1000 1
no
40
1
no
40
1.05 no
40
1
yes
40
1.02 yes
20
> 1 no
20
1.01 yes
10
1.05 yes

different EnKF.
ε̄
0.29
0.44
0.33
0.29
0.28
>1
0.3
0.34

inflation without tapering does reduce the error slightly, the covariance tapering
even yields a better result that the EnKF with N = 1000. Further improvements
are obtained by combining inflation and tapering. Fig. 18 shows the estimation
error xk − x̂k|k for k = 104 , N = 40, c = 1.02, and tapering with ρ. In the
ek|k are illustrated. The estimation error
background, the ensemble deviations X
is mostly contained in the intervals spanned by the ensemble, hence the EnKF
is consistent. Tests on EnKF with N = 20 reveal convergence problems, even
with inflation the initial estimation error persists. With the help of tapering,
however, a competitive error can be achieved. Even further reduction to N = 10
is possible with tapering and inflation. The required inflation factor c must be
increased to counteract the lack of ensemble spread. Similar to Fig. 18, Fig. 19
illustrates the estimation error and deviation ensemble for k = 104 , N = 10,
c = 1.05, and tapering with ρ. Although the obtained error is larger than for
N = 40, the ensemble deviations represent the estimation uncertainty well.
A number of lessons have been learned from related experiments. As alter27

Deviations/estimation error

1
0.5
0
-0.5
-1
5

10

15

20

25

30

35

40

State component

Deviations/estimation error

Figure 18: The estimation error xk − x̄k|k for k = 104 with the deviation ensemek|k in the background for an EnKF with N = 40, covariance localization,
ble X
and inflation factor c = 1.02.

0.5

0

-0.5

5

10

15

20

25

30

35

40

State component

Figure 19: The estimation error xk − x̄k|k for k = 104 with the deviation ensemek|k in the background for an EnKF with N = 10, covariance localization,
ble X
and inflation factor c = 1.05.

28

native to the ρ in Fig. 17, a simpler taper that contains only ones and zeros
to enforce the banded structure was used. Although this ρ was indefinite, a
reduction in ε̄ was achieved without any numerical issues. Hence, the specific
structure of ρ appears secondary. The smooth ρ of Fig. 17 remains preferable in
terms of ε̄, though. Sequential processing of the measurements did not degrade
the performance. Experiments without process noise give the lower errors ε̄
from, e.g., [38, 42].

8

Concluding Remarks

With this paper, we have given a comprehensive and easy to understand introduction to the EnKF for signal processing researchers. The origin of the
EnKF in the KF and its simple implementation have been demonstrated. The
unique literature review provides quick access to the most relevant papers in the
plethora of geoscientific EnKF publications. Furthermore, we have discussed the
challenges related to small ensembles for high-dimensional states, N < n, and
the available solutions such as localization or inflation. Finally, we have tested
the EnKF on signal processing and EnKF benchmark problems.
With its scalability and simple implementation, even for nonlinear and nonGaussian problems, the EnKF stands out as viable candidate for many state
estimation problems. Furthermore, localization ideas and advanced concepts
for estimating covariance matrices and the EnKF gain from the limited information in the ensembles provide new research directions for the EnKF and
high-dimensional filters in general, hopefully with an increased participation
from the signal processing community.

Acknowledgments
This work was supported by the project Scalable Kalman Filters granted by the
Swedish Research Council (VR).

References
[1] E. Kalnay, Atmospheric Modeling, Data Assimilation and Predictability.
New York: Cambridge University Press, Dec. 2002.
[2] R. E. Kalman, “A new approach to linear filtering and prediction problems,” Journal of basic Engineering, vol. 82, no. 1, pp. 35–45, Mar. 1960.
[3] B. D. Anderson and J. B. Moore, Optimal Filtering.
1979.

Prentice Hall, Jun.

[4] S. Julier, J. Uhlmann, and H. Durrant-Whyte, “A new approach for filtering
nonlinear systems,” in Proceedings of the American Control Conference
1995, vol. 3, 1995, pp. 1628–1632.
29

[5] M. Roth, G. Hendeby, and F. Gustafsson, “Nonlinear Kalman filters explained: A tutorial on moment computations and sigma point methods,”
Journal of Advances in Information Fusion, vol. 11, no. 1, pp. 47–70, Jun.
2016.
[6] N. J. Gordon, D. J. Salmond, and A. F. Smith, “Novel approach to
nonlinear/non-Gaussian Bayesian state estimation,” Radar and Signal Processing, IEE Proceedings F, vol. 140, no. 2, pp. 107–113, Apr. 1993.
[7] F. Gustafsson, “Particle filter theory and practice with positioning applications,” IEEE Aerospace and Electronic Systems Magazine, vol. 25, no. 7,
pp. 53–82, 2010.
[8] G. Evensen, “Sequential data assimilation with a nonlinear quasigeostrophic model using Monte Carlo methods to forecast error statistics,”
Journal of Geophysical Research: Oceans, vol. 99, no. C5, pp. 10 143–10 162,
1994.
[9] G. Burgers, P. Jan van Leeuwen, and G. Evensen, “Analysis scheme in
the ensemble Kalman filter,” Monthly Weather Review, vol. 126, no. 6, pp.
1719–1724, Jun. 1998.
[10] H. Durrant-Whyte and T. Bailey, “Simultaneous localization and mapping:
Part I,” IEEE Robotics Automation Magazine, vol. 13, no. 2, pp. 99–110,
Jun. 2006.
[11] M. Baum and U. D. Hanebeck, “Extended object tracking with random hypersurface models,” IEEE Transactions on Aerospace and Electronic Systems, vol. 50, no. 1, pp. 149–159, Jan. 2014.
[12] N. Wahlström and E. Özkan, “Extended target tracking using Gaussian
processes,” IEEE Transactions on Signal Processing, vol. 63, no. 16, pp.
4165–4178, Aug. 2015.
[13] P. L. Houtekamer and H. L. Mitchell, “Data assimilation using an ensemble
Kalman filter technique,” Monthly Weather Review, vol. 126, no. 3, pp.
796–811, Mar. 1998.
[14] ——, “A sequential ensemble Kalman filter for atmospheric data assimilation,” Monthly Weather Review, vol. 129, no. 1, pp. 123–137, Jan. 2001.
[15] A. H. Jazwinski, Stochastic Processes and Filtering Theory.
Press, Mar. 1970.

Academic

[16] G. Evensen, “The ensemble Kalman filter: theoretical formulation and
practical implementation,” Ocean Dynamics, vol. 53, no. 4, pp. 343–367,
Nov. 2003.
[17] ——, Data Assimilation: The Ensemble Kalman Filter, 2nd ed. Dordrecht;
New York: Springer, Aug. 2009.
30

[18] T. M. Hamill, “Ensemble-based atmospheric data assimilation,” in Predictability of Weather and Climate. Cambridge University Press, 2006.
[19] P. L. Houtekamer and H. L. Mitchell, “Ensemble Kalman filtering,” Quarterly Journal of the Royal Meteorological Society, vol. 131, no. 613, pp.
3269–3289, Oct. 2005.
[20] J. S. Whitaker, T. M. Hamill, X. Wei, Y. Song, and Z. Toth, “Ensemble
data assimilation with the NCEP global forecast system,” Monthly Weather
Review, vol. 136, no. 2, pp. 463–482, Feb. 2008.
[21] G. P. Compo, J. S. Whitaker, P. D. Sardeshmukh, N. Matsui, R. J. Allan, X. Yin, B. E. Gleason, R. S. Vose, G. Rutledge, P. Bessemoulin,
S. Brönnimann, M. Brunet, R. I. Crouthamel, A. N. Grant, P. Y. Groisman, P. D. Jones, M. C. Kruk, A. C. Kruger, G. J. Marshall, M. Maugeri,
H. Y. Mok, . Nordli, T. F. Ross, R. M. Trigo, X. L. Wang, S. D. Woodruff,
and S. J. Worley, “The twentieth century reanalysis project,” Quarterly
Journal of the Royal Meteorological Society, vol. 137, no. 654, pp. 1–28,
Jan. 2011.
[22] S. Lakshmivarahan and D. Stensrud, “Ensemble Kalman filter,” IEEE Control Systems, vol. 29, no. 3, pp. 34–46, Jun. 2009.
[23] J. Anderson, “Ensemble Kalman filters for large geophysical applications,”
IEEE Control Systems, vol. 29, no. 3, pp. 66–82, Jun. 2009.
[24] G. Evensen, “The ensemble Kalman filter for combined state and parameter
estimation,” IEEE Control Systems, vol. 29, no. 3, pp. 83–104, Jun. 2009.
[25] J. Mandel, J. Beezley, J. Coen, and M. Kim, “Data assimilation for wildland
fires,” IEEE Control Systems, vol. 29, no. 3, pp. 47–65, Jun. 2009.
[26] R. Furrer and T. Bengtsson, “Estimation of high-dimensional prior and
posterior covariance matrices in Kalman filter variants,” Journal of Multivariate Analysis, vol. 98, no. 2, pp. 227–255, Feb. 2007.
[27] M. Butala, J. Yun, Y. Chen, R. Frazin, and F. Kamalabadi, “Asymptotic
convergence of the ensemble Kalman filter,” in 15th IEEE International
Conference on Image Processing, Oct. 2008, pp. 825–828.
[28] J. Mandel, L. Cobb, and J. D. Beezley, “On the convergence of the ensemble
Kalman filter,” Applications of Mathematics, vol. 56, no. 6, pp. 533–541,
Dec. 2011.
[29] M. Frei, “Ensemble Kalman filtering and generalizations,” Dissertation,
ETH, Zürich, 2013, nr. 21266.
[30] M. Katzfuss, J. R. Stroud, and C. K. Wikle, “Understanding the ensemble
Kalman filter,” The American Statistician, pp. 350–357, Feb. 2016.

31

[31] F. Le Gland, V. Monbet, and V. Tran, “Large sample asymptotics for the
ensemble Kalman filter,” in The Oxford Handbook of Nonlinear Filtering,
D. Crisan and B. Rozovskii, Eds. Oxford University Press, 2011, pp.
598–634.
[32] M. Butala, R. Frazin, Y. Chen, and F. Kamalabadi, “Tomographic Imaging
of dynamic objects with the ensemble Kalman filter,” IEEE Transactions
on Image Processing, vol. 18, no. 7, pp. 1573–1587, Jul. 2009.
[33] J. Dunik, O. Straka, M. Simandl, and E. Blasch, “Random-point-based
filters: analysis and comparison in target tracking,” IEEE Transactions on
Aerospace and Electronic Systems, vol. 51, no. 2, pp. 1403–1421, Apr. 2015.
[34] S. Gillijns, O. Mendoza, J. Chandrasekar, B. De Moor, D. Bernstein, and
A. Ridley, “What is the ensemble Kalman filter and how well does it work?”
in American Control Conference, 2006, Jun. 2006, pp. 4448–4453.
[35] M. Roth, C. Fritsche, G. Hendeby, and F. Gustafsson, “The ensemble
Kalman filter and Its relations to other nonlinear filters,” in European Signal Processing Conference 2015 (EUSIPCO 2015), Nice, France, Aug. 2015.
[36] J. L. Anderson, “An ensemble adjustment Kalman filter for data assimilation,” Monthly Weather Review, vol. 129, no. 12, pp. 2884–2903, Dec.
2001.
[37] C. H. Bishop, B. J. Etherton, and S. J. Majumdar, “Adaptive sampling
with the ensemble transform Kalman filter. Part I: Theoretical aspects,”
Monthly Weather Review, vol. 129, no. 3, pp. 420–436, Mar. 2001.
[38] J. S. Whitaker and T. M. Hamill, “Ensemble data assimilation without
perturbed observations,” Monthly Weather Review, vol. 130, no. 7, pp.
1913–1924, Jul. 2002.
[39] M. K. Tippett, J. L. Anderson, C. H. Bishop, T. M. Hamill, and J. S.
Whitaker, “Ensemble square root filters,” Monthly Weather Review, vol.
131, no. 7, pp. 1485–1490, Jul. 2003.
[40] J. L. Anderson and S. L. Anderson, “A Monte Carlo Implementation of
the nonlinear filtering problem to produce ensemble assimilations and forecasts,” Monthly Weather Review, vol. 127, no. 12, pp. 2741–2758, Dec.
1999.
[41] P. J. van Leeuwen, “Comment on “data assimilation using an ensemble
Kalman filter technique”,” Monthly Weather Review, vol. 127, no. 6, pp.
1374–1377, Jun. 1999.
[42] E. Ott, B. R. Hunt, I. Szunyogh, A. V. Zimin, E. J. Kostelich, M. Corazza,
E. Kalnay, D. J. Patil, and J. A. Yorke, “A local ensemble Kalman filter
for atmospheric data assimilation,” Tellus A, vol. 56, no. 5, Oct. 2004.

32

[43] T. M. Hamill, J. S. Whitaker, and C. Snyder, “Distance-dependent filtering
of background error covariance estimates in an ensemble Kalman filter,”
Monthly Weather Review, vol. 129, no. 11, pp. 2776–2790, Nov. 2001.
[44] P. J. van Leeuwen, “A variance-minimizing filter for large-scale applications,” Monthly Weather Review, vol. 131, no. 9, pp. 2071–2084, Sep. 2003.
[45] C. Snyder, T. Bengtsson, P. Bickel, and J. Anderson, “Obstacles to highdimensional particle filtering,” Monthly Weather Review, vol. 136, no. 12,
pp. 4629–4640, Dec. 2008.
[46] P. J. van Leeuwen, “Particle filtering in geophysical systems,” Monthly
Weather Review, vol. 137, no. 12, pp. 4089–4114, Dec. 2009.
[47] ——, “Nonlinear data assimilation in geosciences: an extremely efficient
particle filter,” Quarterly Journal of the Royal Meteorological Society, vol.
136, no. 653, pp. 1991–1999, Oct. 2010.
[48] M. Frei and H. R. Künsch, “Bridging the ensemble Kalman and particle
filters,” Biometrika, vol. 100, no. 4, pp. 781–800, Dec. 2013.
[49] J. Poterjoy, “A localized particle filter for high-dimensional nonlinear systems,” Monthly Weather Review, vol. 144, no. 1, pp. 59–76, Oct. 2015.
[50] E. N. Lorenz, “Predictability — a problem partly solved,” in Predictability
of Weather and Climate, T. Palmer and R. Hagedorn, Eds. Cambridge
University Press, 2006, pp. 40–58.
[51] D. T. Pham, “Stochastic methods for sequential data assimilation in
strongly nonlinear systems,” Monthly Weather Review, vol. 129, no. 5, pp.
1194–1207, May 2001.
[52] X. Luo and I. Moroz, “Ensemble Kalman filter with the unscented transform,” Physica D: Nonlinear Phenomena, vol. 238, no. 5, pp. 549–562, Mar.
2009.
[53] S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,” Proceedings of the IEEE, vol. 92, no. 3, pp. 401– 422, Mar. 2004.
[54] P. Sakov, “Comment on “ensemble Kalman filter with the unscented transform”,” Physica D: Nonlinear Phenomena, vol. 238, no. 22, pp. 2227–2228,
Nov. 2009.
[55] A. S. Stordal, H. A. Karlsen, G. Nævdal, H. J. Skaug, and B. Vallès,
“Bridging the ensemble Kalman filter and particle filters: the adaptive
Gaussian mixture filter,” Computational Geosciences, vol. 15, no. 2, pp.
293–305, Mar. 2011.
[56] I. Hoteit, X. Luo, and D.-T. Pham, “Particle Kalman filtering: A nonlinear Bayesian framework for ensemble Kalman filters,” Monthly Weather
Review, vol. 140, no. 2, pp. 528–542, Aug. 2011.
33

[57] M. Frei and H. R. Künsch, “Mixture ensemble Kalman filters,” Computational Statistics & Data Analysis, vol. 58, pp. 127–138, Feb. 2013.
[58] L. N. Trefethen and David Bau, III, Numerical Linear Algebra. Philadelphia: SIAM, Jun. 1997.
[59] M. Nørgaard, N. K. Poulsen, and O. Ravn, “New developments in state
estimation for nonlinear systems,” Automatica, vol. 36, no. 11, pp. 1627–
1638, Nov. 2000.
[60] P. Sakov and P. R. Oke, “Implications of the form of the ensemble transformation in the ensemble square root filters,” Monthly Weather Review,
vol. 136, no. 3, pp. 1042–1053, Mar. 2008.
[61] D. M. Livings, S. L. Dance, and N. K. Nichols, “Unbiased ensemble square
root filters,” Physica D: Nonlinear Phenomena, vol. 237, no. 8, pp. 1021–
1028, Jun. 2008.
[62] W. G. Lawson and J. A. Hansen, “Implications of stochastic and deterministic filters as ensemble-based data assimilation methods in varying regimes
of error growth,” Monthly Weather Review, vol. 132, no. 8, pp. 1966–1981,
Aug. 2004.
[63] O. Leeuwenburgh, G. Evensen, and L. Bertino, “The impact of ensemble
filter definition on the assimilation of temperature profiles in the tropical
pacific,” Quarterly Journal of the Royal Meteorological Society, vol. 131,
no. 613, pp. 3291–3300, Oct. 2005.
[64] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with Applications
to Tracking and Navigation. Wiley-Interscience, Jun. 2001.
[65] P. Sakov and L. Bertino, “Relation between two common localisation methods for the EnKF,” Computational Geosciences, vol. 15, no. 2, pp. 225–237,
Jul. 2010.
[66] B. R. Hunt, E. J. Kostelich, and I. Szunyogh, “Efficient data assimilation for
spatiotemporal chaos: A local ensemble transform Kalman filter,” Physica
D: Nonlinear Phenomena, vol. 230, no. 1–2, pp. 112–126, Jun. 2007.
[67] G. Gaspari and S. E. Cohn, “Construction of correlation functions in two
and three dimensions,” Quarterly Journal of the Royal Meteorological Society, vol. 125, no. 554, pp. 723–757, Jan. 1999.
[68] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine
Learning. Cambridge, Mass: The MIT Press, Nov. 2005.
[69] T. Hastie, R. Tibshirani, and J. Friedman, The Elements of Statistical
Learning: Data Mining, Inference, and Prediction, 2nd ed. New York,
NY: Springer, Apr. 2011.

34

[70] Y. Zhang and D. S. Oliver, “Improving the ensemble estimate of the
Kalman gain by bootstrap sampling,” Mathematical Geosciences, vol. 42,
no. 3, pp. 327–345, Feb. 2010.
[71] I. Myrseth, J. Sætrom, and H. Omre, “Resampling the ensemble Kalman
filter,” Computers & Geosciences, vol. 55, pp. 44–53, Jun. 2013.
[72] M. Roth and F. Gustafsson, “Computation and visualization of posterior densities in scalar nonlinear and non-Gaussian Bayesian filtering
and smoothing problems,” in 42nd International Conference on Acoustics,
Speech, and Signal Processing (ICASSP), New Orleans, USA, Mar. 2017.

35

