Ergodicity and Accuracy of Optimal Particle Filters
for Bayesian Data Assimilation

arXiv:1611.08761v1 [math.PR] 26 Nov 2016

March 7, 2018

D. Kelly⋆ and A.M. Stuart†
⋆
†

New York University, Email: dtbkelly@gmail.com
California Institute of Technology, Email: astuart@caltech.edu

Abstract
Data assimilation refers to the methodology of combining dynamical models and observed data with the objective
of improving state estimation. Most data assimilation algorithms are viewed as approximations of the Bayesian
posterior (filtering distribution) on the signal given the observations. Some of these approximations are controlled,
such as particle filters which may be refined to produce the true filtering distribution in the large particle number
limit, and some are uncontrolled, such as ensemble Kalman filter methods which do not recover the true filtering
distribution in the large ensemble limit. Other data assimilation algorithms, such as cycled 3DVAR methods, may
be thought of as approximating the mean of the posterior, but are also uncontrolled in general. For particle filters
and ensemble Kalman filters it is of practical importance to understand how and why data assimilation methods
can be effective when used with a fixed small number of particles, since for many large-scale applications it
is not practical to deploy algorithms close to the large particle limit asymptotic. In this paper we address this
question for particle filters and, in particular, study their accuracy (in the small noise limit) and ergodicity (for
noisy signal and observation) without appealing to the large particle number limit. We first prove the accuracy
and ergodicity properties for the true filtering distribution, working in the setting of conditional Gaussianity for
the dynamics-observation model. We then show that these properties are inherited by optimal particle filters for
any fixed number of particles. For completeness we also prove large particle number consistency results for the
optimal particle filters, by writing the update equations for the underlying distributions as recursions. In addition to
looking at the optimal particle filter with standard resampling, we derive all the above results for the Gaussianized
optimal particle filter and show that the theoretical properties are favorable when compared to the standard optimal
particle filter.

1 Introduction
1.1 Background and Literature Review
Data assimilation describes the blending of dynamical models with data, with the objective of improving state
estimation and forecasts. The use of data assimilation originated in the geophysical sciences, but is now ubiquitous in engineering and the applied sciences. In numerical weather prediction, large scale ocean-atmosphere
models are assimilated with massive data sets, comprising observational data from satellites, ground based
weather stations and underwater sensors for example [3]. Data assimilation is prevalent in robotics; the SLAM
problem seeks to use sensory data made by robots in an unknown environment to create a map of that environment and locate the robot within it [41]. It is used in modelling of traffic flow [48]. And data assimilation
is being used in bio-medical applications such as glucose-insulin systems [37] and the sleep cycle [38]. These
examples serve to illustrate the growth in the use of the methodology, its breadth of applicability and the very
different levels of fidelity present in the models and the data in these many applications.
Although typical data assimilation problems can be understood from a Bayesian perspective, for non-linear
and potentially high dimensional models it is often infeasible to make useful computations with the posterior. To

1

circumvent this problem, practitioners have developed assimilation methods that approximate the true posterior,
but for which computations are more feasible. In the engineering communities, particle filters have been developed for this purpose, providing empirical approximations of non-Gaussian posteriors [11, 10]. In geoscience
communities, methods are typically built on Kalman filtering theory, after making suitable Gaussian approximations [25]. These include variational methods like 3DVAR and 4DVAR [14, 28], the extended Kalman filter
(ExKF) [12] and the ensemble Kalman filter (EnKF) [5, 13] which relies on a Gaussian ansatz rendering it, in
general, invalid as an approximation of the true filtering distribution [26].
Despite their enormous success, many of these algorithms remain mysterious from a theoretical perspective.
At the heart of the mystery is the fact that data assimilation methods are frequently and successfully implemented in regimes where the approximate filter is not provably valid; it is not known which features of the
posterior (the true filtering distribution) are reflected in the approximation and which are not. For example,
the ensemble Kalman filter is often implemented with ensemble size several orders of magnitude smaller than
needed to reproduce large sample size behaviour, and is applied to problems for which the Gaussian ansatz may
not be valid; it nonetheless can still exhibit skillful state estimates, with high correlations between the estimate
and true trajectories [16, 29]. Indeed, the success of the methods in this non-asymptotic regime is the crux of
their success; the methods would often be computationally intractable at large ensemble sizes.
This lack of theory has motivated recent efforts to better understand the properties of data assimilation methods in the practical, non-asymptotic regimes. The 3DVAR algorithm has been investigated in the context of toy
models for numerical weather prediction, including the Lorenz-63 [24], Lorenz-96 [23] and 2d Navier-Stokes
equations [22]; see also [33]. These works focus primarily on the question of accuracy - how well does the state
estimate track the true underlying signal. Accuracy for the EnKF with fixed ensemble size was first investigated
in [19]. Accuracy for EnKF has been further developed in [30] using linear models with random coefficients, but
much more realistic (practical) assumptions on observations than [19], and moreover focussing on covariance
consistency through the Mahalanobis norm. The articles [42, 43] were the first to investigate the stability of
EnKF with fixed ensemble size, by formulating the filter as a Markov chain and applying coupling techniques.
This line of research has been continued in [8] by framing the EnKF as a McKean-Vlasov system. The limitations of the non-practical regimes have also been investigated; in [17] the authors construct simple dissipative
dynamical models for which the EnKF is shown to be highly unstable with respect to initial perturbations. This
was the first theoretical insight into the frequently observed effect of catastrophic filter divergence [15].
For the nonlinear filtering distribution itself, there has been a great deal of research over the last several
decades, particularly on the question of stability. Conditional ergodicity for the filtering distribution for general
nonlinear hidden Markov models has been investigated in [21] and later refined in [45]. Ergodicity for nonlinear
filters has been discussed in [20, 9, 7] and exponential convergence results were first obtained in [2, 4].

1.2 Our Contributions
For particle filters, much of the theoretical literature focuses on the question of consistency in the large ensemble
limit, that is, does the empirical approximation converge to the true posterior as the number of particles in the
ensemble N approaches infinity. In many high dimensional applications, particle filters are implemented in the
non-asymptotic regime, notably in robotics [41] and ocean-atmosphere forecasting [46]. In geosciences, new
filtering algorithms have been proposed to beat the curse of dimensionality and are implemented with ensemble
sizes many orders of magnitude smaller than the state dimension [47]. In this article, we continue the trend
of analyzing algorithms in practical regimes by focusing on the success of particle filters for fixed ensemble
sizes. In this regime, we ask the question of whether particle filters inherit important long-time properties from
the true posterior distributions. In particular, we address the following questions: if it is known that the true
posterior distribution is accurate and conditionally ergodic, can the same be proved of the approximate filter.
We focus our attention on the optimal particle filter (OPF) [1, 49, 27]. The OPF is a sequential importance
sampling procedure whereby particle updates are proposed using an optimal convex combination of the model
prediction and the observational data at the next time step. For details on the OPF, including the justification for
calling it optimal, see [11]. There are two main reasons that we focus our attentions on the OPF. First, the optimal particle filter is known to compare favorably to the standard particle filter, particularly from the perspective
of weight degeneracy in high dimensions [39, 40]. Indeed the optimal particle filter can be considered a special
case of more complicated filters that have been proposed to beat the curse of dimensionality [6, 47]. Second,
2

under natural assumptions on the dynamics-observation model, the optimal particle filters can be formulated as
a random dynamical system, which is very similar to the 3DVAR algorithm. This means that techniques for
proving accuracy for the 3DVAR filter in earlier literature [35] can be leveraged for the OPF.
Throughout the article, we make the assumption of conditional Gaussianity for the dynamics-observation
model. This framework is frequently employed in practice, particularly in geoscience data assimilation problems. Under this assumption, we show that the true posterior, the filtering distribution, satisfies the long-time
properties of stability and accuracy. The stability result is a type of conditional ergodicity: two copies of the
posterior distribution, initialized at different points but conditioned on the same fixed data, will converge to each
other in the long time limit, with exponential rate. The accuracy result states that, if sufficiently many variables
are observed, the posterior will concentrate around the true trajectory in the long time limit. Related conditional
ergodicity results are obtained under quite general assumptions in [44, 45]. However we believe that the approach we adopt in this article is of considerable value, not only here but also in other potential applications,
due to the simplicity of the underlying coupling argument, the formulation through random dynamical systems
and the explicit link to 3DVAR type algorithms.
We then show that, under the same model-observation assumptions, the OPF also exhibits the long-time
properties of stability and accuracy for any fixed ensemble size. For the conditional ergodicity result, we show
that the two copies of the particle ensembles, initialized differently, but updated with the same observational
data, will converge to each other in the long term limit, in a distributional sense. Both results use very similar
arguments to those employed for the posterior.
In addition, we also establish large ensemble consistency results for the OPF. Here we employ a technique
exposed very clearly in [34], which finds a recursion that is approximately satisfied by the bootstrap particle
filter, and leverages this fact to obtain an estimate on the distance between the true posterior and the empirical
approximation. We show that the same idea can be applied to not only the OPF, but a very large class of
sequential importance sampling procedures. We would like to comment that large particle consistency results
for particle filters should not be considered practical results for high dimensional data assimilation problems,
as in practice particle filters are never implemented in this regime. The consistency results are included here
as they are practically informative for low dimensional data assimilation problems and moreover as the results
are natural consequences of the random dynamical system formulation that has been adopted for accuracy
and ergodicity results. For high dimensional data assimilation problems, it may be more practical to look at
covariance consistency, as done in [30].
In addition to obtaining the results of stability, accuracy and consistency for the OPF, for which we perform resampling at the end of each assimilation cycle, we also prove the corresponding results for the so called
Gaussianized OPF. The Gaussianized OPF was introduced in [18] and differs from the OPF only in the implementation of the resampling. Nevertheless, it was shown numerically in [18] that the GOPF compares favorably
to the OPF, particularly when applied to high dimensional models. The analysis in this article lends theoretical
weight to the advantages of the GOPF over the OPF. In particular we find that the upper bound on the convergence rate for conditional ergodicity for the GOPF has favourable dependence on dimension when compared
with the OPF.

1.3 Structure of Article and Notation
The remainder of the article is structured as follows. At the end of this section we introduce some notation
and terminology that will be useful in the sequel. In Section 2, we formulate the Bayesian problem of data
assimilation, introduce the model-observation assumptions under which we work, and prove the stability and
accuracy results for the true posterior. In Section 3, we introduce the bootstrap particle filter, optimal particle
filter and Gaussianized optimal particle filter. In Section 4, we prove the conditional ergodicity results for the
optimal particle filters. In Section 5, we prove the accuracy results for the optimal particle filters. Finally, in
Section 6, we prove the consistency results for the optimal particle filters.
Throughout we let X denote the finite dimensional Euclidean state space and and we let Y denote the finite
dimensional Euclidean observation space. We write M(X ) for the set of probability measures on X . We
denote the Euclidean norm on X by | · | and for a symmetric positive definite matrix A ∈ L(X , X ), we define
P
| · |A = |A−1/2 · |. We define S N : M(X ) → M(X ) to be the sampling operator S N µ = N1 N
n=1 δu(n) where
3

u(n) ∼ µ are i.i.d. random variables.

2 Bayesian Data Assimilation
We describe the set-up which encompasses all the work in this paper, and then study the conditional ergodicity
and accuracy of the true filtering distribution.

2.1 Set-Up
The state model is taken to be a discrete time Markov chain {uk }k≥0 taking values in the state space X . We
assume that the initial condition u0 of the chain is distributed according to µ0 , where µ0 ∈ M(X ). The
transition kernel for the Markov chain is given by P (uk+1 |uk ). For each k ≥ 1, we make an observation of the
Markov chain
yk+1 = h(uk+1 ) + ηk+1 ,
(2.1)
where h : X → Y maps the state space to the observation space, and ηk ∼ N (0, Γ) are centred i.i.d. random
variables representing observational noise. We denote by Yk = (y1 , . . . , yk ) the accumulated observational data
up to time k. We are interested in studying and approximating the filtering distribution µk (·) = P(uk ∈ ·|Yk )
for all k ≥ 1. We will denote the density of µk by P (uk |Yk ).
The focus of this article is analysis of variants of the optimal particle filter, and in this setting we will always
require the following assumptions on the dynamics-observation model:
Assumption 2.1. The dynamics-observation model is given by
uk+1 = ψ(uk ) + ξk ,
yk+1 = Huk+1 + ηk+1 ,

(2.2a)
(2.2b)

where u0 ∼ µ0 , ξk ∼ N (0, σ 2 Σ0 ) i.i.d. and where ηk ∼ N (0, γ 2 Γ0 ) i.i.d. and σ, γ 6= 0. We require that Σ0
and Γ0 are strictly positive-definite, and that the function ψ(·) is bounded.
In order to facilitate comparison with the standard bootstrap filter, for which the theory we describe is
sometimes derived under different conditions, we retain the general setting of a dynamics given by a general
Markov chain with observations (2.1).
For most of the results in this article, we will be interested in properties of conditional distributions P (uk |Yk )
when the observational data Yk is generated by a fixed realization of the model. For this reason, we introduce
the following notation to emphasize that we are considering a fixed realization of the data, generated by a fixed
trajectory of the underlying dynamical system.
Assumption 2.2. Fix u†0 ∈ X and positive semi-definite matrices Σ∗ and Γ∗ on X and Y respectively. Let
{u†k } be a realization of the dynamics satisfying
u†k+1 = ψ(u†k ) + rγξk†
where u†0 ∈ X is fixed and ξk† ∼ N (0, Σ∗ ) i.i.d. and similarly define {yk† } by
†
†
yk+1
= Hu†k+1 + γηk+1

(2.3)

†
where ηk+1
∼ N (0, Γ∗ ) i.i.d. We will refer to {u†k }k≥0 as the true signal and {yk† }k≥1 as the given fixed data.
As above, we use the shorthand Yk† = {yi† }ki=1 .

Note that this data is not necessarily generated from the same statistical model used to define the filtering
distribution since Σ∗ and Γ∗ may differ from Σ and Γ, and the initial condition is fixed. Note also that we have
set σ = rγ and that, when studying accuracy, we will consider families of data sets and truths parameterized by
γ → 0; in this setting it is natural to think of r and the noise sequences {ξk† }k≥0 and {ηk† }k≥0 as fixed, whilst
the truth and data sequences will depend on the value of γ.
4

Assumption 2.1 renders the conditional distribution P (uk |Yk ) tractable as it has a conditionally Gaussian
structure described by the inhomogeneous Markov process with transition kernel


1
1
2
2
(2.4)
pk+1 (uk , duk+1 ) ∝ exp − |yk+1 − Huk+1 |Γ − |uk+1 − ψ(uk )|Σ duk+1
2
2
and initialized at the measure µ0 . A simple completion of the square yields an alternative representation for the
transition kernel, namely


1
1
2
2
(2.5)
pk+1 (uk , duk+1 ) ∝ exp − |yk+1 − Hψ(uk )|S − |uk+1 − mk+1 |C duk+1 ,
2
2
where

C −1 = Σ−1 + H ∗ Γ−1 H
mk+1

S = HΣH ∗ + Γ

= C Σ−1 ψ(uk ) + H ∗ Γ−1 yk+1 .
,

(2.6)

The conditional mean mk+1 is often given in Kalman filter form

mk+1 = (I − KH)ψ(uk ) + Kyk+1 ,

(2.7)

K = ΣH ∗ S −1 .

(2.8)

where K is the Kalman gain matrix
The expression (2.4) is an application of Bayes’ formula in the form
pk+1 (uk , duk+1 ) ∝ P (yk+1 |uk+1 )P (uk+1 |uk )duk+1
whilst (2.5) follows from a second application of Bayes’ formula to derive the identity
P (yk+1 |uk+1 )P (uk+1 |uk )duk+1 = P (yk+1 |uk )P (uk+1 |uk , yk+1 )duk+1 .
We note that



1
2
P (yk+1 |uk ) =
− |yk+1 − Hψ(uk )|S
2


1
−1
2
P (uk+1 |uk , yk+1 ) = ZC exp − |uk+1 − mk+1 |C .
2
ZS−1 exp

(2.9)

These formulae are prevalent in the data assimilation literature; in particular (2.7) describes the evolution of
the mean state estimate in the cycled 3DVAR algorithm, setting uk+1 = mk+1 [25]. We will make use of the
formulae in section 3 when describing optimal particle filters as random dynamical systems.
We now state two theorems concerning the ergodicity of the filtering distribution itself. The first shows that,
when initialized at two different points, the filtering distributions converge towards one another at a geometric
rate. The second shows that, on average over all instances of randomness in the signal process and the observations process, the mean under the filtering distribution is, asymptotically for large k, O(γ) from the truth. We
state and prove these two results in this section. The remainder of the paper is devoted to establishing analogous
results for various particle filters.

2.2 Conditional Ergodicity
The ergodicity result uses a metric on probability measures to quantify convergence of differently initialized
posteriors in the long time limit. To this end, we define the total variation metric on M(X ) by
dT V (µ, ν) =

1
sup |µ(h) − ν(h)|
2 |h|≤1

where theR supremum is taken over all bounded functions h : X → R with |h| ≤ 1, and where we define
µ(h) := X h(x)µ(dx) for any probability measure µ ∈ M(X ) and any real-valued test function h bounded
by 1 on X .
5

Theorem 2.3. Consider the filtering distributions µk , µ′k , defined by Assumption 2.1, and initialized at two
different Dirac masses µ0 = δz0 and µ0 = δz0′ respectively. Assume moreover that the data appearing in each
filtering distribution is the same, and given by {yk† }k≥1 defined in Assumption 2.2. Then there exists z ∈ (0, 1)
such that, almost surely with respect to the randomness generating the data {yk† }k≥1 ,

1/k
lim dTV µk , µ′k
=z.
(2.10)
k→∞

Proof. The proof uses a standard coupling procedure, as can be found for instance in [32]; we follow the
exposition of the methodology in [31]. We divide the proof into three steps; Step A is where we recast the
problem as a coupling problem, Step B is where we obtain a minorization condition and Step C is where we
apply the law of large numbers to obtain a rate. Subsequent proofs will use the same structure and it will be
useful to refer to specific steps in the proof later on.
Step A: We first introduce two Markov chains zk , zk′ whose laws are given by µk , µ′k respectively. To this
end, define zk to be the inhomogeneous Markov chain, with initialization z0 co-inciding with the center of
the Dirac mass defined above, and with inhomogeneous transition kernel pk defined by (2.4). If we define
recursively the kernels pk (z, ·) by
Z
pk+1 (v, ·)pk (z, dv)
(2.11)
pk+1 (z, ·) =
X

0

with p (z, ·) = δz (·), then it is easy to see that µk = pk (z0 , ·) and µ′k = pk (z0′ , ·) for all k ≥ 0. Since pk (z, ·)
is the law of zk , we see that zk is indeed the required Markov chain, and similarly for zk′ .
The main step is to prove a minorization condition for the transition kernel pk+1 . That is, we seek a measure
Q ∈ M(X ) and a sequence of constants ǫk > 0 satisfying
pk+1 (u, A) ≥ ǫk Q(A)

(2.12)

for all u ∈ X and all measurable sets A ⊂ X . Given a minorization condition, we obtain the result via the
following coupling argument. The minorization condition allows us to define a new Markov kernel

pek+1 (x, A) = (1 − ǫk )−1 pk+1 (x, A) − ǫk Q(A) .
We then define a Markov chain z̃k by z̃0 = z0 and

z̃k+1 ∼ pek+1 (z̃k , ·) w.p. (1 − ǫk )
z̃k+1 ∼ Q(·) w.p. ǫk .

(2.13a)
(2.13b)

and similarly for z̃k′ . By evaluating Eϕ(z̃k ) for a suitable class of test functions ϕ, it is easy to verify that the
Markov chain z̃k is equivalent in law to zk .
We derive the minorization condition which allows this coupling below, in Step B. Here in Step A we couple
the two Markov chains z̃k , z̃k′ in such a way that z̃k = z̃k′ for all k ≥ τ , where τ is the first coupling time, that
is, the smallest n such that z̃n = z̃n′ . Importantly the two chains z̃k , z̃k′ share the same random variables and
live on the same probability space; in particular, once a draw from Q(·) is made, the chains can be coupled
and remain identical thereafter. Let Ak be the event that, in (2.13), the first (state dependent) Markov kernel is
picked at times j = 0, · · · , k − 1. Then since pk (z0 , ·) = µk is the law of z̃k , and similarly for z̃k′ , we have that

 1
sup |E f (z̃k ) − f (z̃k′ ) |
dTV pk (z0 , ·), pk (z0′ , ·) =
2 |f |∞ ≤1


 
1
sup E f (z̃k ) − f (z̃k′ ) IAk + f (z̃k ) − f (z̃k′ ) IAck .
=
2 |f |∞ ≤1
Note that for this coupling the second term vanishes, as in the event Ack , the two chains will have coupled for
τ < k. It follows that

dTV pk (z0 , ·), pk (z0′ , ·) ≤ E(IAk ) = P(Ak ) = Πkj=1 (1 − ǫj ) .
6

To obtain the result (2.10), we need to understand the limiting behaviour of the constants ǫj appearing in the
minorization condition (2.12). Hence we turn our attention toward obtaining the minorization condition.
Step B: If we recall (2.4) and let ρ†k,0 = σHξk† + γηk† , then we obtain, via Z ≤ 1, the identity

 1
1
pk+1 (u, dv) ≥ exp − |Hψ(u†k ) + ρ†k,0 − Hv|2Γ − |v − ψ(u)|2Σ dv
2
 2

† 2
† 2
2
≥ exp −2|Hψ(uk )|Γ − |ψ(u)|Σ − 2|ρk,0 |Γ − |Hv|2Γ − |v|2Σ dv



 1
≥ exp −λ2 − 2|ρ†k,0 |2Γ exp − |v|2D dv
2
where

(2.14)
(2.15)
(2.16)
(2.17)



λ2 = sup 2|Hψ(v)|2Γ + |ψ(u)|2Σ
u,v

and

1 −1
D = Σ−1 + H ∗ Γ−1 H .
2
Thus we have a minorization condition of the form

pk+1 (u, A) ≥ ǫk Q(A)


where Q(·) is the Gaussian N (0, D) and ǫk = exp −λ2 − 2|ρ†k,0 |2Γ .
Step C: By the argument above, it follows that
1/k
dTV pk (z0 , ·), pk (z0′ , ·)
≤ zk


1/k
. Since the ǫk are i.i.d. and integrable, by the law of large numbers, almost
where zk = Πkj=1 (1 − ǫj )
surely with respect to the randomness generating the true signal and the data, we have
ln zk =

∞
k
X
1X
1 n
ln(1 − ǫj ) → E ln(1 − ǫ1 ) = −E
ǫ1 .
k j=1
n
n=1

But ǫ1 = c exp(−2|ρ†1,0 |2Γ ) for some c ∈ (0, 1]. Since ρ†1,0 is Gaussian it follows that the nth moment of ǫ1
1
scales like n− 2 so that the limit of ln zk is negative and finite; the result follows.

2.3 Accuracy
We now discuss accuracy of the posterior filtering distribution in the small noise limit γ ≪ 1.
Assumption 2.4. Let r = σ/γ and assume that there is rc > 0 such that, for all r ∈ [0, rc ), the function
(I − KH)ψ(·), with K defined through (2.6) and (2.8), is globally Lipschitz on X with respect to the norm k · k
and with constant α = α(r) < 1.
Theorem 2.5. Suppose Assumption 2.4 holds for some rc > 0. Then for all r ∈ [0, rc ) we have
lim sup Ekuk − Eµk uk k2 ≤ cγ 2 ,
k→∞

where Eµk denotes expectation over the posterior µk defined through Assumption 2.1 and E denotes expectation
over the dynamical model and the observational data.

7

Proof. This follows similarly to Corollary 4.3 in [36], using the fact that the mean of the filtering distribution is
optimal in the sense that
Ekuk − Eµk uk k2 ≤ Ekuk − mk k2
for any Yk −measurable sequence {mk }. We use for mk the 3DVAR filter
mk+1 = (I − KH)ψ(mk ) + Kyk+1 .
Let ek = uk − mk . Following closely Theorem 4.10 of [25] we obtain
Ekek+1 k2k+1 ≤ α2 Ekek k2 + O(γ 2 ).
Application of the Gronwall lemma, plus use of the optimality property, gives the required bound.

3 Particle Filters With Resampling
In this section we introduce the bootstrap particle filter, and the optimal particle filters, in all three cases with
resampling at every step. Assumption 2.1 ensures that the three particle filters have an elegant interpretation
as a random dynamical system (RDS) which, in addition, is useful for our analyses. We thus introduce the
filters in this way before giving the algorithmic definition which is more commonly found in the literature. The
bootstrap particle filter will not be the focus of subsequent theory, but does serve as an important motivation for
the optimal particle filters, and in particular for the consistency results in Section 6.
For each of the three particle filters, we will make frequent use of a resampling operator, which draws a
(n)
(m)
(m) N
sample uk from {b
uk }N
m=1 with weights {wk }m=1 which sum to one. To define this operator, we define
(m)
(m)
(m+1)
(m+1)
(m)
(m)
the intervals Ik = [αk , αk
) where αk
= αk + wk and then set
(n)

uk

=

N
X

(n)

(m)

II (m) (rk )b
uk

m=1

(3.1)

k

(n)

(n)

where rk ∼ U (0, 1) i.i.d. . Since the weights sum to one, rk
(n)
(i )
we will have uk = u
bk ∗ . We also notice that

(i )

will lie in exactly one of the intervals Ik ∗ and

N
N
X
X
1
(m)
δu(m) = S N
wk δub(m)
k
k
N
m=1
m=1

where S N is the sampling operator defined previously.

3.1 The Bootstrap Particle Filter
The bootstrap particle filter (BPF) approximates the filtering distribution µk with an empirical measure
ρN
k =

N
X
1
δu(n) .
k
N
n=1

(3.2)

(n)

The particle positions {uk }N
n=1 are defined as follows.
(n)

(n)

(n)

,

(n)

(m)

u
bk+1 = ψ(uk ) + ξk
(n)

uk+1 =

N
X

m=1

(n)

ξk

II (m) (rk+1 )b
uk+1 ,
k+1

8

∼ N (0, Σ) i.i.d. ,
(3.3)

where the second equation uses the resampling operator defined in (3.1) with weights computed according to
1
(n)
(n),∗
bk+1 |2Γ ) ,
wk+1 = exp(− |yk+1 − H u
2

(n),∗

w
(n)
wk+1 = PN k+1(j),∗ .
j=1 wk+1

(3.4)

Thus, for each particle in the RDS, we propagate them forward using the dynamical model and then re-sample
from the weighted particles to account for the observation likelihood.
As with most particle filters, the motivation for the bootstrap particle filter stems from an importance sampling scheme applied to a particular decomposition of the filtering distribution. By Bayes formula, we have
P (uk+1 |Yk+1 ) = P (uk+1 |Yk , yk+1 )
1
=
P (yk+1 |uk+1 , Yk )P (uk+1 |Yk )
P (yk+1 |Yk )
Z
1
P (yk+1 |uk+1 )P (uk+1 |uk , Yk )P (uk |Yk )duk
=
P (yk+1 |Yk )
Z
1
P (yk+1 |uk+1 )P (uk+1 |uk )P (uk |Yk )duk
=
P (yk+1 |Yk )
Z
1
=
P (yk+1 |uk+1 ) P (uk+1 |uk )P (uk |Yk )duk
P (yk+1 |Yk )

(3.5)

The bootstrap particleR filter approximates the posterior via a sequential application of importance sampling,
using P (uk+1 |Yk ) = P (uk+1 |uk )P (uk |Yk )duk as the proposal and re-weighting according to the likelihood
P (yk+1 |uk+1 ). Thus the method is typically described by the following algorithm for updating the particle
(n)
positions. The particles are initialized with u0 ∼ µ0 and then updated by
(n)

(n)

1. Draw u
bk+1 ∼ P (uk+1 |uk ).
(n)

2. Define the weights wk+1 for n = 1, . . . , N by
(n),∗

(n),∗

(n)

wk+1 = P (yk+1 |b
uk+1 ) ,
(n)

(n)

w
(n)
wk+1 = PN k+1(m),∗ .
m=1 wk+1

(n)

N
3. Draw uk+1 from {b
uk+1 }N
n=1 with weights {wk+1 }n=1 .

Under Assumption 2.1, it is clear that the sampling and re-weighting procedures are consistent with (3.3). Note
that the normalization factor P (yk+1 |Yk ) is not required in the algorithm and is instead approximated via the
normalization procedure in the second step.
It is useful to define the related measure
ρbN
k =

N
X

(n)

wk δub(n) ,

(3.6)

k

n=1

N
N N
with ρbN
bk . As we shall see in Section 6, the
0 = µ0 , which is related to the bootstrap particle filter by ρk = S ρ
N
advantage of ρbk is that it has a recursive definition which allows for elegant proofs of consistency results [34].

3.2 Optimal Particle Filter

The optimal particle filter with resampling can also be defined as a RDS. We once again approximate the filtering
distribution µk with an empirical distribution
µN
k =

N
X
1
δu(n) .
k
N
n=1

9

(3.7)

Under Assumption 2.1 the particles in this approximation are defined as follows. The particle positions are
(n)
initialized with u0 ∼ µ0 and then evolve according to the RDS
(n)

(n)

(n)

u
bk+1 = (I − KH)ψ(uk ) + Kyk+1 + ζk
(n)

uk+1 =

N
X

m=1

(n)

,

(n)

ζk

∼ N (0, C) i.i.d.
(3.8)

(m)

II (m) (rk+1 )b
uk+1 .
k+1

where C, S, K are defined in (2.6), (2.7) and as with the BPF, the second equation uses the resampling operator
defined in (3.1) but now using weights computed by
(n),∗
wk+1

1
(n)
= exp(− |yk+1 − Hψ(uk )|2S ) ,
2

(n),∗

(n)
wk+1

wk+1

= PN

m=1

(m),∗

wk+1

.

(3.9)

In light of the formulae given in (2.9), which are derived under Assumption 2.1, we see that the optimal particle
(n)
filter is updating the particle positions by sampling from P (uk+1 |uk , yk+1 ) and then re-sampling to account
(n)
for the likelihood factor P (yk+1 |uk ). In particular, without necessarily making Assumption 2.1, the optimal
particle filter is a sequential importance sampling scheme applied to the following decomposition of the filtering
distribution
Z
P (uk+1 |Yk+1 ) =
P (uk+1 , uk |Yk+1 )duk
ZX
=
P (uk+1 |uk , yk+1 )P (uk |Yk+1 )duk
(3.10)
ZX
P (yk+1 |uk )
P (uk+1 |uk , yk+1 )P (uk |Yk )duk .
=
X P (yk+1 |Yk )
(n)

In the algorithmic setting, the filter is initialized with u0
(n)

∼ µ0 , then for k ≥ 0

(n)

1. Draw u
bk+1 from P (uk+1 |uk , yk+1 )
(n)

2. Define the weights wk+1 for n = 1, . . . , N by
(n),∗

(n),∗
wk+1

=

(n)
P (yk+1 |b
uk )

,

(n)
wk+1

wk+1

= PN

m=1

(n)

(m)

(m),∗

wk+1

.

(m)

N
3. Draw uk+1 from {b
uk+1 }N
m=1 with weights {wk+1 }m=1 .

It is important to note that, although the OPF is well defined in this general setting for any choice of dynamicsobsevation model, it is only implementable under stringent assumptions on the forward and observation model,
such as those given in Assumption 2.1; under this assumption the step 1 may be implemented using the formulae
given in (2.9) and exploited in the derivation of (3.8). However we emphasize that models satisfying Assumption
2.1 do arise frequently in practice.
As with the BPF, it is beneficial to consider the related particle filter given by
µ
bN
k =

N
X

(n)

wk δub(n)

n=1

k

N
N N
for k ≥ 1 and with µ
bN
bk .
0 = µ0 . Clearly, we have that µk = S µ

10

(3.11)

3.3 Gaussianized Optimal Particle Filter
In [18], an alternative implementation of the OPF is investigated and found to have superior performance on
a range of test problems, particularly with respect to the curse of dimensionality. We refer to this filter as
the Gaussianized optimal particle filter (GOPF). Once again, we approximate the filtering distribution with an
empirical measure
N
X
1
δv(n) .
(3.12)
νkN =
k
N
n=1
As in the previous subsection, we first describe the filter under Assumption 2.1. The filter is initialized with
(n)
v0 ∼ µ0 , with subsequent iterates generated by the RDS
(n)

ṽk

=

N
X

(n)

,

(3.13)

k+1

m=1
(n)
vk+1

(m)

II (m) (rk+1 )vk

= (I −

(n)
KH)ψ(ṽk )

+ Kyk+1 +

(n)
ζk

,

(n)
ζk

∼ N (0, C) i.i.d.

and the weights appearing in the resampling operator are given by
1
(n),∗
(n)
wk+1 = exp(− |yk+1 − Hψ(vk )|2S ) ,
2

(n),∗

w
(n)
wk+1 = PN k+1(m),∗ .
m=1 wk+1

(3.14)

Thus, the update procedure for GOPF is weight-resample-propagate, as opposed to propagate-weight-resample
for the OPF. Hence the only difference between the OPF and GOPF is the implementation of the resampling
step.
In our analysis it is sometimes useful to consider the equivalent RDS
(m,n)

(m)

(m,n)

vbk+1 = (I − KH)ψ(vk ) + Kyk+1 + ζk
(n)

vk+1 =

N
X

m=1

(n)

(m,n)

ζk

∼ N (0, C) i.i.d.
(3.15)

(m,n)

II (m) (rk+1 )b
vk+1 .
k+1

(n)

The sequences vk defined in (3.13) and (3.15) agree because for every n there is exactly one m = m∗ (n) such
(m∗ (n),n)
that vbk+1
survives the resampling step. Writing it this way allows certain parts of our subsequent analysis
to be performed very similarly for both the OPF and GOPF.
For a general dynamics-observation model, the GOPF is described by the algorithm
(n)

1. Define the weights wk+1 for n = 1, . . . , N by
(n),∗

(n),∗
wk+1

=

(n)
P (yk+1 |vk )

,

(n)
wk+1

wk+1

= PN

m=1

(n)

2. Draw ṽk

(n)

(m)

(m),∗

wk+1

.

(m)

N
from {vk }N
m=1 with weights {wk+1 }m=1 .
(n)

3. Draw vk+1 from P (uk+1 |ṽk , yk+1 ) .
Unlike for the previous filters, there is no need to define an associated ‘hatted’ measure, as the GOPF can be
shown to satisfy a very natural recursion. This will be discussed in Section 6.

4 Ergodicity for Optimal Particle Filters
In this section we study the conditional ergodicity of the two optimal particle filters. The proofs are structurally
very similar to that for the filtering distribution itself.
11

4.1 Optimal Filter
(1)

(N )

Before stating the conditional ergodicity result, we first need some notation. Define uk = (uk , . . . , uk ) to
(1)′
(N )′
be particle positions defined by the RDS (3.8) with µ0 = δz0 and similarly with u′k = (uk , . . . , uk ) with
µ0 = δz0′ . Then uk is a Markov chain taking values on X N , whose Markov kernel we denote by qk (z, ·). As in
the proof of Theorem 2.3, the law of uk is given by q k (z0 , ·), defined recursively as in (2.11) and similarly the
law of u′k is given by q k (z0′ , ·). The conditional ergodicity result states that if the two filters uk , u′k are driven
by the same observational data, then the law of uk will converge to the law of u′k exponentially as k → ∞. 1
Theorem 4.1. Consider the OPF particles uk , u′k defined above. Assume moreover that the observational data
used to define each filter is the same, and given by {yk† }k≥1 from Assumption 2.2. Then there exists zN ∈ (0, 1)
such that, almost surely with respect to the randomness generating {yk† }k≥1 ,

1/k
= zN .
lim dTV q k (z0 , ·), q k (z0′ , ·)

(4.1)

k→∞

Proof. We will follow the proof of Theorem 2.3 closely, via Steps A, B and C. Step A proceeds exactly as in
the proof of Theorem 2.3. In particular, if we assume a minorization condition
qk+1 (u, ·) ≥ ǫk Q(·) ,

(4.2)

and repeat the coupling argument from Step A, with uk , u′k in place of zk , zk′ , the we obtain the bound

dTV q k (z0 , ·), q k (z0′ , ·) ≤ Πkj=1 (1 − ǫj ) .

This is indeed how the proof proceeds, with the caveat that the minorization constants ǫk depend on the number
(1)
(N )
of particles N . In Step B, we obtain the minorization condition (4.2). Let u
bk = (b
uk , . . . , u
bk ) as defined in
the RDS (3.8). Before proceeding, we introduce some preliminaries. Using the fact that
†
†
yk+1
= Hψ(u†k ) + γ(rHξk† + ηk+1
)

and defining

N

N
(n)
(n)
ak = (I − KH)ψ(uk ) + KHψ(u†k )
ζk = ζk
n=1
n=1

N
†
†
ρ†k,0 = γ(rHξk† + ηk+1
), ρ†k = γK(rHξk† + ηk+1
)
n=1

we see that

u
bk+1 = ak + ρ†k + ζk .

The next element of the sequence, uk+1 , is then defined by the second identity in (3.8). We are interested in
† ∞
the conditional ergodicity of {uk }∞
k=1 with the sequence {ρk }k=1 fixed. By Assumption 2.1, ak is bounded
uniformly in k. We define the covariance operator C ∈ L(X N , X N ) to be a block diagonal covariance with
each diagonal entry equal to C and then


R = sup |(I − KH)ψ(u) + KHψ(v)|2C ,
(u,v)

which is finite by Assumption 2.1.
Now, let E0 be the event that, upon resampling, every particle survives the resampling. There are N ! such
permutations. We will do the calculation in the case of a trivial permutation, that is, where each particle is
mapped to itself under the resampling. However the bounds which follow work for any permutation of this
1 We abuse notation in this subsection by using u ∈ X N to denote the N particles comprising the optimal particle filter; this differs
k
from the notation uk ∈ X used in the remainder of the paper to denote the underlying dynamical model.

12

because we do not use any information about location of the mean; we simply use bounds on the drift ψ. If each
(n)
(n)
particle is mapped to itself, then uk+1 = u
bk+1 for all n = 1, . . . , N . It follows that
qk+1 (u, A) = P(uk+1 ∈ A|uk = u)

≥ P(uk+1 ∈ A|uk = u, E0 )P(E0 )
= P(b
uk+1 ∈ A|uk = u)P(E0 ) .
We first note that
Z

 1
1
P(b
uk+1 ∈ A|uk = u) = p
exp − |x − ak − ρ†k |2C dx
2
(2π)dN det C A




exp −|ak + ρ†k |2C Z
≥ p
exp −|x|2C dx
(2π)dN det C
A
≥ 2−dN/2 exp(−2|ak |2C ) exp(−2|ρ†k |2C )QC (A)
≥ 2−dN/2 exp(−2N R2 ) exp(−2|ρ†k |2C )QC (A)
where QC (A) is the Gaussian measure N (0, 21 C). Thus we have shown that
P(b
uk+1 ∈ A|uk = u) ≥ δk QC (A)
where

(4.3)

δk = 2−dN/2 exp(−2N R2 ) exp(−2|ρ†k |2C ) .

Moreover, we have that

(n)

P(E0 ) = N !ΠN
n=1 wk+1 .
(n)

(n),∗

(m),∗

Note that we have the bound wk+1 ≥ wk+1 /N for each n = 1, . . . , N because each wk+1 is bounded by 1.
But we have

 1
(n),∗
(n)
wk+1 = exp − |yk+1 − Hψ(uk )|2S
2
 1

(n)
= exp − |Hψ(u†k ) − Hψ(uk ) + ρ†k,0 |2S
 2

≥ exp −r2 − |ρ†k,0 |2S
where
r2 = sup |Hψ(u) − Hψ(v)|2S
u,v

which is finite by Assumption 2.1. From this we see that
P(E0 ) ≥ N !


1
exp −N r2 − N |ρ†k,0 |2S .
NN

Thus we obtain the minorization conditon (4.2) where
ǫk = N !


1
exp −N r2 − N |ρ†k,0 |2S δk
N
N

,

Q = QC .

Finally, Step C follows identically to the proof of Theorem 2.3 and the proof is complete.

13

4.2 Gaussianized Optimal Filter
(1)

(N )

As in the last section, we define vk = (vk , . . . , vk ) and similarly for vk′ using the RDS but now for the
GOPF (3.13) (or alternatively (3.15)) with distinct initializations µ0 = δz0 and µ0 = δz0′ . As with Theorem 4.1,
we let q k (z0 , ·) denote the law of vk .
Theorem 4.2. Consider the GOPF particles vk , vk′ defined above. Assume moreover that the observational data
used to define each filter is the same, and given by {yk† }k≥1 from Assumption 2.2. Then there exists zN ∈ (0, 1)
such that, almost surely with respect to the randomness generating {yk† }k≥1 ,

1/k
lim dTV q k (z0 , ·), q k (z0′ , ·)
= zN .

(4.4)

k→∞

Proof. The proof follows similarly to that of Theorem 4.1, in particular it suffices to obtain a minorization
condition for qk+1 (v, ·). We will use the RDS representation (3.15), which we now recall
(m,n)

(m)

(m,n)

vbk+1 = (I − KH)ψ(vk ) + Kyk+1 + ζk
(n)

vk+1 =

N
X

m=1

(n)

(m,n)

ζk

∼ N (0, C) i.i.d.
(4.5)

(m,n)

vk+1 .
II (m) (rk+1 )b
k+1

(n)

In this formulation, note that for each n there is one and only one m = m∗ (n) such that II (m) (rk ) = 1. We
k
see that
∗
(m (n),n) N
(n)
)n=1 .
vk
vk := (vk )N
n=1 = (b
Using the fact that

†
†
yk+1
= Hψ(u†k ) + γ(rHξk† + ηk+1
)

and defining

N

N
(m∗ (n),n)
(m∗ (n))
ζk = ζk
ak = (I − KH)ψ(vk
) + KHψ(u†k )
n=1
n=1

N
†
†
ρ†k,0 = γ(rHξk† + ηk+1
), ρ†k = γK(rHξk† + ηk+1
)
n=1

we see that

vk+1 = ak + ρ†k + ζk .

Now notice that



(m∗ (n),n) N
qk+1 (v, A) = P(vk+1 ∈ A|vk = v) = P (b
vk+1
)n=1 ∈ A|vk = v
Z
 1

1
p
=
exp − |x − ak − ρ†k |2C dx
2
(2π)dN det C A




exp −|ak + ρ†k |2C Z
≥ p
exp −|x|2C dx
(2π)d det C
A
≥ 2−dN/2 exp(−2|ak |2C ) exp(−2|ρ†k |2C )QC (A)

≥ 2−dN/2 exp(−2N R2 ) exp(−2|ρ†k |2C )QC (A)
where QC is the Gaussian measure N (0, 12 C). Thus we have shown that
qk+1 (v, A) ≥ δk QC (A)
where

δk = 2−dN/2 exp(−2N R2 ) exp(−2|ρ†k |2C ) .

The remainder of the proof (step C) follows identically to Theorem 2.3, Theorem 4.1.
14

(4.6)

Remark 4.3. We can compare our (upper bounds on the) rates of convergence for the true filtering distribution,
and the two optimal filters, using the minorization constants. From the proof of Theorem 2.3, we see that the
minorization constants determine the rate of convergence in the statement of conditional ergodicity. In particular
for k large we have that
dT V (µk , µ′k ) . zk .
where
log z = E log(1 − ǫ1 ) .
The corresponding statements for the law of the particles in each optimal particle filter hold similarly, but ǫ1
depends on the number of particles. For the OPF we have
ǫ1 = N !
where


1
exp −N r2 − N |ρ†1,0 |2S δ1 ,
N
N

δ1 = 2−dN/2 exp(−2N R2 ) exp(−2|ρ†1 |2C ) ;

for the GOPF we simply have ǫ1 = δ1 . The extra N dependence in the OPF clearly leads to a slower (upper
bound on the) rate of convergence for the OPF. Thus, by this simple argument, we obtain a better convergence
rate for the GOPF than for the OPF. This suggests that the GOPF may have a better rate of convergence for fixed
ensemble sizes; further analysis or experimental study of this point would be of interest.

5 Accuracy for Optimal Particle Filters
In this section we study the accuracy of the optimal particle filters, in the small noise limit γ → 0. The
expectation appearing in the theorem statements is with respect to the noise generating the data, and with
respect to the randomness within the particle filter itself. Note that this situation differs from that in the accruacy
result for the filter itself which uses data generated by the statistical model itself. Assumption 2.2 relaxes this
assumption.

5.1 Optimal Particle Filter
(n)

Theorem 5.1. Let Assumption 2.4 hold and consider the OPF with particles {uk }N
n=1 defined by (3.8) with
data {yk† } from Assumption 2.2. It follows that there is constant c such that


(n)
lim sup E max kuk − u†k k2 ≤ cγ 2 .
n

k→∞

Proof. First recall the notation Σ = σΣ0 , Γ = γΓ0 and r = σ/γ. Now define
S0 = r2 HΣ0 H ∗ + Γ0 ,
C0 = r2 (I − KH)Σ0
and note that
S = γ 2 S0 , C = γ 2 C0 , K = r2 Σ0 H ∗ S0−1 .
We will use the RDS representation
(n)

(n)

(n)

†
u
bk+1 = (I − KH)ψ(uk ) + Kyk+1
+ γζ0,k
(n)

uk+1 =

N
X

m=1

(n)

,

(n)

ζk

(n)

II (m) (rk+1 )b
uk+1 .
k+1

(n)

where ζ0,k ∼ N (0, C0 ) i.i.d. . Hence we have
u†k+1 = (I − KH)ψ(u†k ) + KHψ(u†k ) + rγξk† ,
15

∼ N (0, C) i.i.d.
(5.1)

(n)

(n)

(n)

(n)

†
u
bk+1 = (I − KH)ψ(uk ) + K(Hψ(u†k ) + γηk+1
) + γζ0,k .

(n)


(n)
(n)
(n)
u
bk+1 − u†k+1 = (I − KH) ψ(uk ) − ψ(u†k ) + γιk

where ζ0,k ∼ N (0, C0 ) i.i.d. Subtracting, we obtain

where ιk

(5.2)

(n)

†
:= (Kηk+1
+ ζ0,k − rξk† ). Moreover we have the identity

u†k+1 =

N
X

(n)

II (m) (rk+1 )u†k+1 .

(5.3)

k+1

m=1

Thus, defining
(n)

ek

(n)

(n)

= uk − u†k ,

ebk

we have from (5.1) and (5.3)
(n)

ek+1 =

N
X

(n)

(m)

II (m) (rk+1 )b
ek+1 .

m=1

Thus

(n)

=u
bk − u†k

k+1

(n)

(m)

max kek+1 k2 ≤ max kb
ek+1 k2 ,
n

m

where the norm is the one in which we have a contraction. Using (5.2), the Lipschitz property of (I − KH)ψ(·),
taking expectations and using independence, yields




(n)
(n)
E max kuk+1 − u†k+1 k2 ≤ α2 E max kuk − u†k k2 + O(γ 2 )
n

n

and the result follows by Gronwall.

5.2 Gaussianized Optimal Filter
(n)

Theorem 5.2. Let Assumption 2.4 hold and consider the GOPF with particles {vk }N
n=1 defined by (3.13) (or
(3.15)) with data {yk† } from Assumption 2.2. It follows that there is constant c such that


(n)
lim sup E max kvk − u†k k2 ≤ cγ 2 .
k→∞

n

Proof. Recall the notation defined at the beginning of the proof of Theorem 5.1. Recall also the RDS representation of the GOPF (3.15)
(m,n)

(m)

(m,n)

vbk+1 = (I − KH)ψ(vk ) + Kyk+1 + γζ0,k
(n)

vk+1 =

N
X

m=1
(m,n)

where we now have ζ0,k

(n)

(m,n)

II (m) (rk+1 )b
vk+1 .

(5.4)

k+1

∼ N (0, C0 ) i.i.d. , recalling that C = γ 2 C0 . We also have the identity
u†k+1 = (I − KH)ψ(u†k ) + KHψ(u†k ) + rγξk† ,

(5.5)

(n)

where ζ0,k ∼ N (0, C0 ) i.i.d. Subtracting, we obtain
(m,n)

(m)

vbk+1 − u†k+1 = (I − KH) ψ(vk
16


(n)
) − ψ(u†k ) + γιk

(5.6)

(n)

where ιk

(n)

†
:= (Kηk+1
+ ζ0,k − rξk† ). Note that

u†k+1 =

N
X

(n)

II (m) (rk+1 )u†k+1

m=1

(5.7)

k+1

so that, defining
(n)

ek

(n)

= vk

− u†k ,

we have from (5.4), (5.6) and (5.7)
(n)

ek+1 =

N
X

m=1

Thus

(m,n)

ebk

(m,n)

= vbk

(n)

− u†k

(m,n)

II (m) (rk+1 )b
ek+1 .
k+1

(n)

(m,n)

max kek+1 k2 ≤ max kb
ek+1 k2 ,
n

m

where the norm is the one in which we have a contraction. Using (5.6), using the Lipschitz property of (I −
KH)ψ(·), taking expectations and using independence, gives




(n)
(n)
E max kvk+1 − u†k+1 k2 ≤ α2 E max kvk − u†k k2 + O(γ 2 ) .
n

n

The result follows by Gronwall.

6 Consistency in the Large Particle Limit
In this section we state and prove consistency results for the BPF, OPF and GOPF introduced in Section 3. For
the BPF the result is well known but we reproduce it here as it serves as an ideological template for the more
complicated proofs to follow; furthermore we present the clean proof given in [34] (see also [25, Chapter 4]) as
this particular approach to the result generalizes naturally to the OPF and GOPF.

6.1 Bootstrap Particle Filter
In the following, we let fk+1 : X → R be any function with fk+1 (uk+1 ) ∝ P (yk+1 |uk+1 ); any proportionality
constant will suffice, but the normalization constant is of course natural. As in previous sections, we let µk
denote the filtering distribution. The following theorem is stated and then proved through a sequence of lemmas
in the remainder of the subsection.
N
Theorem 6.1. Let ρbN
k , ρk be the BPFs defined by (3.6), (3.2) respectively, and suppose that there exists a
constant κ ∈ (0, 1] such that
κ ≤ fk+1 (uk+1 ) ≤ κ−1
(6.1)

for all uk+1 ∈ X , yk+1 ∈ Y and k ∈ {0, . . . , K − 1}. Then we have
d(b
ρN
K , µK ) ≤

K
X

(2κ−2 )k N −1/2

(6.2)

K
X

(2κ−2 )k N −1/2

(6.3)

k=1

and
d(ρN
K , µK ) ≤

k=0

for all K, N ≥ 1.
Remark 6.2. Note that the constant κ−2 appearing in the estimates above arises as the ratio of the upper and
lower bounds in (6.1). In particular, we cannot optimize κ by choosing a different proportionality constant for
fk+1 .
17

It is straightforward to check that the filtering distribution satisfies the recursion
µk+1 = Lk+1 P µk

(6.4)

as this is nothing more than a statement of the Bayes formula calculation given in (3.5). Here P : M(X ) →
M(X ) is the Markov semigroup
Z
P (uk+1 |uk )ν(duk ),

P ν(A) =

A

which gives the prior, and Lk+1 : M(X ) → M(X ) is the likelihood operator
Z
Lk+1 ν(A) =
Z −1 P (yk+1 |uk+1 )ν(duk+1 )
A

for each A ⊂ X measurable and Z is a normalization constant.
In terms of understanding the approximation properties of the BPF, the key observation is that the measures
{b
ρN
k }k≥0 satisfy the recursion
N
ρbN
bN
ρbN
(6.5)
k
0 = µ0 .
k+1 = Lk+1 S P ρ

where P : M(X ) → M(X ) is the Markov semigroup and, as defined in Section 1.3, S N : M(X ) → M(X )
is the sampling operator. The convergence of the measures is quantified by the metric on random elements of
M(X ) defined by
p
d(µ, ν) = sup|f |∞ ≤1 Eω |µ(f ) − ν(f )|2 ,

where, in our setting, Eω will always denote expectation with respect to the randomness in the sampling operator
S N ; this metric reduces to twice the total variation metric, used in studying ergodicity, when the measures are
not random. The main ingredients for the proof are the following three estimates for the operators P, S N and
Lk+1 with respect to the metric d.
Lemma 6.3. We have the following
1. supν∈M(X ) d(S N ν, ν) ≤ N −1/2 .
2. d(P µ, P ν) ≤ d(µ, ν) for all µ, ν ∈ M(X ).
Proof. See [25, Lemma 4.7, Lemma 4.8].
We state the following Lemma in a slightly more general form than necessary for the BPF, as it will be
applied in different contexts for the optimal particle filters.
Lemma 6.4. Let Z be a finite dimensional Euclidean space. Suppose that gk+1 : Z → [0, ∞) is bounded and
that there exists κ ∈ (0, 1] such that
κ ≤ gk+1 (u) ≤ κ−1
(6.6)
for all u ∈ Z and define Gk+1 : M(Z) → M(Z) by Gk+1 (ν)(ϕ) = ν(gk+1 ϕ)/ν(gk+1 ). Then
d(Gk+1 µ, Gk+1 ν) ≤ (2κ−2 )d(µ, ν)
for all µ, ν ∈ M(Z).
Proof. See [25, Lemma 4.9].
We can now prove the consistency result.

Proof of Theorem 6.1. First note that, taking Z = X and gk+1 = fk+1 in Lemma 6.4, we obtain Gk+1 ν =
Lk+1 ν. Thus, by (6.1), it follows that d(Lk+1 µ, Lk+1 ν) ≤ (2κ−2 )d(µ, ν) for all µ, ν ∈ M(X ). Combining
this fact with the recursions given in (6.5), (6.4) and the estimates given in Lemmas 6.3 we have
N
d(b
ρN
bN
k+1 , µk+1 ) = d(Lk+1 S P ρ
k , Lk+1 P µk )

18

≤ 2κ−2 d(S N P ρbN
k , P µk )

N
−2
d(S P ρbN
bN
bN
≤ 2κ
k ,Pρ
k ) + d(P ρ
k , P µk )

≤ 2κ−2 N −1/2 + 2κ−2 d(b
ρN
k , µk ) .

N N
And since ρbN
bk
0 = µ0 , we obtain (6.2) by induction. Moreover, since ρk = S ρ

and (6.3) follows.

N N
d(ρk , µk ) = d(S N ρbN
bk , ρbN
ρN
k ) + d(b
k , µk )
k , µk ) ≤ d(S ρ

6.2 Sequential Importance Resampler
In this section we will apply the above strategy to prove the corresponding consistency result for the OPF.
Instead of restricting to the OPF, we will obtain results for the sequential importance resampler (SIR), for which
the OPF is a special case. See [11, sections II, III] for background in sequential importance sampling, and on
the use of resampling. As with the OPF, the SIR is an empirical measure
µN
k =

N
X
1
δ (n) .
N uk
n=1

(6.7)

We will abuse notation slightly by keeping the same notation for the OPF and the SIR. The particle positions
are drawn from a proposal distribution π(uk+1 |uk , yk+1 ) and re-weighted accordingly. As usual, the positions
(n)
are initialized with u0 ∼ µ0 and updated by
(n)

(n)

1. Draw u
bk+1 from π(uk+1 |uk , yk+1 )
(n)

2. Define the weights wk+1 for n = 1, . . . , N by
(n)

(n),∗
wk+1

(n)

=

(n)

(n),∗

(n)

P (yk+1 |b
uk+1 )P (b
uk+1 |uk )
(n)

,

(n)

π(b
uk+1 |uk , yk+1 )

(m)

(n)
wk+1

wk+1

= PN

m=1

(m),∗

wk+1

.

(m)

N
3. Draw uk+1 from {b
uk+1 }N
m=1 with weights {wk+1 }m=1 .

Thus, if we take the proposal to be π(uk+1 |uk , yk+1 ) = P (uk+1 |uk , yk+1 ) then we obtain the OPF (3.7).
Without being more specific about the proposal π, it is not possible to represent the SIR as a random dynamical
system in general.
Precisely as with the OPF, for the SIR we define the related filter
µ
bN
k =

N
X

(n)

wk δub(n)

n=1

(6.8)

k

N N
N
bk . The following theorem, and corollary, are proved
with µ
bN
0 = µ0 and note the important identity µk = S µ
in the remainder of the subsection, through a sequence of lemmas.

Theorem 6.5. Let µ
bN , µN be the SIR filters defined by (6.8), (6.7) respectively, with proposal distribution π.
Suppose that there exists fk+1 : X × X → R with
fk+1 (uk+1 , uk ) ∝

P (yk+1 |uk+1 )P (uk+1 |uk )
π(uk+1 |uk , yk+1 )

(6.9)

and satisfying
κ ≤ fk+1 (uk+1 , uk ) ≤ κ−1

19

(6.10)

for all uk+1 , uk ∈ X , k ∈ {0, . . . , K − 1} and some κ ∈ (0, 1]. Then we have
d(b
µN
K , µK ) ≤

K
X

(2κ−2 )k N −1/2

(6.11)

K
X

(2κ−2 )k N −1/2

(6.12)

k=1

and
d(µN
K , µK ) ≤

k=0

for all K, N ≥ 1.
Remark 6.6. As for the boostrap particle filter, the appearance of κ−2 reflects the ratio of the upper and lower
bounds in (6.9); hence there is nothing to be gained from optimizing over the constant of proportionality. If we
let
P (yk+1 |uk+1 )P (uk+1 |uk )
fk+1 (uk+1 , uk ) =
π(uk+1 |uk , yk+1 )
then the estimate (6.10) is equivalent to
κπ(uk+1 |uk , yk+1 ) ≤ P (yk+1 |uk+1 )P (uk+1 |uk ) ≤ κ−1 π(uk+1 |uk , yk+1 ) .
This can thus be interpreted as a quantification of equivalence between measures π and the optimal proposal
P (yk+1 |uk+1 )P (uk+1 |uk ).
Remark 6.7. It is important to note that Assumption 2.1 on the dynamics-observation model is not required
by Theorem 6.5. For the consistency results, Assumption 2.1 is only used to ensure that (6.10) holds. This
observation leads to the following corollary.
Corollary 6.8. Let µ
bN , µN be the OPFs defined in (3.11), (3.7) respectively and satisfying Assumption 2.1.
Then there is κ = κ(YK ) such that we have
d(b
µN
K , µK ) ≤

K
X

(2κ−2 )k N −1/2

(6.13)

K
X

(2κ−2 )k N −1/2

(6.14)

k=1

and
d(µN
K , µK ) ≤

k=0


for all K, N ≥ 1 and where κ−1 = exp max0≤j≤K−1 |yj+1 |2 + supv |Hψ(v)|2S .

Although similar to the argument for the BPF, the recursion argument for the SIR is necessarily more
(n)
(n)
(n)
complicated than that for the BPF, as the weights wk+1 can potentially depend on both uk+1 and uk . This
suggests that we must build a recursion which updates measures on a joint space (uk+1 , uk ) ∈ X × X . This
(n)
would also be necessary if we restricted our attention to the OPF, as the weights are defined using uk and not
(n)
the particle positions u
bk+1 after the proposal.
The recursion is defined using the following three operators.
π
1. First Pk+1
maps probability measures on X to probability measures on X × X by
ZZ
π
Pk+1 µ(A) =
π(uk+1 |uk , yk+1 )µ(duk )duk+1
A

where A is a measurable subset of X × X .

20

2. The reweighting operator Lπk+1 maps probability measures on X × X to probability measures on X × X
and is defined by
ZZ
Lπk+1 Q(A) = Z −1

wk+1 (uk+1 , uk )Q(duk+1 , duk )

A

where Z is the normalization constant of the resulting measure. The weight function is given by
wk+1 (uk+1 , uk ) =

P (yk+1 |uk+1 )P (uk+1 |uk )
.
π(uk+1 |uk , yk+1 )

3. Finally, M maps probability measures on X × X into probability measures on X via marginalization onto
the first component:
Z Z
M Q(B) =

Q(duk+1 , duk )

B×X

.

It is easy to see that the posterior µk satisfies a natural recursion in terms of these operators.
π
Lemma 6.9. µk+1 = M Lπk+1 Pk+1
µk
π
Proof. Let P (uk |Yk ) denote the density of µk , then Pk+1
µk is a measure on X × X with density

π(uk+1 |uk , yk+1 )P (uk |Yk ) .
π
And Lπk+1 Pk+1
µk is a measure on X × X with density

Z −1 wk+1 (uk+1 , uk )π(uk+1 |uk , yk+1 )P (uk |Yk ) = Z −1 P (yk+1 |uk+1 )P (uk+1 |uk )P (uk |Yk ) .
π
Finally, M Lπk+1 Pk+1
µk is a measure on X with density
Z
Z −1 P (yk+1 |uk+1 )P (uk+1 |uk )P (uk |Yk )duk = Z −1 P (yk+1 |uk+1 )P (uk+1 |Yk ) .

(6.15)

X

Similarly, for the normalization factor, we have
Z Z
Z
Z=
P (yk+1 |uk+1 )P (uk+1 |uk )P (uk |Yk )duk duk+1 =
P (yk+1 |uk+1 )P (uk+1 |Yk )duk
X ×X

X

and thus by Bayes’ formula (6.15) is equal to P (uk+1 |Yk+1 ) as required.
We now show that an associated recursion is satisfied by the SIR filter µ
bN .

Lemma 6.10. Let µ
bN be the SIR filter given by (6.8), then

π
N π
µ
bN
bN
k+1 = M Lk+1 S Pk+1 µ
k

for all k ≥ 0 and N ≥ 1, where S N denotes the sampling operator acting on M(X × X ).
PN
(n)
π
Proof. By definition, µ
bN
bN
(n) so that P
k =
k+1 µ
k ∈ M(X × X ) with density
n=1 wk δu
b
k

N
X

n=1

(n)

(n)

(n)

wk π(uk+1 |b
uk , yk+1 )δ(uk − u
bk ) .
(n)

(n)

(n)

π
uk+1 , uk ) obtained as follows: first draw a sample uk
µ
bN
Note that a sample U ∼ Pk+1
k is a pair (b

(n)
(n) N
{b
uk }N
n=1 with weights {wk }n=1 and then draw sample
(n)
π
µ
bN
of the u
bk+1 sequence we see that S N Pk+1
k has density

(n)
u
bk+1

from

(n)
π(uk+1 |uk , yk+1 ).

N
1 X
(n)
(n)
δ(uk+1 − u
bk+1 )δ(uk − uk ) .
N n=1

21

from

Thus, by definition

π
It follows that Lπk+1 S N Pk+1
µN
k has density
N
X

n=1

(n)

(n)

(n)

(n)

Z −1 wk+1 (b
uk+1 , uk )δ(uk+1 − u
bk+1 )δ(uk − uk )

π
and M Lπk+1 S N Pk+1
µN
k has density

N
X

n=1

(n)

(n)

(n)

Z −1 wk+1 (b
uk+1 , uk )δ(uk+1 − u
bk+1 ) .

Lastly, the normalization factor is given by
Z=

Z Z

N
X

X ×X n=1

(n)

(n)

(n)

(n)

wk+1 (b
uk+1 , uk )δ(uk+1 − u
bk+1 )δ(uk − uk )duk duk+1 =

(n)

(n)

N
X

(n)

(n)

wk+1 (b
uk+1 , uk ) ,

n=1

(n)

so that Z −1 wk+1 (b
uk+1 , uk ) = wk+1 and we obtain the result.
In the final step before proving Theorem 6.5, we state some simple properties for the operators appearing in
the recursions. Note that these are similar but not (all) immediately implied by the corresponding results for the
BPF, Lemma 6.3.
Lemma 6.11. We have the following simple estimates:
1. d(M ν, M µ) ≤ d(ν, µ).
π
π
2. d(Pk+1
ν, Pk+1
µ) ≤ d(ν, µ)

3. supν∈M(X ×X ) d(S N ν, ν) ≤ N −1/2
Proof. Let f˜(x, y) = f (x) and let g(x, y) denote an arbitrary function. Then
M ν(f ) − M µ(f ) = ν(f˜) − µ(f˜)
The first inequality follows immediately from taking supremum over all |f | ≤ 1, which is necessarily smaller
than the supremum of ν(g) − µ(g) over all |g| ≤ 1.
We also have
π
π
Pk+1
ν(g) − Pk+1
µ(g) = ν(g π ) − µ(g π )
R
where g π (uk ) = g(uk+1 , uk )π(uk+1 |uk , Yk+1 )duk+1 . And since |g π |∞ ≤ 1, the second inequality follows.
The third inequality is proven in [25, Lemma 4.7], simply replacing X with X × X .
We can now proceed with the main result.
Proof of Theorem 6.5. In the context of Lemma 6.4, take Z = X × X and gk+1 = fk+1 , it follows that
Gk+1 ν = Lπk+1 ν.Indeed, for any ϕ : X × X → R and with gk+1 = Z −1 wk+1 we have
Gk+1 ν(ϕ) =

ν(wk+1 ϕ)
ν(gk+1 ϕ)
=
= Lπk+1 ν(ϕ) .
ν(gk+1 )
ν(wk+1 )

By Assumption 6.10, we therefore obtain from Lemma 6.4 that
d(Lπk+1 µ, Lπk+1 ν) ≤ (2κ−2 )d(µ, ν)
for all µ, ν ∈ M(X × X ).

22

Thus, using the recursions given in Lemmas 6.9, 6.10 and the estimates given in Lemma 6.11, we obtain
π
N π
π
π
d(b
µN
bN
k+1 , µk+1 ) = d(M Lk+1 S Pk+1 µ
k , M Lk+1 Pk+1 µk )
π
π
π
≤ d(Lπk+1 S N Pk+1
µ
bN
k , Lk+1 Pk+1 µk )

π
π
≤ 2κ−2 d(S N Pk+1
µ
bN
k , Pk+1 µk )


π
π
π
π
µ
bN
bN
bN
≤ 2κ−2 d(S N Pk+1
k , Pk+1 µ
k ) + d(Pk+1 µ
k , Pk+1 µk )

≤ 2κ−2 N −1/2 + 2κ−2 d(b
µN
k , µk ) .

N
N N
and since µ
bN
bk
0 = µ0 , we obtain (6.13) by induction. Moreover, since µk = S µ

N N
d(µN
bk , µk ) ≤ d(S N µ
bN
bN
µN
k , µk ) = d(S µ
k ,µ
k ) + d(b
k , µk )

and (6.14) follows.

The corollary follows immediately.
Proof of Corollary 6.8. For the OPF we have
π(uk+1 |uk , yk+1 ) = P (uk+1 |uk , yk+1 ) =

P (yk+1 |uk+1 )P (uk+1 |uk )
P (yk+1 |uk )

where we have applied Bayes formula in the final equality. But under Assumption 2.1 we have that


1
−1
2
P (yk+1 |uk ) = ZS exp − |yk+1 − Hψ(uk )|S .
2
Thus we define fk+1 by


P (yk+1 |uk+1 )P (uk+1 |uk )
1
= exp − |yk+1 − Hψ(uk )|2S
π(uk+1 |uk , yk+1 )
2

and hence (6.10) holds with κ−1 = exp max0≤j≤K−1 |yj+1 |2 + supv |Hψ(v)|2S , which, for each Yk , is finite
by Assumption 2.1. The result follows from Theorem 6.5.
fk+1 (uk+1 , uk ) = ZS

6.3 Gaussianized Optimal Particle Filter
In this section we derive the consistency result for the GOPF.
Theorem 6.12. Let ν N be the GOPF defined by (3.12) and let Assumption 2.1 hold. Then there is κ = κ(YK )
such that
K
X
N
d(νK
, µK ) ≤
(2κ−2 )k N −1/2
k=0

for all K, N ≥ 1, where κ

−1


= exp max0≤j≤K−1 |yj+1 |2 + supv |Hψ(v)|2S .

For the GOPF, the consistency proof uses the same strategy, but turns out to be much more straightforward.
First note that the decomposition of the filtering distribution given in (3.10) gives the recursion formula
µk+1 = Qk+1 Kk+1 µk

(6.16)

where Kk+1 : M(X ) → M(X ) is defined by
Kk+1 µ(A) = Z

−1

Z

P (yk+1 |uk )µ(duk )

A

for all measurable A ⊂ X where Z is the normalization constant, and Qk+1 : M(X ) → M(X ) is the Markov
semigroup with kernel P (uk+1 |uk , yk+1 ).
Moreover, we have the following recursion for the GOPF.
23

Lemma 6.13. The GOPF νkN satsfies the recursion
N
νk+1
= S N Qk+1 Kk+1 νkN

(6.17)

with ν0N = S N µ0 .
Proof. Let νkN =

1
N

PN

(n) ,
n=1 δvk

then Kk+1 ∈ M(X ) with density
N
X

(n)

(n)

Z −1 P (yk+1 |vk )δ(vk − vk ) .

n=1

The normalization constant is given by
Z=

Z X
N

X n=1

(n)

(n)

(n)

P (yk+1 |vk )δ(vk − vk )dvk =

N
X

(n)

P (yk+1 |vk )

n=1

(n)

and thus Z −1 P (yk+1 |vk ) = wk+1 . We then have Qk+1 Lk+1 νkN ∈ M(X ) with density
N
X

(n)

(n)

wk+1 P (vk+1 |vk , yk+1 ) .

n=1

To draw a sample

(n)
vk+1

(n)

(n)

from this mixture model, we draw ṽk

(m)

(m)

N
from {vk }N
m=1 with weights {wk+1 }m=1 and

(n)

N
then draw vk+1 from P (vk+1 |ṽk , yk+1 ). It follows that S N Qk+1 Kk+1 νkN = νk+1
.

Proof. If we let

then gk+1



1
2
gk+1 (vk ) := ZS P (yk+1 |vk ) = exp − |yk+1 − Hψ(vk )|S ,
2
satisfies the assumptions of Lemma 6.4 with


κ−1 = exp
max |yj+1 |2 + sup |Hψ(v)|2S .
0≤j≤K−1

v

In particular, since Gk+1 ν = Kk+1 ν, it follows from Lemma 6.4 that
d(Kk+1 µ, Kk+1 ν) ≤ (2κ−2 )d(µ, ν)
for all µ, ν ∈ M(X ).
Using the recursions (6.16),(6.17) and the estimates from Lemma 6.3, we see that
N
, µk+1 ) = d(S N Qk+1 Kk+1 νkN , Qk+1 Kk+1 µk )
d(νk+1

≤ d(S N Qk+1 Kk+1 νkN , Qk+1 Kk+1 νkN ) + d(Qk+1 Kk+1 νkN , Qk+1 Kk+1 µk )
≤ N −1/2 + d(Kk+1 νkN , Kk+1 µk )
≤ N −1/2 + 2κ−2 d(νkN , µN
k )
by induction, we obtain
N
d(νk+1
, µN
k+1 ) ≤

k
X

(2κ−2 )j N −1/2 + (2κ−2 )k+1 d(ν0N , µ0 ) .

j=0

And the result follows from the fact d(ν0N , µ0 ) = d(S N µ0 , µ0 ) ≤ N −1/2 .
Acknowledgements. DK is supported as a Courant instructor. The work of AMS is supported by EPSRC,
DARPA and ONR.
24

References
[1] Hajime Akashi and Hiromitsu Kumamoto. Random sampling approach to state estimation in switching
environments. Automatica, 13(4):429–434, 1977.
[2] Rami Atar and Ofer Zeitouni. Exponential stability for nonlinear filtering. In Annales de l’IHP Probabilités
et statistiques, volume 33, pages 697–725, 1997.
[3] Peter Bauer, Alan Thorpe, and Gilbert Brunet. The quiet revolution of numerical weather prediction.
Nature, 525(7567):47–55, 2015.
[4] Amarjit Budhiraja and D Ocone. Exponential stability of discrete-time filters for bounded observation
noise. Systems & Control Letters, 30(4):185–193, 1997.
[5] Gerrit Burgers, Peter Jan van Leeuwen, and Geir Evensen. Analysis scheme in the ensemble Kalman filter.
Monthly weather review, 126(6):1719–1724, 1998.
[6] Alexandre J Chorin and Xuemin Tu. Implicit sampling for particle filters. Proceedings of the National
Academy of Sciences, 106(41):17249–17254, 2009.
[7] Dan Crisan and Kari Heine. Stability of the discrete time filter in terms of the tails of noise distributions.
Journal of the London Mathematical Society, 78(2):441–458, 2008.
[8] Pierre Del Moral, Aline Kurtzmann, and Julian Tugaut. On the stability and the uniform propagation of
chaos of extended ensemble Kalman-Bucy filters. arXiv preprint arXiv:1606.08256, 2016.
[9] Randal Douc, Gersende Fort, Eric Moulines, and Pierre Priouret. Forgetting the initial distribution for
hidden markov models. Stochastic processes and their applications, 119(4):1235–1256, 2009.
[10] Arnaud Doucet, Nando De Freitas, and Neil Gordon. An introduction to sequential monte carlo methods.
In Sequential Monte Carlo methods in practice, pages 3–14. Springer, 2001.
[11] Arnaud Doucet, Simon Godsill, and Christophe Andrieu. On sequential monte carlo sampling methods
for Bayesian filtering. Statistics and computing, 10(3):197–208, 2000.
[12] Geir Evensen. Using the extended Kalman filter with a multilayer quasi-geostrophic ocean model. Journal
of Geophysical Research: Oceans, 97(C11):17905–17924, 1992.
[13] Geir Evensen. The ensemble Kalman filter: Theoretical formulation and practical implementation. Ocean
dynamics, 53(4):343–367, 2003.
[14] Michael Ghil and Paola Malanotte-Rizzoli. Data assimilation in meteorology and oceanography. Advances
in geophysics, 33:141–266, 1991.
[15] John Harlim, Andrew J Majda, et al. Catastrophic filter divergence in filtering nonlinear dissipative systems. Communications in Mathematical Sciences, 8(1):27–43, 2010.
[16] Peter L Houtekamer and Herschel L Mitchell. Ensemble kalman filtering. Quarterly Journal of the Royal
Meteorological Society, 131(613):3269–3289, 2005.
[17] David Kelly, Andrew J Majda, and Xin T Tong. Concrete ensemble Kalman filters with rigorous catastrophic filter divergence. Proceedings of the National Academy of Sciences, 112(34):10589–10594, 2015.
[18] David Kelly, Eric Vanden-Eijnden, and Jonathan Weare. Implicit resampling in the optimal particle filter.
2016.
[19] DTB Kelly, KJH Law, and Andrew M Stuart. Well-posedness and accuracy of the ensemble Kalman filter
in discrete and continuous time. Nonlinearity, 27(10):2579, 2014.

25

[20] ML Kleptsyna and A Yu Veretennikov. On discrete time ergodic filters with wrong initial data. Probability
Theory and Related Fields, 141(3-4):411–444, 2008.
[21] Hiroshi Kunita. Asymptotic behavior of the nonlinear filtering errors of markov processes. Journal of
Multivariate Analysis, 1(4):365–393, 1971.
[22] K Law, AM Stuart, KC Zygalakis, et al. Accuracy and stability of the continuous-time 3dvar filter for the
Navier-Stokes equation. Nonlinearity, 26(8):2193, 2013.
[23] KJH Law, D Sanz-Alonso, A Shukla, and AM Stuart. Filter accuracy for the lorenz 96 model: fixed versus
adaptive observation operators. Physica D: Nonlinear Phenomena, 325:1–13, 2016.
[24] KJH Law, Abhishek Shukla, and AM Stuart. Analysis of the 3dvar filter for the partially observed lorenz
’63 model. Discrete and Continuous Dynamical Systems A, 34:1061–1078, 2014.
[25] Kody Law, Andrew Stuart, and Kostas Zygalakis. Data Assimilation. Springer, 2015.
[26] Kody JH Law, Hamidou Tembine, and Raul Tempone. Deterministic mean-field ensemble kalman filtering. arXiv preprint arXiv:1409.0628, 2014.
[27] Jun S Liu and Rong Chen. Blind deconvolution via sequential imputations. Journal of the american
statistical association, 90(430):567–576, 1995.
[28] Andrew C Lorenc. Analysis methods for numerical weather prediction. Quarterly Journal of the Royal
Meteorological Society, 112(474):1177–1194, 1986.
[29] Andrew J Majda and John Harlim. Filtering complex turbulent systems. Cambridge University Press,
2012.
[30] Andrew J Majda and Xin T Tong. Robustness and accuracy of finite ensemble kalman filters in large
dimensions. arXiv, arXiv:1606.09321, 2016.
[31] Jonathan C Mattingly, Andrew M Stuart, and Desmond J Higham. Ergodicity for sdes and approximations: locally lipschitz vector fields and degenerate noise. Stochastic processes and their applications,
101(2):185–232, 2002.
[32] Sean P Meyn and Richard L Tweedie. Markov Chains and Stochastic Stability. Springer Science &
Business Media, 2012.
[33] Alexander JF Moodey, Amos S Lawless, Roland WE Potthast, and Peter Jan Van Leeuwen. Nonlinear
error dynamics for cycled data assimilation methods. Inverse Problems, 29(2):025002, 2013.
[34] Patrick Rebeschini, Ramon Van Handel, et al. Can local particle filters beat the curse of dimensionality?
The Annals of Applied Probability, 25(5):2809–2866, 2015.
[35] Daniel Sanz-Alonso and Andrew M Stuart. Long-time asymptotics of the filtering distribution for partially
observed chaotic dynamical systems. SIAM/ASA Journal on Uncertainty Quantification, 3(1):1200–1220,
2015.
[36] Daniel Sanz-Alonso and Andrew M Stuart. Long-time asymptotics of the filtering distribution for partially
observed chaotic dynamical systems. SIAM/ASA Journal on Uncertainty Quantification, 3(1):1200–1220,
2015.
[37] Madineh Sedigh-Sarvestani, David J Albers, and Bruce J Gluckman. Data assimilation of glucose dynamics for use in the intensive care unit. In 2012 Annual International Conference of the IEEE Engineering in
Medicine and Biology Society, pages 5437–5440. IEEE, 2012.
[38] Madineh Sedigh-Sarvestani, Steven J Schiff, and Bruce J Gluckman. Reconstructing mammalian sleep
dynamics with data assimilation. PLoS Comput Biol, 8(11):e1002788, 2012.
26

[39] Chris Snyder. Particle filters, the “optimal” proposal and high-dimensional systems. In Proceedings of the
ECMWF Seminar on Data Assimilation for atmosphere and ocean, 2011.
[40] Chris Snyder, Thomas Bengtsson, and Mathias Morzfeld. Performance bounds for particle filters using
the optimal proposal. Monthly Weather Review, 143(11):4750–4761, 2015.
[41] Sebastian Thrun. Particle filters in robotics. In Proceedings of the Eighteenth conference on Uncertainty
in artificial intelligence, pages 511–518. Morgan Kaufmann Publishers Inc., 2002.
[42] Xin T Tong, Andrew J Majda, and David Kelly. Nonlinear stability and ergodicity of ensemble based
Kalman filters. Nonlinearity, 29(2):657, 2016.
[43] Xin T Tong, Andrew J Majda, and David Kelly. Nonlinear stability of the ensemble Kalman filter with
adaptive covariance inflation. Communications in Mathematical Sciences, 14(5), 2016.
[44] Xin Thomson Tong and Ramon Van Handel. Ergodicity and stability of the conditional distributions of
nondegenerate Markov chains. The Annals of Applied Probability, pages 1495–1540, 2012.
[45] Ramon Van Handel. The stability of conditional Markov processes and Markov chains in random environments. The Annals of Probability, pages 1876–1925, 2009.
[46] Peter Jan Van Leeuwen. Particle filtering in geophysical systems. Monthly Weather Review, 137(12):4089–
4114, 2009.
[47] Peter Jan van Leeuwen. Nonlinear data assimilation in geosciences: an extremely efficient particle filter.
Quarterly Journal of the Royal Meteorological Society, 136(653):1991–1999, 2010.
[48] Daniel B Work, Sébastien Blandin, Olli-Pekka Tossavainen, Benedetto Piccoli, and Alexandre M Bayen.
A traffic model for velocity data assimilation. Applied Mathematics Research eXpress, 2010(1):1–35,
2010.
[49] VS Zaritskii, VB Svetnik, and LI Šimelevič. Monte-carlo technique in problems of optimal information
processing. Avtomatika i Telemekhanika, (12):95–103, 1975.

27

