Well-Posedness And Accuracy Of The Ensemble Kalman Filter In
Discrete And Continuous Time
October 14, 2013
D.T.B. Kelly, K.J.H. Law, A.M. Stuart
The University of Warwick, Email: David.Kelly@Warwick.ac.uk
Abstract
The ensemble Kalman filter (EnKF) is a method for combining a dynamical model with data in a sequential
fashion. Despite its widespread use, there has been little analysis of its theoretical properties. Many of the
algorithmic innovations associated with the filter, which are required to make a useable algorithm in practice, are
derived in an ad hoc fashion. The aim of this paper is to initiate the development of a systematic analysis of the
EnKF, in particular to do so in the small ensemble size limit. The perspective is to view the method as a state
estimator, and not as an algorithm which approximates the true filtering distribution. The perturbed observation
version of the algorithm is studied, without and with variance inflation. Without variance inflation well-posedness
of the filter is established; with variance inflation accuracy of the filter, with resepct to the true signal underlying the
data, is established. The algorithm is considered in discrete time, and also for a continuous time limit arising when
observations are frequent and subject to large noise. The underlying dynamical model, and assumptions about it,
is sufficiently general to include the Lorenz ’63 and ’96 models, together with the incompressible Navier-Stokes
equation on a two-dimensional torus. The analysis is limited to the case of complete observation of the signal with
additive white noise. Numerical results are presented for the Navier-Stokes equation on a two-dimensional torus
for both complete and partial observations of the signal with additive white noise.
1 Introduction
In recent years the ensemble Kalman filter (EnKF) [Eve06] has become a widely used methodology for combin-
ing dynamical models with data. The algorithm is used in oceanography, oil reservoir simulation and weather
prediction [BVLE98, EVL00, Kal03, ORL08], for example. The basic idea of the method is to propagate an
ensemble of particles to describe the distribution of the signal given data, employing empirical second order
statistics to update the distribution in a Kalman-like fashion when new data is acquired. Despite the widespread
use of the method, its behaviour is not well understood. In contrast with the ordinary Kalman filter, which ap-
plies to linear Gaussian problems, it is difficult to find a mathematical justification for EnKF. The most notable
progress in this direction can be found in [LGMT
+
10, MCB11], where it is proved that, for linear dynamics,
the EnKF approximates the usual Kalman filter in the large ensemble limit. This analysis is however far from
being useful for practitioners who typically run the method with small ensemble size on nonlinear problems.
Furthermore there is an accumulation of numerical evidence showing that the EnKF, and related schemes such
as the extended Kalman filter, can “diverge” with the meaning of “diverge” ranging from simply loosing the
true signal through to blow-up [IKJ02, MH08, GM13]. The aim of our work is to try and build mathematical
foundations for the analysis of the EnKF, in particular with regards to well-posedness (lack of blow-up) and
accuracy (tracking the signal over arbitrarily long time-intervals). To make progress on such questions it is
necessary to impose structure on the underlying dynamics and we choose to work with dissipative quadratic
systems with energy-conserving nonlinearity, a class of problems which has wide applicability [MW06] and
which has proved to be useful in the development of filters [MH12].
In section 3 we derive the perturbed observation form of the EnKF and demonstrate how it links to the
randomized maximum likelihood method (RML) which is widely used in oil reservoir simulation [ORL08].
1
arXiv:1310.3167v1 [math.PR] 11 Oct 2013
We also introduce the idea of variance inflation, widely used in many practical implementations of the EnKF
[And07]. Section 4 contains theoretical analyses of the perturbed observation EnKF, without and with variance
inflation. Without variance inflation we are able only to prove bounds which grow exponentially in the discrete
time increment underlying the algorithm (Theorem 4.2); with variance inflation we are able to prove filter
accuracy and show that, in mean square with respect to the noise entering the algorithm, the filter is uniformly
close to the true signal for all large times, provided enough inflation is employed (Theorem 4.4). These results,
and in particular the one concerning variance inflation, are similar to the results developed in [BLL
+
12] for the
3DVAR filter applied to the Navier-Stokes equation and for the 3DVAR filter applied to the Lorenz ’63 model
in [LSS14], as well as the similar analysis developed in [MLPvL13] for the 3DVAR filter applied to globally
Lipschitz nonlinear dynamical systems. In section 5 we describe a continuous time limit in which data arrives
very frequently, but is subject to large noise. If these two effects are balanced appropriately a stochastic (partial)
differential equation limit is found and it is instructive to study this limiting continuous time process. This
idea was introduced in [BLSZ] for the 3DVAR filter and is here employed for the EnKF filter. The primary
motivation for the continuous time limit is to obtain insight into the mechanisms underlying the EnKF; some of
these mechanisms are more transparent in continuous time. In section 6 we analyze the well-posedness of the
continuous time EnKF (Theorem 6.2). Section 7 contains numerical experiments which illustrate and extend
the theory, and section 8 contains some brief concluding remarks.
Throughout the sequel we use the following notation. Let H be a separable Hilbert space with norm |·| and
inner product h·, ·i. For a linear operator C on H, we will write C ≥ 0 (resp. C > 0) when C is self-adjoint
and positive semi-definite (resp. positive definite). Given C > 0, we will denote
|·|
C
def
=
C
−1/2
(·)
.
2 Set-Up
2.1 Filtering Distribution
We assume that the observed dynamics are governed by an evolution equation
du
dt
= F (u) (2.1)
which generates a one-parameter semigroup Ψ
t
: H → H. We also assume that K ⊂ H is another Hilbert
space, which acts as the observation space. We assume that noisy observations are made in K every h time
units and write Ψ = Ψ
h
. We define u
j
= u(jh) for j ∈ N and, assuming that u
0
is uncertain and modelled as
Gaussian distributed, we obtain
u
j+1
= Ψ(u
j
) , with u
0
∼ N (m
0
, C
0
)
for some initial mean m
0
and covariance C
0
. We are given the observations
y
j+1
= Hu
j+1
+ Γ
1/2
ξ
j+1
, with ξ
j
∼ N (0, I) i.i.d. ,
where H ∈ L(H, K) is the observation operator and Γ ∈ L(H, H) with Γ ≥ 0 is the covariance operator of
the observational noise; the i.i.d. noise sequence {ξ
j
} is assumed independent of u
0
. The aim of filtering is to
approximate the distribution of u
j
given Y
j
= {y
`
}
j
`=1
using a sequential update algorithm. That is, given the
distribution u
j
|Y
j
as well as the observation y
j+1
, find the distribution of u
j+1
|Y
j+1
. We refer to the sequence
P(u
j
|Y
j
) as the filtering distribution.
2.2 Assumptions
To write down the EnKF as we do in section 3, and indeed to derive the continuum limit of the EnKF, as we do in
section 5, we need make no further assumptions about the underlying dynamics and observation operator other
than those made above. However, in order to analyze the properties of the EnKF, as we do in sections 4 and 6,
we will need to make structural assumptions and we detail these here. It is worth noting that the assumptions we
2
make on the underlying dynamics are met by several natural models used to test data assimilation algorithms.
In particular, the 2D Navier-Stokes equations on a torus, as well as both Lorenz ’63 and ’96 models, satisfy
Assumptions 2.3 [MW06, MH12, Tem97].
Assumption 2.3. (Dynamics Model) Suppose there is some Banach space V, equipped with norm k·k, that can
be continuously embedded into H. We assume that (2.1) has the form
du
dt
+ Au + B(u, u) = f , (2.2)
where A : H → H is an unbounded linear operator satisfying
hAu, ui ≥ λ kuk
2
, (2.3)
for some λ > 0, B is a symmetric bilinear operator B : V × V → H and f : R
+
→ H. We furthermore assume
that B satisfies the identity
hB(u, u), ui = 0 , (2.4)
for all u ∈ H and also
hB(u, v), vi ≤ c kuk kvk |v| , (2.5)
for all u, v, w ∈ H, where c > 0 depends only on the bilinear form. We assume that the equation (2.2) has
a unique weak solution for all u(0) ∈ H, and generates a one-parameter semigroup Ψ
t
: V → V which may
be extended to act on H. Finally we assume that there exists a global attractor Λ ⊂ V for the dynamics, and
constant R > 0 such that for any initial condition u
0
∈ Λ, we have that sup
t≥0
ku(t)k ≤ R.
Remark 2.4. In the finite dimensional case the final assumption on the existence of a global attractor does not
need to be made as it is a consequence of the preceding assumptions made. To see this note that
1
2
d|u|
2
dt
+ λkuk
2
≤ hf, ui. (2.6)
The continuous embedding of V, together with the Cauchy-Schwarz inequality, implies the existence of a strictly
positive constant such that
1
2
d|u|
2
dt
+ |u|
2
≤
1
2δ
|f|
2
+
δ
2
|u|
2
(2.7)
for all δ > 0. Choosing δ = gives the existence of an absorbing set and hence a global attractor by Theorem
1.1 in Chapter I of [Tem97]. In infinite dimensions the existence of a global attractor in V follows from the
techniques in [Tem97] for the Navier-Stokes equation by application of more subtle inequalities relating to the
bilinear operator B – see section 2.2 in Chapter III of [Tem97]. Other equations arising in dissipative fluid
mechanics can be treated similarly.
Whilst the preceding assumptions on the underlying dynamics apply to a range of interesting models arising
in applications, the following assumptions on the observation model are rather restrictive; however we have been
unable to extend the analysis without making them. We will demonstrate, by means of numerical experiments,
that our results extend beyond the observation scenario employed in the theory
Assumption 2.5. (Observation Model) The system is completely observed so that K = H and H = I.
Furthermore the i.i.d. noise sequence {ξ
j
} is white so that ξ
1
∼ N (0, Γ) with Γ = γ
2
I.
The following consequence of Assumption 2.3 will be useful to us.
Lemma 2.6. Let Assumptions 2.3 hold. Then there is β ∈ R such that, for any v
0
∈ Λ, h > 0 and w
0
∈ H,
|Ψ
h
(v
0
) − Ψ
h
(w
0
)| ≤ e
βh
|v
0
− w
0
| .
3
Proof. Let v, w denote the solutions of (2.2) with initial conditions v
0
, w
0
respectively; define e = v − w. Then
de
dt
+ Ae + 2B(v, e) − B(e, e) = 0 , (2.8)
with e(0) = v
0
− w
0
. Taking the inner-product with e, using (2.3), (2.4) and (2.5), and choosing δ = λ/(2cR),
gives
1
2
d
dt
|e|
2
+ λkek
2
≤ 2ckvkkek|e|
≤ 2cRkek|e|
≤ cR(δkek
2
+ δ
−1
|e|
2
)
=
λ
2
kek
2
+
2
λ
(cR)
2
|e|
2
.
Thus
d
dt
|e|
2
≤
4
λ
(cR)
2
|e|
2
and the desired result follows from an application of the Gronwall inequality.
3 The Ensemble Kalman Filter
3.1 The Algorithm
The idea of the EnKF is to represent the filtering distribution through an ensemble of particles, to propagate
this ensemble under the model to approximate the mapping P(u
j
|Y
j
) to P(u
j+1
|Y
j
) (refered to as prediction
in the applied literature), and to update the ensemble distribution to include the data point Y
j+1
by using a
Gaussian approximation based on the second order statistics of the ensemble (refered to as analysis in the
applied literature).
The prediction step is achieved by simply flowing forward the ensemble under the model dynamics, that is
bv
(k)
j+1
= Ψ(v
(k)
j
) , for k = 1 . . . K.
The analysis step is achieved by performing a randomised version of the Kalman update formula, and using the
empirical covariance of the prediction ensemble to compute the Kalman gain. There are many variants on the
basic EnKF idea and we will study the perturbed observation form of the method.
The algorithm proceeds as follows.
1. Set j = 0 and draw an independent set of samples {v
(k)
0
}
K
k=1
from N(m
0
, C
0
).
2. (Prediction) Let bv
(k)
j+1
= Ψ(v
(k)
j
) and define
b
C
j+1
as the empirical covariance of {bv
(k)
j+1
}
K
k=1
. That is,
b
C
j+1
=
1
K
K
X
k=1
(bv
(k)
j+1
− ¯v
j+1
) ⊗ (bv
(k)
j+1
− ¯v
j+1
) ,
where ¯v
j+1
=
1
K
P
K
k=1
bv
j+1
denotes the ensemble mean.
3. (Observation) Make an observation y
j+1
= Hu
j+1
+ Γ
1/2
ξ
j+1
. Then, for each k = 1 . . . K, generate an
artificial observation
y
(k)
j+1
= y
j+1
+ Γ
1/2
ξ
(k)
j+1
,
where ξ
(k)
j+1
are N(0, I) distributed and pairwise independent.
4
4. (Analysis) Let v
(k)
j+1
be the minimiser of the functional
J(v) =
1
2
|y
(k)
j+1
− v|
2
Γ
+
1
2
|bv
(k)
j+1
− v|
b
C
j+1
.
5. Set j 7→ j + 1 and return to step 2.
The name “perturbed observation EnKF” follows from the construction of the artificial observations y
(k)
j+1
which
are found by perturbing the given observation with additional noise. The sequence of minimisers v
j+1
can be
written down explicitly by simply solving the quadratic minimization problem. This straightforward exercise
yields the following result.
Proposition 3.2. The sequence {v
(k)
j
}
j≥0
is defined by the equation
(I +
b
C
j+1
H
T
Γ
−1
H)v
(k)
j+1
= bv
(k)
j+1
+
b
C
j+1
H
T
Γ
−1
y
(k)
j+1
,
for each k = 1, . . . , K.
Hence, collecting the ingredients from the preceding, the defining equations of the EnKF are given by
(I +
b
C
j+1
H
T
Γ
−1
H)v
(k)
j+1
= Ψ(v
(k)
j+1
) +
b
C
j+1
H
T
Γ
−1
y
(k)
j+1
(3.1a)
y
(k)
j+1
= y
j+1
+ Γ
1/2
ξ
(k)
j+1
(3.1b)
¯v
j+1
=
1
K
K
X
k=1
Ψ(v
(k)
j+1
) (3.1c)
b
C
j+1
=
1
K
K
X
k=1
Ψ(v
(k)
j+1
) − ¯v
j+1
⊗
Ψ(v
(k)
j+1
) − ¯v
j+1
. (3.1d)
There are other representations of the EnKF that are more algorithmically convenient, but the formulae (3.1)
are better suited to our analysis.
3.3 Connection to Randomized Maximum Likelihood
The analysis step of EnKF can be understood in terms of the Randomised Maximum Likelihood (RML) method
widely used in oil reservoir history matching applications [ORL08]. We will now briefly describe this method.
Suppose that we have a random variable u and that u ∼ N( bm,
b
C). Moreover, let G be some linear operator and
suppose we observe
y = Gu + ξ where ξ ∼ N(0, Γ) .
One can use Bayes’ theorem to write down the conditional density P(u|y). In practice however, it is often
sufficient (or sometimes even better) to simply have a collection of samples {u
(k)
}
K
k=1
from the conditional
distribution, rather than the density itself. RML is a method of taking samples from the prior N( bm,
b
C) and
turning them into samples from the posterior. This is achieved as follows, given bu
(k)
∼ N( bm,
b
C) (samples
from the prior), define u
(k)
for each k = 1 . . . K by u
(k)
= argmin
u
J
(k)
(u) where
J
(k)
(u) =
1
2
|y − Gu + Γ
1/2
ξ
(k)
|
2
Γ
+
1
2
|u − bu
(k)
|
2
b
C
,
where ξ
(k)
∼ N(0, I) and independent of ξ. The u
(k)
are then draws from the posterior distribution of u|y
which is a Gaussian with mean m and covariance C. Since one can explicitly write down (m, C), it may be
checked that the u
(k)
defined as above are independent random variables of the form u
(k)
= m + C
1/2
ζ
(k)
,
where ζ
(k)
∼ N (0, I) i.i.d. and are hence draws from the desired posterior, as we know show.
5