remarkRemark exampleExample Ensemble Kalman Methods: A Mean Field
PerspectiveEdoardo Calvello, Sebastian Reich, Andrew M. Stuart

# Ensemble Kalman Methods: A Mean Field Perspective

Edoardo Calvello California Institute of Technology, Pasadena, CA 91125, USA
<EMAIL_ADDRESS>Sebastian Reich Institut für Mathematik,
Universität Potsdam, D-14476 Potsdam, Germany (). sebastian.reich@uni-
potsdam.de Andrew M. Stuart California Institute of Technology, Pasadena, CA
91125, USA<EMAIL_ADDRESS>
###### Abstract

This paper provides a unifying mean field based framework for the derivation
and analysis of ensemble Kalman methods. Both state estimation and parameter
estimation problems are considered, and formulations in both discrete and
continuous time are employed. For state estimation problems both the control
and filtering approaches are studied; analogously, for parameter estimation
(inverse) problems the optimization and Bayesian perspectives are both
studied. The approach taken unifies a wide-ranging literature in the field,
provides a framework for analysis of ensemble Kalman methods, and suggests
open problems.

##  1 Introduction

The ensemble Kalman methodology is an innovative and flexible set of tools
which can be used for both state estimation in dynamical systems and parameter
estimation for generic inverse problems. It has primarily been developed by
practitioners in the geophysical sciences, with notable impact on the fields
of oceanography, oil reservoir simulation and weather forecasting. Despite its
wide adoption in the geosciences, the methodology is hard to analyze and firm
theoretical foundations are only now starting to emerge. The purpose of this
article is twofold: a) to describe a mathematical framework for the analysis
of ensemble Kalman methods, describing what is known and highlighting the many
open mathematical challenges in the field; b) to provide a literature survey
which bridges the domain-specific development of the methodology with emerging
mathematical analyses. In so doing we will also highlight the flexibility of
the methodology for use in widespread applications, beyond its historical
development in the geosciences.

The material is organized around the two separate themes of state estimation
and inverse problems; within each, both discrete time and continuous time
approaches are explained. The novel perspective which underlies all of this
material is the derivation of ensemble Kalman methods as particle
approximations of carefully designed mean-field models. The relationship of
these mean-field models to exact transport models, for Gaussian problems,
serves to motivate their form.

In Subsection 1.1 we overview the history of ensemble Kalman methods.
Subsection 1.2 describes the organization of the paper. In Subsection 1.3 we
make brief remarks about the pseudo-code that we make available as a
supplementary resource for this paper. The introduction concludes, in
Subsection 1.4, with a summary of the notation that we adopt throughout the
paper.

###  1.1 Historical Context

The Kalman filter (KF) is arguably the first setting in which the integration
of observational data with a dynamical system was considered, leading to both
discrete time (Kalman, 1960) and continuous time (Kalman & Bucy, 1961)
formulations; see Welch et al. (1995) for an overview. The Kalman filter
applies only in the setting of linear Gaussian dynamics and observation; it
computes the distribution of the state of the dynamical system, given
observations, exactly. The extended Kalman filter (historically denoted EKF;
however the acronym ExKF is also used and is useful) was introduced in order
to extend Kalman’s ideas to nonlinear problems; see the texts Jazwinski (2007)
and Anderson & Moore (2012) for overviews. The extended Kalman approach is
based on a linearization approximation; it hence fails to exactly compute the
distribution of the state of the dynamical system, given observations, in
general. Furthermore it requires propagation of covariance matrices which can
be very large for applications arising in the geosciences (Ghil et al., 1981).

The ensemble Kalman filter (EnKF) was introduced in the celebrated paper
Evensen (1994) which made the consequential observation that, if an _ensemble_
of state estimators is employed then it can also be used to make an
approximation of the covariance. In geosciences applications this circumvents
the computation of large covariances, replacing them instead with low rank
approximations, with rank determined by the number of ensemble members. The
original paper developed the idea in the context of ocean models, but was
rapidly and concurrently developed in a variety of geoscience application
domains (Van Leeuwen & Evensen, 1996; Burgers et al., 1998; Houtekamer &
Mitchell, 1998); the paper van Leeuwen (2020) provide a historical overview.
These methods are sometimes referred as the _stochastic EnKF_ : they require
simulation of random variables to implement. A different class of ensemble
methods, known collectively as _ensemble square root filters_ , was
subsequently developed (Anderson, 2001; Whitaker & Hamill, 2002; Bishop et
al., 2001; Hunt et al., 2007; Tippett et al., 2003; Sakov & Oke, 2008); these
methods are a form of _deterministic EnKF_ : they do not require simulation of
random variables to implement.

Central to our mathematical presentation of Kalman-based methods is the
adoption of mean-field and transport perspectives on the subject. The
incorporation of data within filtering constitutes (possibly approximate)
application of Bayes’ theorem; the papers Daum et al. (2010), Reich (2011), El
Moselhy & Marzouk (2012), and the book Cotter & Reich (2013) introduce novel
approaches to Bayesian inversion, rooted in transport and mean-field models.
While El Moselhy & Marzouk (2012), Spantini et al. (2019a) propose a direct
numerical approximation of the underlying optimal transportation problem, Daum
et al. (2010) and Reich (2011) pursue a homotopy approach. We note that the
homotopy approach is closely related to iterative implementations of the EnKF,
as first considered by Li & Reynolds (2009), Gu & Oliver (2007), and Sakov et
al. (2012); these homotopy approaches lead to continuous time formulations of
the EnKF in the limit of infinitely many iterations, as first considered by
Bergemann & Reich (2010a) and Bergemann & Reich (2010b).

The key connection between mean-field models and ensemble methods is that the
latter can be derived as particle approximations of the mean-field limit; this
viewpoint will play a guiding role in our presentation of the subject of
ensemble methods in these notes. In this context it is notable that the field
of optimization, which is linked to Bayesian sampling through MAP estimation
(Kaipio & Somersalo, 2006), has also seen recent development using mean-field
models – see Carrillo et al. (2018) for an overview and unifying mathematical
framework.

The methods covered in this survey provide only approximate solutions to the
underlying filtering, inference and/or optimization problem. The
approximations invoked are based on assuming linear Gaussian structure, but
applying the resulting methodology outside this regime. On the other hand,
alternative methods, such as sequential Monte Carlo (resp. Monte Carlo Markov
chain) can be designed to be consistent with the underlying nonlinear
filtering problem (resp. Bayesian inverse problem). The books Doucet et al.
(2001) and Chopin & Papaspiliopoulos (2020) overview use of the sequential
Monte Carlo methods for general discrete time filtering and inference problems
(resp. (Kaipio & Somersalo, 2006) for Bayesian inverse problems). However,
sequential Monte Carlo methods suffer from the curse of dimensionality and are
presently not directly applicable to high dimensional problems as arising, for
example, from geophysical applications. This issue with the curse of
dimensionality provides the main motivation for the type of methods covered in
this survey. See Snyder et al. (2008), Bickel et al. (2008), Rebeschini & Van
Handel (2015), and Agapiou et al. (2017) for details.

This concludes our chronological overview of the historical context for the
development of ensemble Kalman methods, and the specific mathematical context
which will be our focus. Each section of the paper concludes with a
bibliographic subsection in which a deeper literature review is given.

###  1.2 Overview

Section 2 is devoted to the problem of state estimation for discrete time
(possibly stochastic) dynamical systems subject to noisy observations. We
formulate the problem from the perspectives of both control theory and
probability and we provide a unifying approach to algorithms for these
problems; the approach rests on transport of measures and mean field
stochastic dynamical systems. Ensemble Kalman methods are then derived as
particle approximations of the mean field models. Section 3 adopts a
perspective that parallels the previous section, but works in a continuous
time setting. Stochastic differential equations (SDEs) are used to describe
the state and its observation process, and mean field SDEs, and related
stochastic partial differential equations (SPDEs), are used to provide the
underpinnings of algorithms; particle approximations of the mean field systems
give interacting systems of SDEs which describe ensemble Kalman methods. The
formulation in continuous time is useful both because in some applications
state estimation problems are most naturally formulated this way, and because
they provide insight into discrete time algorithms, giving rise to cleaner
analysis of phenomena present in both discrete and continuous time. Section 4
is devoted to the use of ensemble Kalman methods for parameter estimation and
inverse problems, demonstrating how a useful change of perspective opens up
their use in this broader setting; both discrete and continuous time are
considered in parallel.

Figure 1: Organizational flow of sections 2, 3 and 4.

Sections 2, 3 and 4 are all organized according to the flow of ideas displayed
in Figure 1. Indeed, starting from a probabilistic perspective, that describes
the evolution of the filtering distribution, it is possible to formulate mean
field maps whose sample paths on state space are equal in law to the filtering
evolution. Since it is typically not tractable to identify these maps
explicitly, it is of interest to determine mean field maps whose sample paths
are only approximately equal in law to the filtering evolution. This leads to
the idea of second order mean field models whose sample paths have law with
first and second order moments which match those of the filtering
distribution, conditional on the prediction step in the filtering cycle.
Particle approximations of these second order mean field maps are then used to
derive implementable numerical algorithms which form the class of ensemble
Kalman methods.

Section 5 concludes, and in particular highlights open problems in the area
which might be of particular interest to the mathematical community. The
appendices in Sections 6 and 7 are devoted to pseudo-code for the algorithms
introduced and to background on the underlying example that we use throughout
the paper. The appendices in Sections 8 and 9 provide underpinning material on
mean-field maps and on stochastic integration.

###  1.3 Pseudo-Code

Pseudo-code describing several of the algorithms that we present and deploy in
this paper are described in Appendix 6. The reader is encouraged to consult
Algorithms 1 and 2, 3DVAR and the ensemble Kalman filter (EnKF) respectively,
in the context of the problem of state estimation for discrete time dynamical
systems presented in Section 2. The scheme 3DVAR is employed in Examples 2.1,
2.2, 2.3 and 2.13. The ensemble Kalman filter is applied in the context of
Example 2.13. Ensemble Kalman methods for inversion, as shown in Algorithms 3,
4 and 5, are presented in Section 4 and applied in the context of Example
4.33.

###  1.4 Notation

Throughout we denote the positive integers and non-negative integers
respectively by $\mathbb{N}=\\{1,2,\cdots\\}$ and
$\mathbb{Z}^{+}=\\{0,1,2,\cdots\\},$ and the notation
$\mathbb{R}=(-\infty,\infty)$ and $\mathbb{R}^{+}=[0,\infty)$ for the reals
and the non-negative reals. We let $\langle\cdot,\cdot\rangle,|\cdot|$ denote
the Euclidean inner-product and norm, noting that $|v|^{2}=\langle
v,v\rangle.$ We also use $|\cdot|$ to denote the resulting induced norm on
matrices. We use $:$ to denote the Frobenius inner-product between matrices,
and $|\cdot|_{\rm F}$ the induced norm on matrices. For any function
$g:\mathbb{R}^{d_{1}}\mapsto\mathbb{R}^{d_{2}}$ we denote by $Dg(\cdot)$ the
derivative.

We use $\mathbb{E}$ and $\mathbb{P}$ to denote expectation and probability
under the prevailing probability measure; if we wish to make clear that
measure $\pi$ is the prevailing probability measure then we write
$\mathbb{E}^{\pi}$ and $\mathbb{P}^{\pi}$. We use
$\textit{Law}\,(\mathsf{rv})$ to denote the law of random variable
$\mathsf{rv}$. We let $\delta_{u}$ denote the Dirac mass on $\mathbb{R}^{d}$,
centered at point $u\in\mathbb{R}^{d}.$ We let $\mathsf{N}(m,C)$ denote the
distribution of a Gaussian random variable with mean $m$ and covariance $C$.

We write $A\succ B$ when $A-B$ is positive definite and $A\succeq B$ when
$A-B$ is positive semi-definite. We will also write $A\prec B$ when $B-A\succ
0$ and $A\preceq B$ when $B-A\succeq 0.$

For $C\succ 0$ (and therefore a covariance matrix) we define
$\langle\cdot,\cdot\rangle_{C},|\cdot|_{C}$ the covariance-weighted Euclidean
inner-product and norm, by $\langle u,v\rangle_{C}=\left\langle
u,C^{-1}v\right\rangle$ and $|v|_{C}^{2}=\langle v,v\rangle_{C}.$

##  2 Discrete Time

In Subsection 2.1 we provide the set-up for the problem of state estimation in
discrete time. Subsection 2.2 introduces an algorithm for this problem based
on a control theoretic perspective. Subsection 2.3 describes the Bayesian
probabilistic perspective; algorithms are not presented but foundations for
the creation of algorithms are laid through the decomposition of the filtering
cycle into an iteration which alternates prediction and the assimilation of
data. In Subsection 2.4 we introduce the important idea of Gaussian
projection, and the resulting Gaussian projected filter. Subsection 2.5
introduces various mean-field dynamical systems which approximate the
filtering cycle; a unifying transport map perspective is adopted. This leads,
in Subsection 2.6, to definitions of the ensemble Kalman filter, the ensemble
adjustment filter and the ensemble transform filter, all derived as particle
approximations of specific mean-field models. We conclude with bibliographic
notes in Subsection 2.7 in which we review relevant literature, and include
discussion of a variety of other algorithms for state estimation, and relate
them to the perspective we adopt here.

###  2.1 Set-Up

The objective of sequential data assimilation is to iteratively update the
state of a (possibly stochastic) dynamical system based on observations and
knowledge of the dynamical and observational processes; we refer to this as
state estimation. The typical setting is one in which the initial condition is
uncertain, but this uncertainty is compensated for by (typically noisy)
observations of a (possibly nonlinear) function of the state. These
observations often live in a space of lower dimension than the dimension of
the state space meaning that the goal of state estimation goes beyond
denoising, and into the realm of control theoretic considerations such as
observability.

A useful starting point is to consider a stochastic dynamical system in which
the evolution of the state, and the relationship between the observations
(which we also refer to as data) and the state, are defined, respectively, by
the equations

(1a) $\displaystyle v_{n+1}$ $\displaystyle=\Psi(v_{n})+\xi_{n},$ (1b)
$\displaystyle y_{n+1}$ $\displaystyle=h(v_{n+1})+\eta_{n+1},$

taken to hold for all $n\in\mathbb{Z}^{+}.$ We assume that, for each fixed
$n\in\mathbb{Z}^{+}$, the _state_ $v_{n}\in\mathbb{R}^{d_{v}}$ and the
_observations_ $y_{n}\in\mathbb{R}^{d_{y}}.$ The maps $\Psi(\cdot)$ and
$h(\cdot)$ describe the systematic, deterministic components of the dynamics
and observation processes, and are assumed known. The initial condition for
$v_{0}$ is assumed random and the systematic components of the model are
subjected to mean zero noise, $\xi_{n}$ and $\eta_{n+1}$. To be concrete we
assume $v_{0},\\{\xi_{n}\\}_{n\in\mathbb{Z}^{+}}$ and
$\\{\eta_{n}\\}_{n\in\mathbb{N}}$ are mutually independent Gaussians defined
by

(2)
$v_{0}\sim\mathsf{N}(m_{0},C_{0}),\quad\xi_{n}\sim\mathsf{N}(0,\Sigma)\,\,\text{i.i.d.},\quad\eta_{n}\sim\mathsf{N}(0,\Gamma)\,\,\text{i.i.d.}\,\,.$

In practice we will have available to us the observation coordinates of a
specific true realization of the random dynamical system (1), from which we
wish to recover the specific true realization of the state that gave rise to
these observations. We denote the true realizations of the state of the system
by sequence $\\{{v}^{\dagger}_{n}\\}_{n\in\mathbb{Z}^{+}}$, and the observed
data by $\\{{y}^{\dagger}_{n}\\}_{n\in\mathbb{N}}.$ These are generated by
${v}^{\dagger}_{0},\\{{\xi}^{\dagger}_{n}\\}_{n\in\mathbb{Z}^{+}}$ and
$\\{{\eta}^{\dagger}_{n}\\}_{n\in\mathbb{N}}$, specific realizations of the
initial condition and state and observational noise from the distribution
defined by (2).

We let ${Y}^{\dagger}_{n}=\\{{y}^{\dagger}_{\ell}\\}_{1\leq\ell\leq n}.$ Our
objective is to estimate the state ${v}^{\dagger}_{n}$ at time $n$ from data
${Y}^{\dagger}_{n}$. Specifically, it is natural to think about this design
objective in two different ways:

  * •

Objective 1: design an algorithm producing output $v_{n}$ from
${Y}^{\dagger}_{n}$ so that $\\{v_{n}\\}_{n\in\mathbb{Z}^{+}}$ estimates
$\\{{v}^{\dagger}_{n}\\}_{n\in\mathbb{Z}^{+}}$, the true state generated by
(1a);

  * •

Objective 2: design an algorithm which estimates the distribution of random
variable $v_{n}|{Y}^{\dagger}_{n}$.

In both cases we are interested in Markovian formulations which update the
estimate $v_{n}$, or the distribution $v_{n}|{Y}^{\dagger}_{n}$, sequentially
as the data is acquired. In the next two subsections we describe control
theoretic and probabilistic approaches to this problem which, respectively,
provide the basis for algorithms addressing Objectives 1 and 2. We note that
fulfilling Objective 2 immediately implies resolution of Objective 1, for
example by taking the mean of $v_{n}|{Y}^{\dagger}_{n}$; but the reverse is
not typically true. However, Objective 1 is easier to address, is especially
relevant when noise levels are small, and more generally serves to motivate
approaches which are used to address Objective 2.

###  2.2 Control Theory Perspective

A very natural idea from control theory underlies ensemble Kalman filtering
and it is encapsulated in the following algorithmic approach. This way of
attacking state estimation is most appropriate when $|C_{0}|,|\Gamma|$ and
$|\Sigma|$ are small so that the state and observations are close to
deterministic. To understand this setting we will first study algorithms in
which the covariances $\Gamma$ and $\Sigma$ are set to zero, and then return
to the inclusion of noise.

The algorithmic idea works as follows: from current state estimate $v_{n}$,
given ${Y}^{\dagger}_{n}$, predict the outcome of the model and data, denoted
by $(\widehat{v}_{n+1},\widehat{h}_{n+1})$, using the update equations (1),
but ignoring the noise; then correct the state estimate by nudging the
prediction using the mismatch between observed and predicted data
$({y}^{\dagger}_{n+1},\widehat{h}_{n+1})$. This results in a deterministic map
$v_{n}\mapsto v_{n+1}$, assumed to hold for all $n\in\mathbb{Z}^{+}$, of the
following form:

(3a) $\displaystyle\widehat{v}_{n+1}$ $\displaystyle=\Psi(v_{n}),$ (3b)
$\displaystyle\widehat{h}_{n+1}$ $\displaystyle=h(\widehat{v}_{n+1}),$ (3c)
$\displaystyle v_{n+1}$
$\displaystyle=\widehat{v}_{n+1}+K({y}^{\dagger}_{n+1}-\widehat{h}_{n+1}).$

Equations (3a, 3b) create predicted state and data. The difference between the
predicted data and the true data
$\\{{y}^{\dagger}_{n+1}\\}_{n\in\mathbb{Z}^{+}}$, the latter found from a
fixed realization of (1), is then used to correct the predicted state
resulting in (3c). We can write the algorithm compactly in the form

(4)
$v_{n+1}=\Psi(v_{n})+K\Bigl{(}{y}^{\dagger}_{n+1}-h\bigr{(}\Psi(v_{n})\bigr{)}\Bigr{)}.$

Choice of the gain matrix $K$ will complete definition of an algorithm which
we refer to as 3DVAR.111The nomenclature “3DVAR” was coined in the geophysics
community, and stands for three dimensional variational data assimilation.
This is natural for algorithms which incorporate spatially distributed data
sequentially in time, in the context where the state variable $v$ varies
across three spatial coordinates. In our setting the state variable $v$ is not
required to have any spatial structure, but the control formulation (4)
reproduces the 3DVAR algorithm from the geophysics community when it does.
Hence we still refer to it as 3DVAR. Pseudo-code for 3DVAR may be found as
Algorithm 1 in Appendix 6. We note that the difference between the observed
value ${y}^{\dagger}_{n+1}$ and its estimator
$\widehat{h}_{n+1}=h\bigr{(}\Psi(v_{n})\bigr{)}$ is often referred to as the
innovation, and for this we introduce the notation

(5) $I_{n}={y}^{\dagger}_{n+1}-\widehat{h}_{n+1}.$

In the following example we illustrate this algorithm. The example serves to
demonstrate the viability of the proposed control theoretic way of thinking
about state estimation, and addresses Objective 1. In subsequent subsections
we will show how the ideas can be generalized to an adaptive gain matrix
$K_{n}$, leading us to the ensemble Kalman methodology and to addressing both
Objectives 1 and 2.

######  Example 2.1.

To illustrate the problem of state estimation we consider the Lorenz ’96
(singlescale) model for unknown $v\in C(\mathbb{R}^{+},\mathbb{R}^{L})$
satisfying the equations

(6a) $\displaystyle\dot{v}_{\ell}$
$\displaystyle=-v_{\ell-1}(v_{\ell-2}-v_{\ell+1})-v_{\ell}+F+h_{v}m\bigl{(}v_{\ell}\bigr{)},\quad\ell=1\dots
L\,,$ (6b) $\displaystyle v_{\ell+L}$ $\displaystyle=v_{\ell},\quad\ell=1\dots
L\,.$

Here we set $L=9,h_{v}=-0.8$ and $F=10.$ Function $m$ is shown in Figure 2.
222This specific choice of $m$ does not make any material difference to what
is presented in this example; any $m$ for which the equation is well-posed
will lead to similar conclusions. However the specific choice of $m$ shown in
Figure 2 is relevant within Example 2.2. We note that setting $h_{v}=0$ in (6)
leads to the standard singlescale Lorenz’96 model.

Figure 2: Plot of the approximation $m$.

Let $\Psi$ denote the solution operator for (6) over the observation time
interval $\tau$. We emphasize that, at the parameter values we have chosen,
(6) is chaotic and exhibits sensitivity to perturbations. We consider
observations $\\{{y}^{\dagger}_{n}\\}_{n\in\mathbb{Z}^{+}}$ arising from the
model

$\displaystyle{v}^{\dagger}_{n+1}$
$\displaystyle=\Psi({v}^{\dagger}_{n})+{\xi}^{\dagger}_{n},$
$\displaystyle{y}^{\dagger}_{n+1}$
$\displaystyle=h({v}^{\dagger}_{n+1})+{\eta}^{\dagger}_{n+1},$

where $\\{{\xi}^{\dagger}_{n}\\}_{n\in\mathbb{Z}^{+}}$,
$\\{{\eta}^{\dagger}_{n}\\}_{n\in\mathbb{N}}$ are mutually independent
Gaussian sequences defined by

${\xi}^{\dagger}_{n}\sim\mathsf{N}(0,\sigma^{2}I)\,\,\text{i.i.d.},\quad{\eta}^{\dagger}_{n}\sim\mathsf{N}(0,\gamma^{2}I)\,\,\text{i.i.d.}\,\,.$

Because of the chaotic nature of the dynamical system defined by iteration of
$\Psi$, a key question concerning the problem of determining the state
${v}^{\dagger}_{n}$ from ${Y}^{\dagger}_{n}$ is whether the observations
compensate for the sensitive dependence of the state evolution, enabling
accurate recovery of the state.

We assume that the observation function is linear: $h(v)=Hv$ for matrix
$H:\mathbb{R}^{9}\to\mathbb{R}^{6}$ defined by

(7) $Hv=(v_{1},v_{2},v_{4},v_{5},v_{7},v_{8})^{\top}.$

The 3DVAR algorithm (3) reduces, in the setting of this example, to the
mapping

(8) $v_{n+1}=(I-KH)\Psi(v_{n})+K{y}^{\dagger}_{n+1}.$

We define the filter by choosing gain $K:\mathbb{R}^{6}\to\mathbb{R}^{9}$ to
be

(9) $Kw=(w_{1},w_{2},0,w_{3},w_{4},0,w_{5},w_{6},0)^{\top}.$

To interpret the algorithm, and motivate the choice of $K$, given $H$, notice
that

(10a) $\displaystyle KHv$
$\displaystyle=(v_{1},v_{2},0,v_{4},v_{5},0,v_{7},v_{8},0)^{\top},$ (10b)
$\displaystyle(I-KH)v$ $\displaystyle=(0,0,v_{3},0,0,v_{6},0,0,v_{9})^{\top},$
(10c) $\displaystyle HK$ $\displaystyle=I.$

Note that, applying the observation map $H$ to the recursion (8) and using
(10c) we find that

$Hv_{n+1}={y}^{\dagger}_{n+1}=H\bigl{(}\Psi({v}^{\dagger}_{n})+{\xi}^{\dagger}_{n}\bigr{)}+{\eta}^{\dagger}_{n+1}.$

Assuming that $\sigma^{2}$ and $\gamma^{2}$ are small and neglecting the noise
contributions shows that

${y}^{\dagger}_{n+1}\approx H\Psi({v}^{\dagger}_{n})\approx
H{v}^{\dagger}_{n+1}.$

Thus, ignoring small noise perturbations, ${y}^{\dagger}_{n+1}\approx
H{v}^{\dagger}_{n+1}.$ Then, using (10a) and (10b) we see that the algorithm
(8) has the very natural approximate interpretation of iterating using the
model and then overwriting using the observed true state in the observed
components. This explains why the specific choice of $K$ is reasonable.

Because of this interpretation it is natural to study the 3DVAR algorithm, for
this example, by displaying the output of 3DVAR on one of the unobserved
components, and comparing with the truth; the key question is whether the
observed components induce synchronization of 3DVAR with the truth in the
unobserved components. Thus, in the following numerical experiments, we
display component $v_{3}.$

(a)

(b)

Figure 3: In both (a) and (b) the estimates of $v_{3}$ in time produced by
3DVAR using $\tau=10^{-3}$ are displayed. In (a) $\sigma^{2}$ and $\gamma^{2}$
are set to $10^{-3}$, while in (b) to $10^{-1}$. The acronym “dns” refers to
direct numerical simulation and, in both cases, it is noteworthy that 3DVAR is
able to synchronize with the true trajectory of the chaotic dynamical system,
even though it is initialized differently from the direct numerical
simulation.

(a)

(b)

Figure 4: In both (a) and (b) the noise variances $\sigma^{2}$ and
$\gamma^{2}$ are set to $10^{-3}$. In (a) we display the estimates of $v_{3}$
in time produced by 3DVAR against the true dynamics using $\tau=5\cdot
10^{-1}$, while in (b) using $\tau=10^{0}$. Again the acronym “dns” refers to
direct numerical simulation. 3DVAR successfully synchronizes with the direct
numerical simulation at the smaller value of $\tau$ but fails to do so when
$\tau$ is larger.

Figure 3 illustrates the foregoing intuition about the behavior of 3DVAR in
experiments conducted with the choice $\tau=10^{-3}$. For small noise with
variances of size $10^{-3}$ we observe in Figure 3(a) the phenomenon of near
perfect synchronization of the 3DVAR output with the truth. Although not shown
here, this synchronization occurs in all components of the solution, observed
and unobserved, and thus the entire state of 3DVAR synchronizes with the
truth, up to a small error on the scale of the noise. The algorithm thus
produces an accurate estimate of the true state. In the scenario of larger
state and observational noise, of variance $10^{-1}$, displayed in Figure
3(b), 3DVAR still captures the correct trend of the true dynamics, but there
are clear overshoots and undershoots in the estimates because the noise is
larger and not accounted for in the 3DVAR algorithm.

It is intuitive that the synchronization phenomena studied above will depend
not only on the size of the noise, but also on the inter-observation time
$\tau.$ Figure 4 illustrates the effect of varying $\tau$, in the setting
where the variance is $10^{-3}.$ The simulations show that for $\tau$ of order
$\mathcal{O}(10^{-1})$, 3DVAR cannot accurately estimate the true state.
$\blacksquare$

A common source of error in data assimilation arises from the fact that the
data is not produced by the mathematical model used for assimilation. This is
known as _model misspecification_. To illustrate the fact that data
assimilation algorithms can be used to address this issue we study it in the
context of the Lorenz ’96 model. Specifically we will generate data using the
Lorenz ’96 multiscale model, and assimilate data using the Lorenz ’96
singlescale model (6). The relationship between the multiscale and singlescale
models is detailed in Appendix 7.

######  Example 2.2.

To define the multiscale model, let $v\in C(\mathbb{R}^{+},\mathbb{R}^{L})$
and $w\in C(\mathbb{R}^{+},\mathbb{R}^{L\times J})$. Each variable
$v_{\ell}\in\mathbb{R}$ is coupled to a subgroup of fast variables
$\\{w_{\ell,j}\\}_{j=1}^{J}\in\mathbb{R}^{J}$. For $\ell=1\dots L$ and
$j=1\dots J$, we write

(11a) $\displaystyle\dot{v}_{\ell}$
$\displaystyle=f_{\ell}(v)+h_{v}\overline{w}_{\ell},\quad\overline{w}_{\ell}=\frac{1}{J}\sum_{j=1}^{J}w_{\ell,j},$
(11b) $\displaystyle\dot{w}_{\ell,j}$
$\displaystyle=\frac{1}{\epsilon}r_{j}(v_{\ell},w_{\ell})$

where

(12a) $\displaystyle f_{\ell}(v)$
$\displaystyle=-v_{\ell-1}(v_{\ell-2}-v_{\ell+1})-v_{\ell}+F$ (12b)
$\displaystyle r_{j}(v_{\ell},w_{\ell})$
$\displaystyle=-w_{\ell,j+1}(w_{\ell,j+2}-w_{\ell,j-1})-w_{\ell,j}+h_{w}v_{\ell},$

and we impose the boundary conditions

(13) $v_{\ell+L}=v_{\ell},\quad w_{\ell+L,j}=w_{\ell,j},\quad
w_{\ell,j+J}=w_{\ell+1,j}.$

Here $\epsilon>0$ is a scale-separation parameter, $h_{v},h_{w}\in\mathbb{R}$
govern the couplings between the fast and slow systems, and $F>0$ provides a
constant forcing. Note that the equation (11a) has the form of the singlescale
Lorenz model (6) with a coupling to the fast variables replacing the function
$m(\cdot).$ We set $L=9,\ J=8,\ h_{v}=-0.8,\ h_{w}=1,\ F=10$ and
$\epsilon=2^{-7}$ in what follows. In Figure 5 we display the dynamics of a
slow variable and one associated fast variable of this Lorenz ’96 multiscale
model; at the parameter values chosen the system is chaotic.

We consider observations $\\{{y}^{\dagger}_{n}\\}_{n\in\mathbb{Z}^{+}}$
arising from the model

$\displaystyle({v}^{\dagger}_{n+1},{w}^{\dagger}_{n+1})$
$\displaystyle=\Psi_{\rm{mult}}({v}^{\dagger}_{n},{w}^{\dagger}_{n}),$
$\displaystyle{y}^{\dagger}_{n+1}$ $\displaystyle=h({v}^{\dagger}_{n+1}),$

where $\Psi_{\rm{mult}}$ is the solution operator to the multiscale model
(11)–(13) over the observation time interval $\tau$. We then take the data
from the multiscale model and assimilate it into the singlescale model (6),
recalling the relationship between them detailed in Appendix 7, and the
importance of the specific function $m$ shown in Figure 2, to the elimination
of the fast variables in favour of a simple closure, using the averaging
principle. As in the preceding Example 2.1, we assume that the observation
function is linear: $h(v)=Hv$ for matrix $H:\mathbb{R}^{9}\to\mathbb{R}^{6}$
defined by (7). As before we choose the gain
$K:\mathbb{R}^{6}\to\mathbb{R}^{9}$ to be defined by (9) and employ the 3DVAR
algorithm (8) with $\Psi$ the solution operator over time-interval $\tau$ for
the singlescale model.

Figure 5: Dynamics of a slow variable and an associated fast variable. Here
“dns”, in blue, labels the slow variable computed by direct numerical
simulation; in grey is the related fast variable.

Again we display the results of 3DVAR on the unobserved component $v_{3}.$
Figure 6(a) shows the behavior of the algorithm for $\tau=10^{-3}.$ Despite
the fact that the data is generated from the multiscale model, whilst
assimilation is conducted using the singlescale model, 3DVAR produces an
accurate estimate of the true state with, after synchronization, the only
discernible errors being slight under and overshoots. Figure 6(b) shows the
effect of varying $\tau$, the time between observations. As in Example 2.1 the
assimilation breaks down when $\tau$ is larger. $\blacksquare$

(a)

(b)

Figure 6: In (a) the estimates of $v_{3}$ in time produced by 3DVAR using
$\tau=10^{-3}$ are displayed against the true dynamics. In (b) the estimates
of 3DVAR are shown for $\tau=10^{0}$. Again the acronym “dns” refers to direct
numerical simulation. 3DVAR successfully synchronizes with the direct
numerical simulation at the smaller value of $\tau$ but fails to do so when
$\tau$ is larger. It is noteworthy that the synchronization takes place here
in the context of model misspecification: the data is generated by the
multiscale model, but 3DVAR is applied using the singlescale model.

The control theoretic approach encapsulated in 3DVAR addresses Objective 1.
However, when $|C_{0}|,|\Gamma|$ and $|\Sigma|$ are no longer necessarily
small, so that the state and observations are subject to noise, it is natural
to try and generalize the approach to address Objective 2, and include non-
zero covariances; Example 2.1 shows that this may be needed when noise is
larger. A natural stochastic generalization of the observer approach is as
follows: from current state estimate $v_{n}$, given ${Y}^{\dagger}_{n}$,
predict the outcome of the model and data from the update equations (1), which
we denote by $(\widehat{v}_{n+1},\widehat{y}_{n+1})$; then correct the state
estimate by nudging the mean of the prediction using the mismatch between
observed and predicted data $({y}^{\dagger}_{n+1},\widehat{y}_{n+1})$. Given
$v_{n}$ computed from ${Y}^{\dagger}_{n}$ this results in state estimate
$v_{n+1}$ from ${Y}^{\dagger}_{n+1}$ defined through the following stochastic
dynamical system, assumed to hold for all $n\in\mathbb{Z}^{+}:$

(14a) $\displaystyle\widehat{v}_{n+1}$ $\displaystyle=\Psi(v_{n})+\xi_{n},$
(14b) $\displaystyle\widehat{y}_{n+1}$
$\displaystyle=h(\widehat{v}_{n+1})+\eta_{n+1},$ (14c) $\displaystyle v_{n+1}$
$\displaystyle=\widehat{v}_{n+1}+K_{n}I_{n};$

here $v_{0},\xi_{n}$ and $\eta_{n+1}$ are random variables given by the known
distributions in (2) and $I_{n}$ is modified from (5) to read

(15) $I_{n}={y}^{\dagger}_{n+1}-\widehat{y}_{n+1}.$

The data $\\{{y}^{\dagger}_{n+1}\\}_{n\in\mathbb{Z}^{+}}$ is a fixed
realization of (1). Given this data, equations (14) then define a random map
$v_{n}\mapsto v_{n+1}$ which uses knowledge of the model and the observed data
to update our state estimate. The predicted state and data
$(\widehat{v}_{n+1},\widehat{y}_{n+1})$ are also referred to as the simulated
state and data. Choice of the gain matrix $K_{n}$ will complete definition of
an algorithm. The key question we proceed to study in subsequent sections
concerns the choice of $\\{K_{n}\\}_{n\in\mathbb{Z}^{+}}.$ We do this in the
context of addressing Objective 2. Before doing so we continue our running
example, studying the effect of noise on the 3DVAR algorithm, where $K_{n}=K$
is constant; this example demonstrates the need for an adaptive choice of
$K_{n}$ when noise is significant.

######  Example 2.3.

We return to the setting of Example 2.1 and consider the effect of including
noise in 3DVAR as in (14), (15). We again study the Lorenz ’96 singlescale
model (6) for unknown $v\in C(\mathbb{R}^{+},\mathbb{R}^{L})$, with
$L=9,h_{v}=-0.8$ and $F=10.$ We let $\Psi$ denote the solution operator for
(6) over the observation time interval $\tau$ and consider observations
$\\{{y}^{\dagger}_{n}\\}_{n\in\mathbb{Z}^{+}}$ defined by

$\displaystyle{v}^{\dagger}_{n+1}$
$\displaystyle=\Psi({v}^{\dagger}_{n})+{\xi}^{\dagger}_{n},$
$\displaystyle{y}^{\dagger}_{n+1}$
$\displaystyle=h({v}^{\dagger}_{n+1})+{\eta}^{\dagger}_{n+1},$

where $\\{{\xi}^{\dagger}_{n}\\}_{n\in\mathbb{Z}^{+}}$,
$\\{{\eta}^{\dagger}_{n}\\}_{n\in\mathbb{N}}$ are mutually independent
Gaussian sequences

${\xi}^{\dagger}_{n}\sim\mathsf{N}(0,\sigma^{2}I)\,\,\text{i.i.d.},\quad{\eta}^{\dagger}_{n}\sim\mathsf{N}(0,\gamma^{2}I)\,\,\text{i.i.d.}\,\,,$

with $\sigma^{2}=0.1$ and $\gamma^{2}=0.1$. We again assume that the
observation function is linear: $h(v)=Hv$ for matrix
$H:\mathbb{R}^{9}\to\mathbb{R}^{6}$ defined by (7). As in Example 2.1 we
choose fixed gain $K_{n}\equiv K$ with $K:\mathbb{R}^{6}\to\mathbb{R}^{9}$
defined by (9).

Figure 7 illustrates that this version of noisy 3DVAR does not perform better
than the noise-free 3DVAR (8), in the setting where true state and true
observation noise levels are high. To quantitatively demonstrate this, we
compute the mean squared error between the estimates arising from the 3DVAR
algorithm (8) and the true states, and the mean squared error between the
estimates obtained using noisy 3DVAR algorithm (14), (15) and the true states.
Given either method the error $e$ is computed using the following formula, in
which, recall, $\\{{v}^{\dagger}_{n}\\}$ is the truth and $\\{v_{n}\\}$ is the
output of the 3DVAR or noisy 3DVAR algorithm:

(16) $e=\frac{1}{N\cdot
d_{v}}\sum_{n=1}^{N}\big{|}{v}^{\dagger}_{n+t/\tau}-v_{n+t/\tau}\big{|}^{2}$

where $N$ is such that $t+N\tau=T$. Time $t$ is chosen to remove error from
the incorrect initialization, and focus on quantifying error in statistical
steady state; $T$ is chosen to allow time-averaging over a long enough window
to capture this statistical steady state. In this case, we compute this
average for the estimates obtained after time $t=5$, up to $T=30$, and report
$e_{\text{3DVAR}}=1.47\cdot 10^{0}$ and $e_{\text{noisy3DVAR}}=3.07\cdot
10^{0}$. Clearly use of noisy 3DVAR does not improve the error, in comparison
with noise-free 3DVAR. $\blacksquare$

Figure 7: In this experiment we set the noise levels
$\sigma^{2}=10^{-1},\gamma^{2}=10^{-1}$. Again the acronym “dns” refers to
direct numerical simulation. We display the estimates of $v_{3}$ in time
produced by 3DVAR and noisy 3DVAR against the true dynamics using
$\tau=10^{-3}$.

Example 2.3 highlights the need to quantify uncertainty and pass to a
probabilistic interpretation (Objective 2) of the filtering problem; and, in
particular, to make an informed choice of adaptive gain matrices $K_{n}.$ We
turn to the probabilistic interpretation in the next subsection; in subsequent
subsections we derive algorithms, leading in particular to a specific choice
of adaptive gain matrices.

###  2.3 Probabilistic Perspective

We have shown that the 3DVAR methodology can recover the state of a (possibly
chaotic) dynamical system, even though the initial condition is not known, by
exploiting the observations. However it does not quantify uncertainty in the
state estimate. We now introduce a probabilistic perspective which enables us
to address this issue. In Subsection 2.3.1 we define the filtering
distribution and in Subsection 2.3.2 we establish some notation used
throughout the sequel.

####  2.3.1 The Filtering Distribution

To understand how uncertainty may be expressed mathematically we introduce the
filtering distribution $\mu_{n}$: the distribution of the conditioned random
variable $v_{n}|{Y}^{\dagger}_{n}$. This captures the knowledge of the state
of the system, and uncertainties in the state, given the observations. It is
important to understand how this distribution evolves with $n$. In particular,
the evolution has a _nonlinear_ Markov structure; we emphasize that standard
Markov processes are linear, and defined by the linear transition kernel;
nonlinearity arises in filtering through the incorporation of data. The map
$\mu_{n}\mapsto\mu_{n+1}$ is most easily described by defining
$\widehat{\mu}_{n+1}$, the distribution of $v_{n+1}|{Y}^{\dagger}_{n}$ and
$\nu_{n+1}$, the distribution of $(v_{n+1},y_{n+1})|{Y}^{\dagger}_{n}.$ The
map from $\mu_{n}$ to $\widehat{\mu}_{n+1}$ is defined by equation (3a) and is
linear as a map from the space of probability measures defined on
$\mathbb{R}^{d_{v}}$ into itself; the map from $\widehat{\mu}_{n+1}$ to
$\nu_{n+1}$ is defined by (3b), and is linear as a map from the space of
probability measures defined on $\mathbb{R}^{d_{v}}$ into the space of
probability measures defined on $\mathbb{R}^{d_{v}+d_{y}}$; the map from
$\nu_{n+1}$ to $\mu_{n+1}$ is defined by conditioning $\nu_{n+1}$ on
${y}^{\dagger}_{n+1}$ and is nonlinear as a map from the space of probability
measures defined on $\mathbb{R}^{d_{v}+d_{y}}$ into the space of probability
measures defined on $\mathbb{R}^{d_{v}}.$ We thus have

(17a) $\displaystyle\widehat{\mu}_{n+1}$ $\displaystyle=P\mu_{n},$ (17b)
$\displaystyle\nu_{n+1}$ $\displaystyle=Q\widehat{\mu}_{n+1},$ (17c)
$\displaystyle\mu_{n+1}$ $\displaystyle=B_{n}(\nu_{n+1}),$

where $P$ and $Q$ are linear maps on measures, whilst $B_{n}(\cdot)$ is
nonlinear and time dependent. Concatenating we find that

(18) $\mu_{n+1}=B_{n}(QP\mu_{n}),\quad\mu_{0}=\mathsf{N}(m_{0},C_{0}).$

This map defines an inhomogeneous nonlinear Markov process on the space of
probability measures on $\mathbb{R}^{d_{v}}.$ Markovian because it maps
probability measures at time $n$ into probability measures at time $n+1$;
nonlinear because of the conditioning operator $B_{n}(\cdot)$; and
inhomogeneous because map $B_{n}$ depends on time $n$ through the time-varying
observation ${y}^{\dagger}_{n+1}.$ The map $B_{n}(QP\cdot)$ may be decomposed
into two maps: the prediction $P$, which represents application of the
dynamical model (1a), and application of Bayes Theorem 333Often referred to as
the analysis step in the geophysical data assimilation community. through
operator $L_{n}\cdot:=B_{n}(Q\cdot)$, which corresponds to use of likelihood
defined by the observation model (1b). With this notation we thus obtain

(19a) $\displaystyle\widehat{\mu}_{n+1}$ $\displaystyle=P\mu_{n},$ (19b)
$\displaystyle\mu_{n+1}$ $\displaystyle=L_{n}(\widehat{\mu}_{n+1}).$

We refer to iteration of (19) as the filtering cycle. The cycle involves
iterative interleaving of prediction, using the dynamical model, a linear
operation on measures, and Bayes Theorem, using the observation model, a
nonlinear operation on measures.

It is useful, notationally, to work with measures
$\mu_{n},\widehat{\mu}_{n},\nu_{n}$ rather than their Lebesgue densities; this
facilitates clean presentation of particle approximations which employ Dirac
measures. However, to define the maps $P,Q$ and $L_{n}$, we use probability
densities arising from the Gaussian nature of the noise appearing in (1),
summarized in (2). To this end, note that (1a) is defined by Markov transition
kernel $p(u,v)$, and the probability of $y|v$ underlying (1b) is defined by
kernel $q(v,y)$. When $\Sigma\succ 0,\Gamma\succ 0$, these are given by

(20a) $\displaystyle p(u,v)$
$\displaystyle=\frac{1}{(2\pi)^{d_{v}/2}\sqrt{{\rm
det}(\Sigma)}}\exp\Bigl{(}-\frac{1}{2}|v-\Psi(u)|_{\Sigma}^{2}\Bigr{)},$ (20b)
$\displaystyle q(v,y)$ $\displaystyle=\frac{1}{(2\pi)^{d_{y}/2}\sqrt{{\rm
det}(\Gamma)}}\exp\Bigl{(}-\frac{1}{2}|y-h(v)|_{\Gamma}^{2}\Bigr{)}.$

Then we have

$\displaystyle(P\mu)(dv)$
$\displaystyle=\Bigl{(}\int_{u\in\mathbb{R}^{d_{v}}}p(u,v)\mu(du)\Bigr{)}dv,$
$\displaystyle(Q\mu)(dv,dy)$ $\displaystyle=q(v,y)\mu(dv)dy,$
$\displaystyle(B_{n}\nu)(dv)$
$\displaystyle=\int_{y\in\mathbb{R}^{d_{y}}}\delta_{{y}^{\dagger}_{n+1}}(y)\nu(dv,dy)\Big{/}\Bigl{(}\int_{(v,y)\in\mathbb{R}^{d_{v}}\times\mathbb{R}^{d_{y}}}\delta_{{y}^{\dagger}_{n+1}}(y)\nu(dv,dy)\Bigr{)}$
$\displaystyle L_{n}(\mu)(dv)$
$\displaystyle=q(v,{y}^{\dagger}_{n+1})\mu(dv)\Big{/}\Bigl{(}\int_{v\in\mathbb{R}^{d_{v}}}q(v,{y}^{\dagger}_{n+1})\mu(dv)\Bigr{)}.$

It is important to appreciate that there is, in general, no closed form
expression for $\mu_{n}$ defined by the iteration (18); thus (18) does not
constitute an algorithm. However if $\Psi$ and $h$ are linear then, since
$\mu_{0}=\mathsf{N}(m_{0},C_{0})$ is Gaussian, it follows that
$\widehat{\mu}_{n+1},\nu_{n+1},\mu_{n+1}$ are all Gaussian for all
$n\in\mathbb{Z}^{+}$, and closed form expressions, based on dynamical updates
of means and covariances, are available; this is the Kalman filter discussed
in Subsection 2.4.2. Furthermore, the idea of making a Gaussian approximation
will play a central role in ensemble Kalman methods, and this perspective is
developed in Subsection 2.4.

####  2.3.2 Important Repeatedly Used Notation

Throughout this subsection let $v_{n}$ be a random variable distributed
according to $\textit{Law}\,(v_{n})$, let $\widehat{v}_{n+1}$ be random
variable with $\textit{Law}\,(\widehat{v}_{n+1})=P\textit{Law}\,(v_{n})$ and
let $(\widehat{v}_{n+1},\widehat{y}_{n+1})$ be random variable with law
$Q\textit{Law}\,(\widehat{v}_{n+1}).$ Using Law avoids a proliferation of
notation for the different measures arising from the use of various different
algorithms to approximate the filtering cycle.

All of the algorithms are based on a prediction of state, and possibly
observation, from $\textit{Law}\,(v_{n})$. To this end recall (14) and for all
$n\in\mathbb{Z}^{+}$,

(21a) $\displaystyle\widehat{v}_{n+1}$ $\displaystyle=\Psi(v_{n})+\xi_{n},$
(21b) $\displaystyle\widehat{y}_{n+1}$
$\displaystyle=h(\widehat{v}_{n+1})+\eta_{n+1}.$

The distribution of $(\widehat{v}_{n+1},\widehat{y}_{n+1})$ given by these
equations defines $P\textit{Law}\,(v_{n})$ and
$Q\textit{Law}\,(\widehat{v}_{n+1}).$ A common theme in this paper is to
augment these equations for the predicted state and data with a final map

(22) $\displaystyle\widehat{v}_{n+1}\mapsto v_{n+1},\,\,{\rm
or}\,\,(\widehat{v}_{n+1},\widehat{y}_{n+1})\mapsto v_{n+1}$

in such a way that, if $v_{n}\sim\mu_{n}$ then $v_{n+1}\sim\mu_{n+1}$; the
algorithms we study will be based on approximations of such maps. We refer to
equations (21,22) as the sample path perspective, designed to exactly or
approximately reproduce the probabilistic perspective encapsulated in the
three steps of (17). Similarly we will also introduce equations of the form
(21a,22) to provide sample path exact, or approximate, reproduction of the
probabilistic perspective encapsulated in the two steps of (19).

In service of designing these approximations, it is useful to define various
first and second order statistics computed under the law of
$(\widehat{v}_{n+1},\widehat{y}_{n+1})$. First define mean and covariance of
$\widehat{v}_{n+1}:$

(23a) $\displaystyle\widehat{m}_{n+1}$
$\displaystyle=\mathbb{E}\widehat{v}_{n+1},$ (23b)
$\displaystyle\widehat{C}_{n+1}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}\widehat{v}_{n+1}-\widehat{m}_{n+1})\otimes\bigl{(}\widehat{v}_{n+1}-\widehat{m}_{n+1}\bigr{)}\Bigr{)}.$

Then define mean of the predicted data, cross-covariance from predicted data
to state and covariance of the data:

(24a) $\displaystyle\widehat{o}_{n+1}$
$\displaystyle=\mathbb{E}\widehat{y}_{n+1},$ (24b)
$\displaystyle\widehat{C}_{n+1}^{vy}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}\widehat{v}_{n+1}-\widehat{m}_{n+1}\bigr{)}\otimes\bigl{(}\widehat{y}_{n+1}-\widehat{o}_{n+1}\bigr{)}\Bigr{)},$
(24c) $\displaystyle\widehat{C}_{n+1}^{yy}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}\widehat{y}_{n+1}-\widehat{o}_{n+1}\bigr{)}\otimes\bigl{(}\widehat{y}_{n+1}-\widehat{o}_{n+1}\bigr{)}\Bigr{)}.$

From these covariances we define the matrix

(25) $\displaystyle
K_{n}=\widehat{C}^{vy}_{n+1}(\widehat{C}^{yy}_{n+1})^{-1}.$

This particular choice of $K_{n}$, known as the _Kalman gain_ , is required
for definition of the Kalman filter, and will play a central role in the mean
field maps which underpin ensemble Kalman methods through their particle
approximations. The origin of this key concept will be described in the next
subsection. Note that, if $\widehat{C}^{yy}_{n+1}$ is not invertible, then its
action may still be defined through a pseudo-inverse.

We also note that, in view of (21b),
$\widehat{o}_{n+1}=\mathbb{E}h(\widehat{v}_{n+1}).$ It is sometimes useful to
express the Kalman gain $K_{n}$ in terms of the variable
$\widehat{h}_{n+1}=h(\widehat{v}_{n+1})$. For this purpose we define the
associated correlation matrices

(26a) $\displaystyle\widehat{C}_{n+1}^{vh}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}\widehat{v}_{n+1}-\widehat{m}_{n+1}\bigr{)}\otimes\bigl{(}\widehat{h}_{n+1}-\widehat{o}_{n+1}\bigr{)}\Bigr{)},$
(26b) $\displaystyle\widehat{C}_{n+1}^{hh}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}\widehat{h}_{n+1}-\widehat{o}_{n+1}\bigr{)}\otimes\bigl{(}\widehat{h}_{n+1}-\widehat{o}_{n+1}\bigr{)}\Bigr{)},$

and note that

(27a) $\displaystyle\widehat{C}_{n+1}^{vy}$
$\displaystyle=\widehat{C}_{n+1}^{vh},\quad\widehat{C}_{n+1}^{yy}=\widehat{C}_{n+1}^{hh}+\Gamma,$
(27b) $\displaystyle K_{n}$
$\displaystyle=\widehat{C}_{n+1}^{vh}\left(\widehat{C}_{n+1}^{hh}+\Gamma\right)^{-1}.$

Note that, if $\Gamma\succ 0$, $K_{n}$ is well-defined without recourse to the
use of pseudo-inverse.

###  2.4 Gaussian Projected Filtering Distribution

The Gaussian projected filtering distribution is an approximation of the true
filtering which leads to a filter which is Gaussian. It is defined by the
following steps: (i) taking input Gaussian at time $n$ as
$\textit{Law}\,(v_{n})$ and pushing this measure forward under (21) to find a
(typically non-Gaussian) measure for
$\textit{Law}\,(\widehat{v}_{n+1},\widehat{y}_{n+1})$; (ii) projecting this
joint measure onto the nearest Gaussian (in a sense that we will make
precise); (iii) then finally conditioning this Gaussian on the data
${y}^{\dagger}_{n+1}$ to find output Gaussian at time $n+1.$ We note that
conditioning a Gaussian random variable on linear functionals of the random
variable returns another Gaussian.

The resulting approximation of the filtering distribution plays an important
role in motivating the mean field maps, introduced in Section 2.5, whose
particle approximations lead to various ensemble Kalman methods introduced in
Section 2.6. In Subsection 2.4.1 we introduce the Gaussian projected
approximation to the evolution (18). In the case where $\Psi(\cdot)$ and
$h(\cdot)$ are linear the resulting formulae deliver exact solutions of the
filtering cycle (18), leading to the Kalman filter as derived in Subsection
2.4.2. The Gaussian projected filter leads to a derivation of the Kalman gain
(25); an alternative derivation, using the minimum variance approach, may be
found in Appendix C, Subsection 8.3.

####  2.4.1 Gaussian Projected Filtering

To describe the Gaussian projected approximation to the filtering
distribution, we define the following map $G$:

######  Definition 2.4.

Let $\mathfrak{P}=\mathfrak{P}(\mathbb{R}^{d})$ denote the set of probability
measures on $\mathbb{R}^{d}$ and let
$\mathfrak{G}=\mathfrak{G}(\mathbb{R}^{d})$ denote the set of Gaussian
probability measures on $\mathbb{R}^{d}$ (including indefinite covariances and
hence all Dirac masses) and define $G:\mathfrak{P}\mapsto\mathfrak{G}$, for
$u\sim\mu$, by

$G\mu=N(m^{\mu},C^{\mu}),\qquad m^{\mu}=\mathbb{E}^{\mu}u,\quad
C^{\mu}=\mathbb{E}^{\mu}\bigl{(}(u-\mathbb{E}^{\mu}u)\otimes(u-\mathbb{E}^{\mu}u)\bigr{)}.$

Thus the map $G$ applied to measure $\mu$ simply computes the Gaussian with
mean and covariance under $\mu.$ Notice that $G$ is the identity on Gaussians.
Furthermore $G\circ G=G.$ We refer to this as a projection onto Gaussians
because it corresponds to finding the closest point to given measure $\mu$,
with respect to a Kullback-Leibler divergence444For details see the
bibliography Subsection 2.7:

$G\mu={\rm argmin}_{\pi\in\mathfrak{G}}d_{\rm{KL}}(\mu||\pi).$

We now use mapping $G$ to find an approximation to the evolution (17) which
generates measures which remain Gaussian; intuitively this will be a good
approximation whilst the measures $\\{\mu_{n}\\}$ evolving under (18) remain
close to Gaussian. To this end we consider random variable
$v_{n}\sim\mu^{G}_{n}$ where the probability measure $\mu^{G}_{n}$ evolves
according to

(28)
$\mu^{G}_{n+1}=B_{n}(GQP\mu^{G}_{n}),\quad\mu^{G}_{0}=\mathsf{N}(m_{0},C_{0}).$

This may be decomposed as follows:

(29a) $\displaystyle\widehat{\mu}^{G}_{n+1}$ $\displaystyle=P\mu^{G}_{n},$
(29b) $\displaystyle\nu^{G}_{n+1}$ $\displaystyle=Q\widehat{\mu}^{G}_{n+1},$
(29c) $\displaystyle\mu^{G}_{n+1}$ $\displaystyle=B_{n}(G\nu^{G}_{n+1}).$

Map (28) defines a nonlinear Markov process, similarly to (18). It also maps
Gaussians into Gaussians. This fact follows from the fact that the nonlinear
map $B_{n}(\cdot)$, which represents conditioning, maps Gaussians into
Gaussians. The map $\mu^{G}_{n}\mapsto\mu^{G}_{n+1}$ hence defines a
deterministic mapping from the mean $m_{n}$ and covariance $C_{n}$ of
$\mu^{G}_{n}$ into the mean $m_{n+1}$ and covariance $C_{n+1}$ of
$\mu^{G}_{n+1}.$ We now identify this map explicitly.

For $v_{n}\sim\mu^{G}_{n}$ we introduce the random variables
$\widehat{v}_{n+1},\widehat{y}_{n+1}$ defined by (21). It then follows that
$\widehat{v}_{n+1}\sim\widehat{\mu}^{G}_{n+1}=P\mu^{G}_{n}$ and that
$(\widehat{v}_{n+1},\widehat{y}_{n+1})\sim\nu^{G}_{n+1}=QP\mu^{G}_{n}.$ Note
also that $\widehat{\mu}^{G}_{n+1}$ and $\nu^{G}_{n+1}$ are not Gaussian, but
are defined by the Gaussian $\mu^{G}_{n}$ and hence completely determined by
$m_{n}$ and $C_{n}.$ Thus the mean and covariance under
$\widehat{\mu}^{G}_{n+1}$ are $(\widehat{m}_{n+1},\widehat{C}_{n+1})$ given by
(23). Furthermore, $G\nu^{G}_{n+1}$ is defined by

$G\nu^{G}_{n+1}=\mathsf{N}\Bigl{(}\begin{bmatrix}\widehat{m}_{n+1}\\\
\widehat{o}_{n+1}\end{bmatrix},\begin{bmatrix}\widehat{C}_{n+1}&\widehat{C}_{n+1}^{vy}\\\
\bigl{(}\widehat{C}_{n+1}^{vy}\bigr{)}^{\top}&\widehat{C}_{n+1}^{yy}\end{bmatrix}\Bigr{)},$

where all relevant quantities are defined in Subsection 2.3.2.

We now condition the Gaussian $G\nu^{G}_{n+1}$, on the second component of the
vector taking value ${y}^{\dagger}_{n+1}$. From this we find the Gaussian
measure $\mu^{G}_{n+1}=B_{n}(G\nu^{G}_{n+1})$ characterized by mean $m_{n+1}$
and covariance $C_{n+1}$ given by the following lemma:

######  Lemma 2.5.

Assume that $\Gamma\succ 0.$ Let $m_{n}$ and $C_{n}$ denote the mean and
covariance under the Gaussian projected filter. Consider
$\widehat{m}_{n+1},\widehat{C}_{n+1}$ given by equations (21), (23), and the
mean of the observed data and covariances given by (24). Then
$\widehat{C}_{n+1}^{yy}\succ 0$ for all $n\in\mathbb{Z}^{+}$ and

(30a) $\displaystyle m_{n+1}$
$\displaystyle=\widehat{m}_{n+1}+\widehat{C}_{n+1}^{vy}(\widehat{C}_{n+1}^{yy})^{-1}\bigl{(}{y}^{\dagger}_{n+1}-\widehat{o}_{n+1}\bigr{)},$
(30b) $\displaystyle C_{n+1}$
$\displaystyle=\widehat{C}_{n+1}-\widehat{C}_{n+1}^{vy}(\widehat{C}_{n+1}^{yy})^{-1}\bigl{(}\widehat{C}_{n+1}^{vy}\bigr{)}^{\top},$

where $\\{{y}^{\dagger}_{n}\\}$ arises from a fixed realization of (1).

######  Proof 2.6.

We first note that $\widehat{C}_{n+1}^{yy}\succ 0$. Indeed, since by
assumption $\Gamma\succ 0$ and by definition $\widehat{C}_{n+1}^{hh}\succeq
0$, then $\widehat{C}_{n+1}^{hh}+\Gamma\succ 0$. Hence by (27) we have that
$\widehat{C}_{n+1}^{yy}\succ 0$. By explicitly formulating the distribution of
the Gaussian $G\nu^{G}_{n+1}$ conditioned on ${y}^{\dagger}_{n+1}$, it is
possible to conclude that from standard formulae for conditioned Gaussians
that $m_{n+1}$ and $C_{n+1}$ are given by the expressions in (30).

Equations (21), (23), (24) and (30) define a mapping from $\mu^{G}_{n}$,
characterized by $(m_{n},C_{n})$, into $\mu^{G}_{n+1}$, characterized by
$(m_{n+1},C_{n+1})$. They comprise an explicit set of formulae for the mapping
(28): since Gaussians are determined by mean and covariance, the map on
measures is completely determined by the map from $(m_{n},C_{n})$ to
$(m_{n+1},C_{n+1})$.

We note that (30) can also be rewritten as

(31a) $\displaystyle m_{n+1}$
$\displaystyle=\widehat{m}_{n+1}+\widehat{C}_{n+1}^{vh}(\widehat{C}_{n+1}^{hh}+\Gamma)^{-1}\bigl{(}{y}^{\dagger}_{n+1}-\widehat{o}_{n+1}\bigr{)},$
(31b) $\displaystyle C_{n+1}$
$\displaystyle=\widehat{C}_{n+1}-\widehat{C}_{n+1}^{vh}(\widehat{C}_{n+1}^{hh}+\Gamma)^{-1}\bigl{(}\widehat{C}_{n+1}^{vh}\bigr{)}^{\top}.$

Equations (21), (23), (26) and (31) then also define the mapping from
$\mu^{G}_{n}$ into $\mu^{G}_{n+1}$ and also comprise an explicit set of
formulae for the updates of the mean and covariance which characterize mapping
(28).

Finally we note that, using the definition (25) of Kalman gain, we can rewrite
(30) as

$\displaystyle m_{n+1}$
$\displaystyle=\widehat{m}_{n+1}+K_{n}\bigl{(}{y}^{\dagger}_{n+1}-\widehat{o}_{n+1}\bigr{)},$
$\displaystyle C_{n+1}$
$\displaystyle=\widehat{C}_{n+1}-K_{n}\bigl{(}\widehat{C}_{n+1}^{vy}\bigr{)}^{\top}.$

######  Remark 2.7.

The preceding explicit formulae for (28) involve the computation of
expectations under the non-Gaussian measure
$\textit{Law}\,(\widehat{v}_{n+1},\widehat{y}_{n+1})$. For this reason they do
not constitute an algorithm. A possible approach to algorithmic
implementations involves quadrature to approximate the expectations, leading,
for example, to the _unscented Kalman filter_ approach; details may be found
in the bibliography Subsection 2.7. However the formulation of explicit maps
on Gaussians plays a different, important, role in this paper: we use it as a
way of explaining the sense in which the distribution of our mean-field models
approximate evolution of measures under the true filtering distribution; see
Subsection 2.5.5. The explicit map on means and covariances also enables us to
derive the Kalman filter, which applies in the linear Gaussian setting, the
topic of the next subsection. $\blacksquare$

####  2.4.2 The Kalman Filter

In the setting where $\Psi(\cdot)=M\cdot$ and $h(\cdot)=H\cdot$ are both
linear, (1) becomes

(32a) $\displaystyle v_{n+1}$ $\displaystyle=Mv_{n}+\xi_{n},$ (32b)
$\displaystyle y_{n+1}$ $\displaystyle=Hv_{n+1}+\eta_{n+1}.$

Furthermore, the mappings (18) and (28) are identical and the Gaussian
projected filter is exact. Then, equations (23) give

(33a) $\displaystyle\widehat{m}_{n+1}$ $\displaystyle=Mm_{n},$ (33b)
$\displaystyle\widehat{C}_{n+1}$ $\displaystyle=MC_{n}M^{\top}+\Sigma;$

and (24) and (30) give

(34a) $\displaystyle m_{n+1}$
$\displaystyle=\widehat{m}_{n+1}+\widehat{C}_{n+1}H^{\top}(H\widehat{C}_{n+1}H^{\top}+\Gamma)^{-1}({y}^{\dagger}_{n+1}-H\widehat{m}_{n+1}),$
(34b) $\displaystyle C_{n+1}$
$\displaystyle=\widehat{C}_{n+1}-\widehat{C}_{n+1}H^{\top}(H\widehat{C}_{n+1}H^{\top}+\Gamma)^{-1}H\widehat{C}_{n+1},$

where $\\{{y}^{\dagger}_{n}\\}$ arises from a fixed realization of (32). The
update equations (33), (34) constitute the Kalman filter, which applies to (1)
in the linear setting, mapping $(m_{n},C_{n})$ to $(m_{n+1},C_{n+1}).$ The
equations as stated are well-defined if $\Gamma\succ 0.$

###  2.5 Mean Field Maps

Note that elements of $\mathfrak{P}(\mathbb{R}^{d})$ are infinite dimensional
objects. The preceding subsection provides explicit finite dimensional maps
which characterize the Gaussian projected filter
$\mu^{G}_{n}\mapsto\mu^{G}_{n+1}.$ This is possible because Gaussians are
completely determined by a finite amount of information – the mean and
covariance. In this section we have a more ambitious aim: to find maps on
finite dimensional spaces with the property that (possibly only approximately)
their output is equal in law to the map on measures $\mu_{n}\mapsto\mu_{n+1}$
given by the filtering cycle. We achieve this by studying transport maps that
achieve (22), or approximate it; this leads us to the subject of mean field
maps, namely random maps that depend on the law of the state being mapped.

Consider probability measures $\nu$ and $\nu^{\prime}$ on $\mathbb{R}^{d}$ and
$\mathbb{R}^{d^{\prime}}$ respectively, and recall the definition of
pushforward of a measure under under
$T:\mathbb{R}^{d}\to\mathbb{R}^{d^{\prime}}$: the statement
$\nu^{\prime}=T^{\sharp}\nu$ is a succinct way of stating that, if
$\textit{Law}\,(v)=\nu$ and $v^{\prime}=T(v)$, then
$\textit{Law}\,(v^{\prime})=\nu^{\prime}.$ Given probability measures $\pi$
and $\pi^{\prime}$ on $\mathbb{R}^{d}$ and $\mathbb{R}^{d^{\prime}}$
respectively, a transport $T:\mathbb{R}^{d}\to\mathbb{R}^{d^{\prime}}$ is a
map with property that the pushforward of probability measure $\pi$ under $T$,
$T^{\sharp}\pi$, is equal to probability measure $\pi^{\prime}.$ In the
following we refer to $\pi$ as the source measure and $\pi^{\prime}$ as the
target measure defining the transport. In our setting $\pi^{\prime}$ will be
uniquely determined by $\pi$, together with an observed piece of finite
dimensional data. Thus $T$ depends on $\pi$ and so we may view $T$ as a
mapping $\mathbb{R}^{d}\times\mathfrak{P}\to\mathbb{R}^{d},$ suppressing, for
the moment, explicit dependence on the observed data. We can compute the
pushforward of $T$ on any measure in $\mathfrak{P}$; but when we compute the
pushforward on $\pi$ we obtain $\pi^{\prime}.$ We emphasize that transport
maps are not uniquely defined by their source and target measures and require
certain conditions for their existence, which we assume here to be satisfied.
The underlying mathematical concept is that of coupling of measures. See
Subsection 2.7 for discussion of transport, optimal transport and coupling.

We will also consider classes of _approximate transport maps_ which do not
achieve transport from $\pi$ to $\pi^{\prime}$, but instead match second order
moments information (we will be precise below). Such maps will also depend on
$\pi.$ To be completely clear, the dependence of (possibly approximate)
transport $T$ on a measure in $\mathfrak{P}$ does not affect the definition of
pushforward; we employ the following general definition of pushforward for
measure-dependent maps, taken to hold for all $\pi_{1},\pi_{2}$ regardless of
any assumed relationship between them:

(35)
$T(\cdot;\pi_{1})^{\sharp}\pi_{2}=T(\cdot;\widetilde{\pi})^{\sharp}\pi_{2}\Bigr{|}_{\widetilde{\pi}=\pi_{1}};\\\
$

in particular,

(36a) $\displaystyle T(\cdot;\pi)^{\sharp}\pi$
$\displaystyle=T(\cdot;\widetilde{\pi})^{\sharp}\pi\Bigr{|}_{\widetilde{\pi}=\pi},$
(36b) $\displaystyle T(\cdot;\pi)^{\sharp}(G\pi)$
$\displaystyle=T(\cdot;\widetilde{\pi})^{\sharp}(G\pi)\Bigr{|}_{\widetilde{\pi}=\pi}.$

In the preceding, pushforward under $T(\cdot,\widetilde{\pi})$ denotes regular
pushforward with no relationship assumed between $\widetilde{\pi}$ and measure
being pushed forward.

The (approximate) transport maps just identified can be recast as mean field
maps when used in the context where source $\pi$ is the distribution of the
input to $T$. In this section we identify mean field maps which effect exact
filtering (transport maps), or which effect forms of approximate filtering
(approximate transport maps). Subsection 2.5.1 introduces two distinct
transport approaches which effect exact filtering: one based on the prior to
posterior map that constitutes the analysis (Bayesian inference) step of
filtering, a transport between probability measures on the same space; and the
other based on the conditioning component of the analysis step, a transport
between probability measures on different spaces. We refer to these maps which
effect exact filtering as perfect transport.555Perfect transport should not be
confused with _optimal transport_ which identifies among all (perfect)
transport maps the one minimizing a certain cost functional such as the
Wasserstein distance. See the bibliography Subsection 2.7 for more details. We
use _perfect_ here to distinguish from the _approximate_ transport maps, based
only on matching first and second moment; these approximate methods underpin
ensemble Kalman methods. Subsections 2.5.3 and 2.5.4 are concerned with
approximations of these two perfect transports. In these approximations the
pushforward under the transport map is designed to match only the first and
second order moments of the target measure. For this reason these
approximations are closely related to the previously defined Gaussian
projected filter; we will elaborate on this connection. Finally, Subsection
2.5.6 applies the ideas of this section to find a mean field formulation of
the Kalman filter for linear Gaussian systems.

####  2.5.1 Perfect Transport

Consider the idea of finding a transport map that acts on the joint space of
state and data, to effect conditioning with respect to the observed data. To
this end we consider the dynamical system, assumed to hold for all
$n\in\mathbb{Z}^{+}$:

(37a) $\displaystyle\widehat{v}_{n+1}$ $\displaystyle=\Psi(v_{n})+\xi_{n},$
(37b) $\displaystyle\widehat{y}_{n+1}$
$\displaystyle=h(\widehat{v}_{n+1})+\eta_{n+1},$ (37c) $\displaystyle v_{n+1}$
$\displaystyle=T^{S}(\widehat{v}_{n+1},\widehat{y}_{n+1};\nu_{n+1},{y}^{\dagger}_{n+1}),$

where $\\{{y}^{\dagger}_{n}\\}$ arises from a fixed realization of (1). The
first two equations, which coincide with (21), effect the mappings from
$\mu_{n}$ to $\widehat{\mu}_{n+1}$ and from $\widehat{\mu}_{n+1}$ to
$\nu_{n+1}.$ Map
$T^{S}_{n}(\cdot,\cdot):=T^{S}(\cdot,\cdot;\nu_{n+1},{y}^{\dagger}_{n+1})$ is
then an explicit example of (22) defined to effect the desired conditioning of
$\nu_{n+1}$ on ${y}^{\dagger}_{n+1}$ in order to obtain $\mu_{n+1}.$ The
letter $S$ in $T^{S}$ connotes the dependence of the map on the stochastic
data $\widehat{y}_{n+1}.$ This is a mean-field stochastic dynamical system
mapping $v_{n}$ to $v_{n+1}$: stochastic because of the noise in (37a, 37b),
mean-field because the map $T^{S}$ in (37c) depends on the law $\nu_{n+1}$ of
$(\widehat{v}_{n+1},\widehat{y}_{n+1})$, and hence on $\mu_{n}$. The three
update steps in this mean field stochastic dynamical system lead to the
following maps on measures:

(38a) $\displaystyle\widehat{\mu}_{n+1}$ $\displaystyle=P\mu_{n},$ (38b)
$\displaystyle\nu_{n+1}$ $\displaystyle=Q\widehat{\mu}_{n+1},$ (38c)
$\displaystyle\mu_{n+1}$ $\displaystyle=(T^{S}_{n})^{\sharp}\nu_{n+1}.$

This is simply a restatement of (17), noting that $(T^{S}_{n})^{\sharp}$ has
been chosen so that pushforward corresponds to conditioning $\nu_{n+1}$ on
data ${y}^{\dagger}_{n+1}$ to obtain $\mu_{n+1}.$ Note that the implied map
from $\mu_{n}$ to $\mu_{n+1}$ is a nonlinear Markov process, because of the
dependence of $T^{S}_{n}$ on $\nu_{n+1}$ and hence on $\mu_{n}.$ Furthermore,
we have that $\textit{Law}\,(v_{\ell})=\mu_{\ell}$ for all
$\ell\in\mathbb{Z}^{+}$.

Consider now a different approach to transport for filtering: we seek a
transport map that acts on the state space only to effect Bayes Theorem, i.e.
mapping prior $\widehat{\mu}_{n+1}$ to posterior $\mu_{n+1}$. To this end we
consider the dynamical system

(39a) $\displaystyle\widehat{v}_{n+1}$ $\displaystyle=\Psi(v_{n})+\xi_{n},$
(39b) $\displaystyle v_{n+1}$
$\displaystyle=T^{D}(\widehat{v}_{n+1};\widehat{\mu}_{n+1},{y}^{\dagger}_{n+1}),$

again assumed to hold for all $n\in\mathbb{Z}^{+}$. The first equation maps
$v_{n}\sim\mu_{n}$ to $\widehat{v}_{n+1}\sim\widehat{\mu}_{n+1}$, thus giving
a sample path realization of (19a). In the second equation, the map
$T^{D}_{n}(\cdot):=T^{D}(\cdot;\widehat{\mu}_{n+1},{y}^{\dagger}_{n+1})$ is
chosen so that, if $\widehat{v}_{n+1}\sim\widehat{\mu}_{n+1}$ then
$v_{n+1}\sim\mu_{n+1}$, thus giving a sample path realization of (19b). This
is also a mean-field stochastic dynamical system: stochastic because of the
noise in (39a); mean field because the map $T^{D}$ in (39b) depends on the law
of $\widehat{v}_{n+1}$ itself, and hence on $\mu_{n}.$ The symbol $D$
distinguishes map $T^{D}$ from map $T^{S}$: map $T^{D}$ is deterministic in
the sense that it does not require stochastic data $\widehat{y}_{n+1}$, in
contrast to $T^{S}.$ Therefore, we again have that
$\textit{Law}\,(v_{\ell})=\mu_{\ell}$ for all $\ell\in\mathbb{Z}^{+}$, where
the probability measure $\mu_{n}$ evolves according to

(40a) $\displaystyle\widehat{\mu}_{n+1}$ $\displaystyle=P\mu_{n},$ (40b)
$\displaystyle\mu_{n+1}$
$\displaystyle=(T^{D}_{n})^{\sharp}\widehat{\mu}_{n+1}.$

This is simply a restatement of (19), noting that $(T^{D}_{n})^{\sharp}$ has
been chosen so that the pushforward corresponds to the application of Bayes
Theorem to incorporate data ${y}^{\dagger}_{n+1}.$ The evolution (40) is
another nonlinear Markov process, now because of the dependence of $T^{D}_{n}$
on $\widehat{\mu}_{n+1}.$

The two transport maps $T^{S}$ and $T^{D}$ introduce an important conceptual
approach to algorithms for filtering, but determining the maps can be as hard,
or harder, than solving the filtering problem itself. Thus, in the next two
subsections, we turn to relaxations of the perfect transport effected by
$T^{S}$ and $T^{D}$. We instead seek mean field maps which match only first
and second order moment information; this relaxation allows for approximate
transport maps with simple affine (in senses to be made precise) forms. The
perspective of matching first and second order moments naturally suggests
working with Gaussians and hence we also relate the approximate transport to
the Gaussian projected filter.

####  2.5.2 Second Order Transport – Motivation

To motivate the more general ideas behind second order transport, we study an
explicit example. It is well known how to transform samples from a unit
centered Gaussian random variable on $\mathbb{R}$ into samples from a Gaussian
random variable with mean $m\not=0$ and variance $\sigma\not=1$ by a simple
scaling and shifting operation. An appropriate generalization of such a
procedure leads to the map

(41a) $\displaystyle\widehat{v}_{n+1}$ $\displaystyle=\Psi(v_{n})+\xi_{n},$
(41b) $\displaystyle\widehat{y}_{n+1}$
$\displaystyle=h(\widehat{v}_{n+1})+\eta_{n+1},$ (41c) $\displaystyle v_{n+1}$
$\displaystyle=m_{n+1}+C_{n+1}^{\frac{1}{2}}\widehat{C}_{n+1}^{-\frac{1}{2}}(\widehat{v}_{n+1}-\mathbb{E}\widehat{v}_{n+1}),$

where $m_{n+1}$, $\widehat{C}_{n+1}$, and $C_{n+1}$ are determined by (23),
(24) and (30), using (41a) and (41b). Here $\\{{y}^{\dagger}_{n}\\}$ arises
again from a fixed realization of (1). The map $v_{n}\mapsto v_{n+1}$ defined
by (41) is a mean-field map because of the dependence of $m_{n+1},C_{n+1}$ and
$\widehat{C}_{n+1}$ on $Q\textit{Law}\,(v_{n})$.

It is clear from (41c) that $\mathbb{E}v_{n+1}=m_{n+1}$ and also that

$\displaystyle\mathbb{E}\bigl{(}(v_{n+1}-m_{n+1})\otimes(v_{n+1}-m_{n+1})\bigr{)}$
$\displaystyle\quad\quad\quad\quad\quad=C_{n+1}^{\frac{1}{2}}\widehat{C}_{n+1}^{-\frac{1}{2}}\mathbb{E}\bigl{(}(v_{n+1}-m_{n+1})\otimes(v_{n+1}-m_{n+1})\bigr{)}\widehat{C}_{n+1}^{-\frac{1}{2}}C_{n+1}^{\frac{1}{2}},$
$\displaystyle\quad\quad\quad\quad\quad=C_{n+1}^{\frac{1}{2}}\widehat{C}_{n+1}^{-\frac{1}{2}}\widehat{C}_{n+1}\widehat{C}_{n+1}^{-\frac{1}{2}}C_{n+1}^{\frac{1}{2}},$
$\displaystyle\quad\quad\quad\quad\quad=C_{n+1}.$

It is important to recognize that $v_{n+1}$, defined by (41c), will not be
Gaussian distributed since $\widehat{v}_{n+1}$, defined by (41a), will not be
Gaussian either, in general. However, although (41) does not provide a closed
iteration on Gaussians, the map from $\widehat{v}_{n+1}$ to $v_{n+1}$ agrees
with the same step in the Gaussian projected filter, at the level of first and
second order moments.

Whilst the mean field map from (41) is a relatively transparent way to achieve
the goal of matching first and second order moments there is an uncountable
set of ways of achieving this objective; the next two subsections demonstrate
this, identifying all mean field maps effecting approximate transport from
within two specific classes of affine transformations. We will then highlight
a small subset that have been used in practice, each of which is useful in
certain specific contexts.

####  2.5.3 Second Order Transport – Stochastic Case

The first class of approximate filters determined by mean field maps have the
form

(42a) $\displaystyle\widehat{v}_{n+1}$ $\displaystyle=\Psi(v_{n})+\xi_{n},$
(42b) $\displaystyle\widehat{y}_{n+1}$
$\displaystyle=h(\widehat{v}_{n+1})+\eta_{n+1},$ (42c) $\displaystyle v_{n+1}$
$\displaystyle=\widetilde{T}^{S}(\widehat{v}_{n+1},\widehat{y}_{n+1};\nu_{n+1},{y}^{\dagger}_{n+1}),$

where $\\{{y}^{\dagger}_{n}\\}$ arises from a fixed realization of (1). We
define666We temporarily drop explicit notational dependence on $n+1$ in $\nu$
and in ${y}^{\dagger}$. This should not cause confusion as the approximate map
we derive is concerned simply with finding a push forward which approximates
conditioning of $\textit{Law}\,(\widehat{v}_{n+1},\widehat{y}_{n+1})$ on
$\widehat{y}_{n+1}={y}^{\dagger}$.
$\nu=\textit{Law}\,(\widehat{v}_{n+1},\widehat{y}_{n+1})$ and then identify
$\widetilde{T}^{S}_{n}:\mathbb{R}^{d_{v}}\times\mathbb{R}^{d_{y}}\to\mathbb{R}^{d_{v}}$,
where
$\widetilde{T}^{S}_{n}(\cdot)=\widetilde{T}^{S}(\cdot;\nu,{y}^{\dagger})$
approximates an exact transport map
$T^{S}_{n}(\cdot)=T^{S}(\cdot;\nu,{y}^{\dagger})$ by matching the first and
second order moments. To be precise we assume that the exact and approximate
transport maps satisfy, respectively,

(43a) $\displaystyle(T^{S}_{n})^{\sharp}\nu$ $\displaystyle=B_{n}(\nu),$ (43b)
$\displaystyle G\bigl{(}(\widetilde{T}^{S}_{n})^{\sharp}\nu\bigr{)}$
$\displaystyle=B_{n}(G\nu),$

for all measures $\nu$ on the product space
$\mathbb{R}^{d_{v}}\times\mathbb{R}^{d_{y}}.$ Of course $T^{S}_{n}$ and
$\widetilde{T}^{S}_{n}$ will depend on $\nu$ and then pushforward is to be
interpreted as in (35), (36). We ask that for all measures $\nu$ on the space
$\mathbb{R}^{d_{v}}\times\mathbb{R}^{d_{y}}$ either (43a) holds (perfect
transport) or (43b) holds (second order transport). Whilst achieving (43a) may
be harder than solving the filtering problem directly, achieving (43b) is
straightforward and leads to computationally tractable methods. This gain in
tractability comes at the price of only achieving (43b), rather than (43a);
however it is intuitive that this price will not be high for settings in which
the filtering distribution, and the predictive distribution on state and data,
is not too far from Gaussian; we flesh out this idea in Subsection 2.5.5
below.

To find tractable approximate transport we seek $\widetilde{T}^{S}$ in the
form

(44)
$\widetilde{T}^{S}(\widehat{v}_{n+1},\widehat{y}_{n+1};\nu,{y}^{\dagger}):=A\widehat{v}_{n+1}+B\widehat{y}_{n+1}+a.$

We allow the matrices/vectors $A,B,a$ to depend on $(\nu,{y}^{\dagger})$;
however, they are assumed to be independent of
$(\widehat{v}_{n+1},\widehat{y}_{n+1})$. Making this assumption ensures that
the transport map is affine with respect to the realization of
$(\widehat{v}_{n+1},\widehat{y}_{n+1})$ (but not their law). This in turn
leads to tractable computations to determine $A,B,a$ on the basis of matching
second order moments of perfect transport. In addition to computational
tractability, the affine form of the transport map $\widetilde{T}^{S}$ is
motivated by the following which shows that the approximate transport is
perfect when applied to Gaussian source:

######  Lemma 2.8.

Let
$\widetilde{T}^{S}_{n}:=\widetilde{T}^{S}(\widehat{v}_{n+1},\widehat{y}_{n+1};\nu,{y}^{\dagger})$
have the form (44) and satisfy (43b). Then $\widetilde{T}^{S}$ depends on
$\nu$ only through $G\nu$. Furthermore,

$(\widetilde{T}^{S}_{n})^{\sharp}(G\nu)=B_{n}(G\nu);$

thus, if $\nu$ is Gaussian, equation (43b) implies (43a).

######  Proof 2.9.

We first note that (43b) is equivalent to insisting that

(45) $G\bigl{(}(\widetilde{T}^{S}_{n})^{\sharp}G\nu\bigr{)}=B_{n}(G\nu)$

for all measures $\nu;$ this follows because, noting the definition (35) and
consequence (36), the first and second moments of
$(\widetilde{T}^{S}_{n})^{\sharp}G\nu$ and
$(\widetilde{T}^{S}_{n})^{\sharp}\nu$ agree, because of the affine form (44)
assumed for $\widetilde{T}^{S}_{n}$. From the identity (45), it is clear that
$\widetilde{T}^{S}_{n}$ only depends on $\nu$ through $G\nu$ because changing
$\nu\to G\nu$ leaves the identity invariant, as $G\circ G=G.$ Now note that,
because Gaussians are preserved under affine transformations,

$\bigl{(}\widetilde{T}^{S}(\widehat{v}_{n+1},\widehat{y}_{n+1};\nu,{y}^{\dagger})\bigr{)}^{\sharp}(G\nu)=G\Bigl{(}\bigl{(}\widetilde{T}^{S}(\widehat{v}_{n+1},\widehat{y}_{n+1};\nu,{y}^{\dagger})\bigr{)}^{\sharp}\nu\Bigr{)},$

or, in compact notational form,

$(\widetilde{T}^{S}_{n})^{\sharp}(G\nu)=G\bigl{(}(\widetilde{T}^{S}_{n})^{\sharp}\nu\bigr{)}.$

The desired display in the lemma is then immediate from (43b).

An affine transport map of the form (44), when combined with particle
approximations, leads to practical implementable algorithms and achieves (43b)
by ensuring that $(\widetilde{T}^{S}_{n})^{\sharp}\nu$ has first and second
moments which agree with those of the Gaussian projected filter, given by
equations (23), (24) and (30) when $\nu$ is the law of
$(\widehat{v}_{n+1},\widehat{y}_{n+1})$.

In Appendix C, Subsection 8.1, we identify the (uncountable) set of all
possible $A,B,a$ which achieve the desired matching of first and second order
moments. Here we focus on the two specific choices given in Example 8.6 from
that Appendix. The first that we highlight corresponds to the choice

$\widetilde{T}^{S}(\widehat{v}_{n+1},\widehat{y}_{n+1};\widetilde{\nu}_{n+1},{y}^{\dagger}_{n+1}):=\widehat{v}_{n+1}+K_{n}({y}^{\dagger}_{n+1}-\widehat{y}_{n+1}),$

with $K_{n}=K_{n}(\nu_{n+1})$ given by (25). Thus we obtain the following
mean-field dynamical system, which corresponds to (14) with the Kalman gain
$K_{n}$ defined by (25):

(46a) $\displaystyle\widehat{v}_{n+1}$ $\displaystyle=\Psi(v_{n})+\xi_{n},$
(46b) $\displaystyle\widehat{y}_{n+1}$
$\displaystyle=h(\widehat{v}_{n+1})+\eta_{n+1},$ (46c) $\displaystyle v_{n+1}$
$\displaystyle=\widehat{v}_{n+1}+\widehat{C}^{vy}_{n+1}(\widehat{C}^{yy}_{n+1})^{-1}({y}^{\dagger}_{n+1}-\widehat{y}_{n+1}),$

where $\\{{y}^{\dagger}_{n}\\}$ arises from a fixed realization of (1). We
refer to this as Kalman transport, noting that it serves as a derivation of
the Kalman gain, beyond the linear Gaussian setting.

The second transport map from Example 8.6 corresponds to the choice

$\widetilde{T}^{S}(\widehat{v}_{n+1},\widehat{y}_{n+1};\nu_{n+1},{y}^{\dagger}_{n+1}):=m_{n+1}+C_{n+1}^{\frac{1}{2}}\widehat{C}_{n+1}^{-\frac{1}{2}}(\widehat{v}_{n+1}-\mathbb{E}\widehat{v}_{n+1}),$

leading to the mean field map (41). A key difference between the mean field
models (46) and (41) is that the former involves inversion of matrices in data
space, and the latter in state space. The relative dimensions of the two
spaces will determine which mean field model is more appropriate as the basis
of algorithms. Note also that the mean field model (41) does not require
generation of the stochastic data $\widehat{y}_{n+1}$ because we may employ
the identity $\mathbb{E}\widehat{y}_{n+1}=\mathbb{E}h(\widehat{v}_{n+1})$ and
use (26), (27) to compute $m_{n+1},C_{n+1}$ and $\widehat{C}_{n+1}$. Motivated
by this observation, the next subsection studies a wide class of approximate
transport maps with the property that they do not require generation of
stochastic data.

####  2.5.4 Second Order Transport – Deterministic Case

We now turn our attention to approximate filters determined by deterministic
mean field maps. We define, for $\widehat{v}_{n+1}$ given by (21b),
$\widehat{\mu}=\textit{Law}\,(\widehat{v}_{n+1})$, as in the previous
subsection dropping explicit $n$-dependence on the measure $\widehat{\mu}$ and
on the data ${y}^{\dagger}$. We seek to approximate the exact transport (39)
by mean field maps in the form

(47a) $\displaystyle\widehat{v}_{n+1}$ $\displaystyle=\Psi(v_{n})+\xi_{n},$
(47b) $\displaystyle v_{n+1}$
$\displaystyle=\widetilde{T}^{D}(\widehat{v}_{n+1};\widehat{\mu}_{n+1},{y}^{\dagger}_{n+1}),$

where $\\{{y}^{\dagger}_{n}\\}$ arises from a fixed realization of (1).
Analogously to the previous subsection, we seek an approximation
$\widetilde{T}^{D}$ which, in comparison with the true transport map $T^{D}$,
satisfies

(48a) $\displaystyle(T^{D}_{n})^{\sharp}\widehat{\mu}$
$\displaystyle=B_{n}(Q\widehat{\mu}),$ (48b) $\displaystyle
G\bigl{(}(\widetilde{T}^{D}_{n})^{\sharp}\widehat{\mu}\bigr{)}$
$\displaystyle=B_{n}(GQ\widehat{\mu}),$

for all measures $\widehat{\mu}$ on the state space $\mathbb{R}^{d_{v}}.$ Here
$\widetilde{T}^{D}_{n}:\mathbb{R}^{d_{v}}\to\mathbb{R}^{d_{v}}$ is defined by
$\widetilde{T}^{D}_{n}(\cdot)=\widetilde{T}^{D}(\cdot;\widehat{\mu},{y}^{\dagger}).$
and similarly for $T^{D}_{n}.$ As in the previous subsection, where we studied
approximate stochastic transport, we again seek maps with a specific affine
form. Concretely, the maps are assumed to be affine in the pair
$\bigl{(}\widehat{v}_{n+1},\widehat{h}_{n+1}\bigr{)}$, with
$\widehat{h}_{n+1}=h(\widehat{v}_{n+1})$, leading to the assumed form
$\widetilde{T}^{D}_{n}(\cdot)=\widetilde{T}^{D}(\cdot;\mu,{y}^{\dagger})$ with

$\widetilde{T}^{D}(\widehat{v}_{n+1},\widehat{h}_{n+1};\mu,{y}^{\dagger}):=R\widehat{v}_{n+1}+S\widehat{h}_{n+1}+r,$

for $(\widehat{\mu},{y}^{\dagger})-$dependent matrices/vectors $R,S,r$ of
appropriate dimensions. Note, however, that $R,S,r$ are assumed to be
independent of the realization $(\widehat{v}_{n+1},\widehat{h}_{n+1}$),
depending only on its law, so that the transport map is affine in
$(\widehat{v}_{n+1},\widehat{h}_{n+1})$. With this restriction, which will
lead to practical implementable algorithms, we simply ask that (48b) holds:
the first and second moments of the output map agree with those of the
Gaussian projected filter, given by equations (23), (26) and (31).

As in the previous subsection, there are uncountably many choices of $R,S,r$
which we identify in Appendix C, Subsection 8.2; Example 8.15 highlights two
important cases. The first coincides with (41) since $S=0$, but the second
leads to a new mean field map. To formulate this new map we first define
$\widetilde{K}_{n}=\widetilde{K}_{n}(\widehat{\mu})$ by

(49)
$\widetilde{K}_{n}=\widehat{C}^{vh}_{n+1}\Bigl{(}(\widehat{C}^{hh}_{n+1}+\Gamma)+\Gamma^{1/2}(\widehat{C}^{hh}_{n+1}+\Gamma)^{1/2}\Bigr{)}^{-1}.$

We then make the choice

$\widetilde{T}^{D}(\widehat{v}_{n+1},\widehat{h}_{n+1};\widehat{\mu},{y}^{\dagger}):=\widehat{v}_{n+1}-\widetilde{K}_{n}(\widehat{h}_{n+1}-\mathbb{E}\widehat{h}_{n+1})+K_{n}({y}^{\dagger}-\mathbb{E}\widehat{h}_{n+1}),$

with $K_{n}$ given by (27) and $\widetilde{K}_{n}$ by (49). This leads to the
mean field map

(50a) $\displaystyle\widehat{v}_{n+1}$ $\displaystyle=\Psi(v_{n})+\xi_{n},$
(50b) $\displaystyle\widehat{h}_{n+1}$ $\displaystyle=h(\widehat{v}_{n+1}),$
(50c) $\displaystyle v_{n+1}$
$\displaystyle=\widehat{v}_{n+1}-\widetilde{K}_{n}(\widehat{h}_{n+1}-\mathbb{E}\widehat{h}_{n+1})+K_{n}({y}^{\dagger}_{n+1}-\mathbb{E}\widehat{h}_{n+1}),$

where $K_{n}$ and $\widetilde{K}_{n}$ are computed under
$\textit{Law}\,(\widehat{v}_{n+1})$.

######  Remark 2.10.

If the ensemble spread is such that the size of $\widehat{C}_{n+1}^{hh}$ is
much smaller than the size of the observational covariance $\Gamma$ then we
may invoke the approximation $\widehat{C}_{n+1}^{hh}+\Gamma\approx\Gamma$.
With this approximation it follows that
$\widetilde{K}_{n}\approx\frac{1}{2}K_{n}$ in (49). Some deterministic
ensemble Kalman filters are derived from mean-field dynamics which exploit
this approximation by setting $\widetilde{K}_{n}=\frac{1}{2}K_{n}$ in (50). We
then replace (50c) by the compact updating step

$\displaystyle
v_{n+1}=\widehat{v}_{n+1}+K_{n}\left({y}^{\dagger}_{n+1}-\frac{1}{2}\bigl{(}\mathbb{E}\widehat{h}_{n+1}+\widehat{h}_{n+1})\bigr{)}\right).$

Such a formulation corresponds to the control-theoretic perspective of (14),
with $K_{n}$ given by (27) and the innovation (15) replaced by

(51)
$I_{n}={y}^{\dagger}_{n+1}-\frac{1}{2}\bigl{(}\mathbb{E}\widehat{h}_{n+1}+\widehat{h}_{n+1}\bigr{)}.$

Filters based on this mean-field dynamics thus invoke an additional
approximation of perfect transport, over and above that stemming from matching
only first and second moments, assume further that the observational noise
dominates ensemble variation. However we will see that, in the continuous time
limit described in Section 3, this form of the innovation arises naturally and
does not constitute an additional approximation. $\hfill\blacksquare$

####  2.5.5 Second Order Transport – Summary

It is helpful at this point to take stock of two approximations to filtering
that we have introduced, Gaussian projected filtering and approximate
transport, and discuss their inter-relations. For simplicity we do this in the
context of mean field stochastic maps, but similar considerations extend to
mean field deterministic maps.

With this goal in mind, let $\mu^{MF}$ denote the measure associated with
using the mean-field map $\widetilde{T}^{S}_{n}$ to approximate the
conditioning step in (17). Notice that

(52a) $\displaystyle\mu_{n+1}$
$\displaystyle=B_{n}(QP\mu_{n}),\quad\mu_{0}=\mathsf{N}(m_{0},C_{0}),$ (52b)
$\displaystyle\mu^{G}_{n+1}$
$\displaystyle=B_{n}(GQP\mu^{G}_{n}),\quad\mu^{G}_{0}=\mathsf{N}(m_{0},C_{0}),$

and that, by (43b),

(53a) $\displaystyle\mu^{G}_{n+1}$
$\displaystyle=G\bigl{(}(\widetilde{T}^{S}_{n})^{\sharp}(QP\mu^{G}_{n})\bigr{)},\quad\mu^{G}_{0}=\mathsf{N}(m_{0},C_{0}),$
(53b) $\displaystyle\mu^{MF}_{n+1}$
$\displaystyle=(\widetilde{T}^{S}_{n})^{\sharp}(QP\mu^{MF}_{n}),\quad\mu^{MF}_{0}=\mathsf{N}(m_{0},C_{0}).$

These maps show that $\\{\mu^{MF}_{n}\\}$ is close to $\\{\mu^{G}_{n}\\}$, if
the Gaussian projection in (53a) is close to the identity where it acts, and
that $\\{\mu^{G}_{n}\\}$ is close to $\\{\mu_{n}\\}$ if the Gaussian
projection in (52b) is close to the identity where it acts. Together these two
facts suggest that $\\{\mu^{MF}_{n}\\}$, $\\{\mu^{G}_{n}\\}$ and
$\\{\mu_{n}\\}$ are all close to one another if the two Gaussian projections
can be viewed as being close to the identity map, where they appear in (52),
(53). This provides a potential path for analysis of the mean field model,
away from the linear setting where it is exact.

####  2.5.6 Mean Field Formulation of the Kalman Filter

Application of the Kalman transport map (46) in the linear Gaussian setting
leads to the following:

######  Theorem 2.11.

Assume that $v_{0}\sim\mathsf{N}(m_{0},C_{0}),$ $\Psi(\cdot)=M\cdot$ and
$h(\cdot)=H\cdot$ and that $\Gamma\succ 0.$ The mean field stochastic
dynamical system (46) takes the form

(54a) $\displaystyle\widehat{v}_{n+1}$ $\displaystyle=Mv_{n}+\xi_{n},$ (54b)
$\displaystyle\widehat{y}_{n+1}$
$\displaystyle=H\widehat{v}_{n+1}+\eta_{n+1},$ (54c) $\displaystyle v_{n+1}$
$\displaystyle=\widehat{v}_{n+1}+\widehat{C}_{n+1}H^{\top}(H\widehat{C}_{n+1}H^{\top}+\Gamma)^{-1}({y}^{\dagger}_{n+1}-\widehat{y}_{n+1}),$

where $\\{{y}^{\dagger}_{n}\\}$ arises from a fixed realization of (32) and
where $C_{n}$ is the covariance of $v_{n}$ and
$\widehat{C}_{n+1}=MC_{n}M^{\top}+\Sigma$ is the covariance of
$\widehat{v}_{n+1}.$ The resulting dynamics are identical in law to the Kalman
filter in that $v_{n}\sim\mathsf{N}(m_{n},C_{n})$ where $m_{n},C_{n}$ are as
given in Subsection 2.4.2.

######  Proof 2.12.

The map defined by (2.5.3) is well-defined since $\Gamma\succ 0.$ Lemma 2.8
shows that the approximate transport is exact in this Gaussian setting.

Similar ideas can be applied to (41) and (50) to determine other mean field
models with law equal to that of the Kalman filter. Furthermore we observe
that the formulation based on (41) can be symmetrized to obtain, in the linear
setting of Theorem 2.11,

(55a) $\displaystyle\widehat{v}_{n+1}$
$\displaystyle=Mv_{n}+\xi_{n},\>n\in\mathbb{Z}^{+},$ (55b) $\displaystyle
v_{n+1}$ $\displaystyle=m_{n+1}+A_{n}(\widehat{v}_{n+1}-\widehat{m}_{n+1}),$
(55c) $\displaystyle A_{n}$
$\displaystyle=(C_{n+1})^{1/2}\left[(C_{n+1})^{1/2}\widehat{C}_{n+1}(C_{n+1})^{1/2}\right]^{-1/2}(C_{n+1})^{1/2},$

with
$(\widehat{m}_{n+1},\widehat{C}_{n+1},\widehat{C}_{n+1}^{vh},\widehat{C}_{n+1}^{hh})$
defined by (2.4.1) and $(m_{n+1},C_{n+1})$ defined by (30). We note that the
second component of the map may be written in gradient form and corresponds to
an optimal transport from $\widehat{\mu}_{n+1}$ into $\mu_{n+1}$ in the sense
of the Euclidean Wasserstein distance of optimal transportation (see
Subsection 2.7 for details); indeed this is true of for an entire family of
weighted Wasserstein distances – see Example 8.9. To recognize the gradient
structure define

$\Phi_{n}(v):=\langle m_{n+1},v\rangle+\frac{1}{2}\langle
A_{n}(v-\widehat{m}_{n+1}),v-\widehat{m}_{n+1}\rangle$

and note that then

(56a) $\displaystyle\widehat{v}_{n+1}$
$\displaystyle=Mv_{n}+\xi_{n},\>n\in\mathbb{Z}^{+},$ (56b) $\displaystyle
v_{n+1}$ $\displaystyle=\nabla\Phi_{n}(\widehat{v}_{n+1}).$

###  2.6 Ensemble Kalman Methods

In this subsection we take the mean field models from Subsection 2.5 and use
particle approximations to derive implementable numerical algorithms. We
start, in Subsection 2.6.1, with the setting in which transport is perfect;
these algorithms are not, in general, implementable since determining perfect
transport is itself a difficult computational task and the subject of ongoing
research; see the bibliography Subsection 2.7. Thus we turn to particle
approximations of the transports designed to match first and second order
statistics. This leads to the (stochastic) ensemble Kalman filter in
Subsection 2.6.2, and to (deterministic) ensemble square root filters in
Subsection 2.6.3. The methods derived in this subsection involve approximating
the Kalman gain by computing covariances under the empirical measure defined
by the ensemble of particles. To avoid overloading notation, in the rest of
this subsection $C_{n+1}$, $\widehat{C}_{n+1}$, $\widehat{C}_{n+1}^{vy}$ and
$\widehat{C}_{n+1}^{yy}$ will denote covariances computed with expectation
under the empirical measure. With this notation in place, in the rest of this
subsection $K_{n}$ will directly refer to the particle approximation of the
Kalman gain, computed using the covariances with respect to the empirical
measure, without further specification needed. Throughout this section
$\\{{y}^{\dagger}_{n}\\}$ arises from a fixed realization of (1).

####  2.6.1 Perfect Particle Filters

The mean-field equations (37) could, in principle, be approximated through a
particle approximation of the mean-field leading to the following conceptual
(because map $T^{S}_{n}$ is not known explicitly) algorithm: let
$\mathsf{J}=\\{1,\cdots,J\\}$ and consider, for
$(n,j)\in\mathbb{Z}^{+}\times\mathsf{J}$, the interacting particle dynamical
system

(57a) $\displaystyle\widehat{v}_{n+1}^{(j)}$
$\displaystyle=\Psi(v_{n}^{(j})+\xi_{n}^{(j)},$ (57b)
$\displaystyle\widehat{y}_{n+1}^{(j)}$
$\displaystyle=h(\widehat{v}_{n+1}^{(j})+\eta_{n+1}^{(j)},$ (57c)
$\displaystyle v_{n+1}^{(j)}$
$\displaystyle=T^{S}(\widehat{v}_{n+1}^{(j)},\widehat{y}_{n+1}^{(j)};\nu_{n+1}^{\mathsf{J}},{y}^{\dagger}_{n+1}),$
(57d) $\displaystyle\nu_{n+1}^{\mathsf{J}}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\delta_{(\widehat{v}_{n+1}^{(j)},\widehat{y}_{n+1}^{(j)})}.$

This evolves the particles $\\{{v}^{(j)}_{n}\\}_{j\in\mathsf{J}}$ into
$\\{{v}^{(j)}_{n+1}\\}_{j\in\mathsf{J}}$. Here the $\\{{\xi}^{(j)}_{n}\\}$
are, for each $j$, random variables given by the known distribution of
$\xi_{n}$ specified in (2) and, furthermore, are drawn independently with
respect to each $(n,j)\in\mathbb{Z}^{+}\times\mathsf{J}$. Similar
considerations apply to the $\\{{\eta}^{(j)}_{n}\\}$ which, additionally, are
independent of the $\\{{\xi}^{(j)}_{n}\\}$. It is intuitive that the large $J$
limit of this system recovers the mean-field dynamics (37) and, in particular,

(58)
$\mu_{n}^{\mathsf{J}}=\frac{1}{J}\sum_{j=1}^{J}\delta_{v_{n}^{(j)}}\approx\mu_{n}.$

Applying a similar idea to (39) leads to the following conceptual (because map
$T^{D}_{n}$ is not known explicitly) algorithm. Consider, for
$(n,j)\in\mathbb{Z}^{+}\times\mathsf{J}$, the interacting particle dynamical
system

(59a) $\displaystyle\widehat{v}_{n+1}^{(j)}$
$\displaystyle=\Psi(v_{n}^{(j})+\xi_{n}^{(j)},$ (59b) $\displaystyle
v_{n+1}^{(j)}$
$\displaystyle=T^{D}(\widehat{v}_{n+1}^{(j)};\widehat{\mu}_{n+1}^{\mathsf{J}},{y}^{\dagger}_{n+1}),$
(59c) $\displaystyle\widehat{\mu}_{n+1}^{\mathsf{J}}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\delta_{\widehat{v}_{n+1}^{(j)}}.$

This evolves the particles $\\{{v}^{(j)}_{n}\\}_{j\in\mathsf{J}}$ into
$\\{{v}^{(j)}_{n+1}\\}_{j\in\mathsf{J}}$. The same assumptions are made about
the $\\{{\xi}^{(j)}_{n}\\}$ as for the preceding interacting particle
dynamical system. It is again intuitive that the large $J$ limit of this
system recovers the mean field dynamics (39), and the evolution (40). In
particular it is intuitive that (58) holds for this particle approximation
too. We reiterate that in practice these algorithms are, in general, not easy
to use. More specifically, finding particle based approximations $T^{S}_{J}$
and $T^{D}_{J}$ to the desired transport maps such that

$\lim_{J\to\infty}T^{S}_{J}=T^{S},\qquad\lim_{J\to\infty}T^{D}_{J}=T^{D}$

in an appropriate sense is a computationally challenging task and the subject
of ongoing research. This leads to the next two subsections in which we
replace $T^{S}$ and $T^{D}$, in the interacting particles systems (57) and
(59), by the previously introduced approximate transports $\widetilde{T}^{S}$
and $\widetilde{T}^{D},$ respectively.

####  2.6.2 Stochastic Ensemble Kalman Filters

Particle approximation of the mean field dynamical system (46), effecting
Kalman transport, bring us to the stochastic EnKF (ensemble Kalman filter).
This method may be derived by writing down a particle approximation of the
mean field stochastic dynamics defined by (46). We evolve the particles
$\\{{v}^{(j)}_{n}\\}_{j\in\mathsf{J}}$ into
$\\{{v}^{(j)}_{n+1}\\}_{j\in\mathsf{J}}$ according to the following stochastic
interacting particle system, holding for
$(n,j)\in\mathbb{Z}^{+}\times\mathsf{J}$:

(60a) $\displaystyle\widehat{v}^{(j)}_{n+1}$
$\displaystyle=\Psi({v}^{(j)}_{n})+{\xi}^{(j)}_{n},\>n\in\mathbb{Z}^{+},$
(60b) $\displaystyle\widehat{y}^{(j)}_{n+1}$
$\displaystyle=h(\widehat{v}^{(j)}_{n+1})+{\eta}^{(j)}_{n+1},\>n\in\mathbb{Z}^{+},$
(60c) $\displaystyle{v}^{(j)}_{n+1}$
$\displaystyle=\widehat{v}^{(j)}_{n+1}+K_{n}\bigl{(}{y}^{\dagger}_{n+1}-\widehat{y}^{(j)}_{n+1}\bigr{)},$
(60d) $\displaystyle\nu_{n+1}^{\mathsf{J}}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\delta_{(\widehat{v}_{n+1}^{(j)},\widehat{y}_{n+1}^{(j)})}.$

Here the Kalman gain from (25) is approximated using the empirical measure
$\nu_{n+1}^{\mathsf{J}}$ and but is still denoted by $K_{n}$; details follow
below. The same assumptions regarding $\\{{\xi}^{(j)}_{n}\\}$ and
$\\{{\eta}^{(j)}_{n+1}\\}$ are made as for (57). We let
$\mathbb{E}^{\mathsf{J}}_{n}$ denote expectation under $\nu_{n}^{\mathsf{J}}.$
For the basic implementation of EnKF (60) the desired covariance matrices, and
Kalman gain (25) are then approximated by expectation under
$\nu_{n+1}^{\mathsf{J}}$, so that 777The empirical covariance computations are
often modified to accommodate the widely adopted convention of scaling by
$1/(J-1)$, instead of $1/J$, in view of the matrix being computed from $J-1$
independent increments about the mean.

$\displaystyle\widehat{C}_{n+1}^{vy}=$
$\displaystyle\,\mathbb{E}^{\mathsf{J}}_{n+1}\Bigl{(}\bigl{(}\widehat{v}_{n+1}-\mathbb{E}^{\mathsf{J}}_{n+1}\widehat{v}_{n+1}\bigr{)}\otimes\bigl{(}\widehat{y}_{n+1}-\mathbb{E}^{\mathsf{J}}_{n+1}\widehat{y}_{n+1}\bigr{)}\Bigr{)},$
$\displaystyle\widehat{C}_{n+1}^{yy}=$
$\displaystyle\,\mathbb{E}^{\mathsf{J}}_{n+1}\Bigl{(}\bigl{(}\widehat{y}_{n+1}-\mathbb{E}^{\mathsf{J}}_{n+1}\widehat{y}_{n+1}\bigr{)}\otimes\bigl{(}\widehat{y}_{n+1}-\mathbb{E}^{\mathsf{J}}_{n+1}\widehat{y}_{n+1}\bigr{)}\Bigr{)}.$

An alternative is to use a particle approximation in formula (27b), leading to
the use of

$\displaystyle\widehat{C}_{n+1}^{vh}=$
$\displaystyle\,\mathbb{E}_{n+1}^{\mathsf{J}}\Bigl{(}\bigl{(}\widehat{v}_{n+1}-\mathbb{E}^{\mathsf{J}}_{n+1}\widehat{v}_{n+1}\bigr{)}\otimes\bigl{(}\widehat{y}_{n+1}-\mathbb{E}_{n+1}^{\mathsf{J}}h(\widehat{v}_{n+1})\bigr{)}\Bigr{)},$
$\displaystyle\widehat{C}_{n+1}^{hh}=$
$\displaystyle\,\mathbb{E}_{n+1}^{\mathsf{J}}\Bigl{(}\bigl{(}\widehat{y}_{n+1}-\mathbb{E}^{\mathsf{J}}_{n+1}h(\widehat{v}_{n+1})\bigr{)}\otimes\bigl{(}\widehat{y}_{n+1}-\mathbb{E}_{n+1}^{\mathsf{J}}h(\widehat{v}_{n+1})\bigr{)}\Bigr{)}+\Gamma$

in computation of the gain $K_{n}.$ The advantage of this latter formulation
is that it ensures positivity, and hence invertibility, of the covariance in
data space, if $\Gamma$ is assumed positive-definite.

Pseudo-code for the stochastic EnKF may be found as Algorithm 2 in Appendix 6.

######  Example 2.13.

We return to the set-up of Example 2.1, and now demonstrate performance of the
stochastic EnKF on the same Lorenz ’96 model. Indeed, we again study the
Lorenz ’96 (singlescale) model for unknown $v\in
C(\mathbb{R}^{+},\mathbb{R}^{L})$ satisfying the equations (6) with
$L=9,h_{v}=-0.8$ and $F=10.$ We consider observations
$\\{{y}^{\dagger}_{n}\\}_{n\in\mathbb{Z}^{+}}$ arising from the model

$\displaystyle{v}^{\dagger}_{n+1}$
$\displaystyle=\Psi({v}^{\dagger}_{n})+{\xi}^{\dagger}_{n},$
$\displaystyle{y}^{\dagger}_{n+1}$
$\displaystyle=h({v}^{\dagger}_{n+1})+{\eta}^{\dagger}_{n+1},$

where $\Psi$ is the solution operator for (6) over the observation time
interval $\tau$, and $\\{{\xi}^{\dagger}_{n}\\}_{n\in\mathbb{Z}^{+}}$,
$\\{{\eta}^{\dagger}_{n}\\}_{n\in\mathbb{N}}$ are mutually independent
Gaussian sequences defined by

${\xi}^{\dagger}_{n}\sim\mathsf{N}(0,\sigma^{2}I)\,\,\text{i.i.d.},\quad{\eta}^{\dagger}_{n}\sim\mathsf{N}(0,\gamma^{2}I)\,\,\text{i.i.d.}\,\,,$

with $\sigma^{2}=0.1$ and $\gamma^{2}=0.1$. We again assume that the
observation function is linear: $h(v)=Hv$ for matrix
$H:\mathbb{R}^{9}\to\mathbb{R}^{6}$ defined by (7).

Figures 8(a) and 8(b) demonstrate the performance of stochastic EnKF in this
experimental setting with $\tau=10^{-3}$ and using $J=10^{2}$ and $J=10^{3}$,
respectively, against the performance of 3DVAR with no noise; note that the
EnKF uses a time-varying estimate of the gain $K_{n}$, whilst 3DVAR uses the
fixed $K$ given in Example 2.1. These experiments illustrate that using
sufficiently large ensembles, the ensemble Kalman filter outperforms 3DVAR on
such a non-linear filtering problem where the true state and observational
noise levels are high. To quantitatively demonstrate this improvement, we
compute the mean squared error between the estimates yielded by 3DVAR and the
true states and the mean squared error between the estimates obtained by
computing the ensemble mean of the stochastic EnKF and comparing with the true
states. In particular we report time-averaged mean squared errors obtained
from both 3DVAR and stochastic EnKF given by use of formula (16) from Example
2.3 using $t=5$ and $T=20$. An ensemble size of $J=10^{2}$ yields
$e_{\text{EnKF}}=8.20\cdot 10^{-1}$, while for $J=10^{3}$ we obtain
$e_{\text{EnKF}}=4.91\cdot 10^{-1}$. For comparison, the error obtained using
3DVAR is $e_{\text{3DVAR}}=1.45\cdot 10^{0}$.

(a) EnKF with ensemble size $J=10^{2}$.

(b) EnKF with ensemble size $J=10^{3}$.

Figure 8: In this experiment we set the noise levels
$\sigma^{2}=10^{-1},\gamma^{2}=10^{-1}$. We display the estimates of $v_{3}$
in time produced by EnKF and 3DVAR against the true dynamics using
$\tau=10^{-3}$. Again “dns” refers to direct numerical simulation.

####  2.6.3 Ensemble Square Root Filters

The variants on the EnKF described in this subsection are known as ensemble
square root filters; they are based on mean field maps (41) and (50),
approximated by interacting particle systems. They are also sometimes referred
to as deterministic ensemble Kalman filters, to distinguish them from their
stochastic counterparts described in the preceding subsection. We introduce
two families of square root filters: of adjustment and transform type. We
emphasize again that the choice of which method to use in practice is
determined by implementations details such as the number of particles $J$, the
dimension of state space $d_{v}$, and the dimension of the observation space
$d_{y}$.

##### Ensemble Adjustment Kalman Filters

As in the previous subsection the methods evolve particle ensemble
$\\{{v}^{(j)}_{n}\\}_{j\in\mathsf{J}}$ into
$\\{{v}^{(j)}_{n+1}\\}_{j\in\mathsf{J}}$, via the predictive ensemble
$\\{\widehat{v}^{(j)}_{n+1}\\}_{j\in\mathsf{J}}.$ However they do not employ
simulated data $\\{\widehat{y}^{(j)}_{n+1}\\}_{j\in\mathsf{J}}$, rather they
make use of $\\{\widehat{h}^{(j)}_{n+1}\\}_{j\in\mathsf{J}}$, where
$\widehat{h}^{(j)}_{n}=h(\widehat{v}^{(j)}_{n}).$ To define the methods it
helps to introduce new notation. Slightly modifying the notation in the
preceding subsection, we now let $\mathbb{E}^{\mathsf{J}}_{n}$ denote
expectation with respect to the empirical measure

$\widehat{\mu}_{n}^{\mathsf{J}}=\frac{1}{J}\sum_{j=1}^{J}\delta_{\widehat{v}^{(j)}_{n}}.$

We let $\widehat{v}_{n+1}$ denote the random variable with this distribution
and, as before, $\widehat{h}_{n+1}=h(\widehat{v}_{n+1}).$

Next we define matrix $\widehat{V}_{n}$ comprising scaled ensemble deviations
in state space:

$\widehat{V}_{n}=\frac{1}{\sqrt{J}}\left(\widehat{v}_{n}^{(1)}-\mathbb{E}^{\mathsf{J}}_{n}\widehat{v}_{n},\widehat{v}_{n}^{(2)}-\mathbb{E}^{\mathsf{J}}_{n}\widehat{v}_{n},\cdots,\widehat{v}_{n}^{(J)}-\mathbb{E}^{\mathsf{J}}_{n}\widehat{v}_{n}\right)\in\mathbb{R}^{d_{v}\times
J},$

we then define the analogous matrix $\widehat{H}_{n}$ in observation space

$\widehat{H}_{n}=\frac{1}{\sqrt{J}}\left(\widehat{h}_{n}^{(1)}-\mathbb{E}^{\mathsf{J}}_{n}\widehat{h}_{n},\widehat{h}_{n}^{(2)}-\mathbb{E}^{\mathsf{J}}_{n}\widehat{h}_{n},\cdots,\widehat{h}_{n}^{(J)}-\mathbb{E}^{\mathsf{J}}_{n}\widehat{h}_{n}\right)\in\mathbb{R}^{d_{y}\times
J}.$

For later use we also define

$V_{n}=\frac{1}{\sqrt{J}}\left(v_{n}^{(1)}-\mathbb{E}_{*,n}^{\mathsf{J}}v_{n},v_{n}^{(2)}-\mathbb{E}_{*,n}^{\mathsf{J}}v_{n},\cdots,v_{n}^{(J)}-\mathbb{E}_{*,n}^{\mathsf{J}}v_{n}\right)\in\mathbb{R}^{d_{v}\times
J},$

where, here, expectation $\mathbb{E}_{*,n}^{\mathsf{J}}$ is with respect to
the empirical measure (58) and $v_{n}$ is a random variable with this
distribution.

With these notations in hand, we have, with expectations computed under
$\mathbb{E}^{\mathsf{J}}_{n}$,

(61a) $\displaystyle\widehat{C}^{vh}_{n}$
$\displaystyle=\widehat{V}_{n}\widehat{H}_{n}^{\top},$ (61b)
$\displaystyle\widehat{C}^{hh}_{n}$
$\displaystyle=\widehat{H}_{n}\widehat{H}_{n}^{\top}$

and the ensemble-based approximation of the Kalman gain matrix is 888Here too,
the empirical covariance computations are often modified to accommodate the
widely adopted convention of scaling by $1/(J-1)$, instead of $1/J$.

(62)
$K_{n}=\widehat{V}_{n+1}\widehat{H}_{n+1}^{\top}(\widehat{H}_{n+1}\widehat{H}_{n+1}^{\top}+\Gamma)^{-1}.$

This is a linear algebraic reformulation of the Kalman gain approximation
resulting from (2.6.2). By making a particle approximation of (41) we obtain

(63a) $\displaystyle\widehat{v}^{(j)}_{n+1}$
$\displaystyle=\Psi({v}^{(j)}_{n})+{\xi}^{(j)}_{n},\>n\in\mathbb{Z}^{+},$
(63b) $\displaystyle\widehat{m}_{n+1}$
$\displaystyle=\mathbb{E}^{\mathsf{J}}_{n+1}\widehat{v}_{n+1},$ (63c)
$\displaystyle m_{n+1}$
$\displaystyle=\widehat{m}_{n+1}+K_{n}({y}^{\dagger}_{n+1}-\mathbb{E}^{\mathsf{J}}_{n+1}\widehat{h}_{n+1}),$
(63d) $\displaystyle v_{n+1}^{(j)}$
$\displaystyle=m_{n+1}+C_{n+1}^{\frac{1}{2}}\widehat{C}_{n+1}^{-\frac{1}{2}}\left(\widehat{v}_{n+1}^{(j)}-\widehat{m}_{n+1}\right),$

where $K_{n}$ is given by (62), whilst $\widehat{C}_{n+1}$ and $C_{n+1}$ are
computed empirically using

(64) $\widehat{C}_{n+1}=\widehat{V}_{n+1}\widehat{V}_{n+1}^{\top},\quad
C_{n+1}=\widehat{V}_{n+1}\left(I+\widehat{H}_{n+1}^{\top}\Gamma^{-1}\widehat{H}_{n+1}\right)^{-1}\widehat{V}_{n+1}^{\top}.$

The first of these two formulae follows similarly to (61); the second of these
two formulae is derived as follows: 999Using, in the last line, the identity
$I-W^{\top}(WW^{\top}+I)^{-1}W=(I+W^{\top}W)^{-1}$ which holds for all (not
necessarily square) matrices $W$.

(65a) $\displaystyle C_{n+1}$
$\displaystyle=\widehat{C}_{n+1}-\widehat{C}^{vh}_{n+1}(\widehat{C}^{hh}_{n+1}+\Gamma)^{-1}(\widehat{C}^{vh}_{n+1})^{\top}$
(65b)
$\displaystyle=\widehat{V}_{n+1}\widehat{V}_{n+1}^{\top}-\widehat{V}_{n+1}\widehat{H}_{n+1}^{\top}(\widehat{H}_{n+1}\widehat{H}_{n+1}^{\top}+\Gamma)^{-1}\widehat{H}_{n+1}\widehat{V}_{n+1}^{\top},$
(65c)
$\displaystyle=\widehat{V}_{n+1}\left(I-\widehat{H}_{n+1}^{\top}(\widehat{H}_{n+1}\widehat{H}_{n+1}^{\top}+\Gamma)^{-1}\widehat{H}_{n+1}\right)\widehat{V}_{n+1}^{\top}$
(65d)
$\displaystyle=\widehat{V}_{n+1}\left(I+\widehat{H}_{n+1}^{\top}\Gamma^{-1}\widehat{H}_{n+1}\right)^{-1}\widehat{V}_{n+1}^{\top}.$

Similarly, we may apply a particle approximation to the mean field dynamical
system (50) to obtain

(66a) $\displaystyle\widehat{v}^{(j)}_{n+1}$
$\displaystyle=\Psi({v}^{(j)}_{n})+{\xi}^{(j)}_{n},\>n\in\mathbb{Z}^{+},$
(66b) $\displaystyle\widehat{h}_{n+1}^{(j)}$
$\displaystyle=h(\widehat{v}_{n+1}^{(j)}),$ (66c)
$\displaystyle\widehat{m}_{n+1}$
$\displaystyle=\mathbb{E}^{\mathsf{J}}_{n+1}\widehat{v}_{n+1},$ (66d)
$\displaystyle m_{n+1}$
$\displaystyle=\widehat{m}_{n+1}+K_{n}({y}^{\dagger}_{n+1}-\mathbb{E}^{\mathsf{J}}_{n+1}\widehat{h}_{n+1}),$
(66e) $\displaystyle v_{n+1}^{(j)}$
$\displaystyle=m_{n+1}+(\widehat{v}_{n+1}^{(j)}-\widehat{m}_{n+1})-\widetilde{K}_{n}(\widehat{h}_{n+1}^{(j)}-\mathbb{E}^{\mathsf{J}}_{n+1}\widehat{h}_{n+1}).$

By making an empirical approximation of the formula (49) the matrix
$\widetilde{K}_{n}$ is defined using the identification

$\widetilde{K}_{n}=\widehat{V}_{n+1}\widehat{H}_{n+1}^{\top}\left[(\widehat{H}_{n+1}\widehat{H}_{n+1}^{\top}+\Gamma)+(\widehat{H}_{n+1}\widehat{H}_{n+1}^{\top}+\Gamma)^{1/2}\Gamma^{1/2}\right]^{-1}$

The key difference between (63) and (66) is that the former involves inversion
in state space, and the latter in data space. The relative size of the two
dimensions dictates which is preferable. Both methods are examples of a
general class of algorithms known as ensemble adjustment Kalman filters: EAKF.

##### Ensemble Transform Kalman Filters

The ensemble adjustment Kalman filters involve application, and inversion, of
matrices which are applied on the left and act on state space. A different
class of algorithms, known as ensemble transform Kalman filters (ETKF),
involve matrix multiplication from the right and consequently inversions take
place in the ensemble space of dimension $J$. In many applications this is far
smaller than the dimension of the state or data spaces, and then use of this
version of the methodology is preferred. The aim is to determine matrix
$Z_{n}\in\mathbb{R}^{J\times J}$, and to derive Kalman gain $K_{n}$ from
$Z_{n}$, so that the following interacting particle system produces an
ensemble of particles $\\{{v}^{(j)}_{n+1}\\}_{j\in\mathsf{J}}$ with empirical
covariance $C_{n+1}$ defined by the second item in display (64):

(67a) $\displaystyle\widehat{v}^{(j)}_{n+1}$
$\displaystyle=\Psi({v}^{(j)}_{n})+{\xi}^{(j)}_{n},\>n\in\mathbb{Z}^{+},$
(67b) $\displaystyle\widehat{m}_{n+1}$
$\displaystyle=\mathbb{E}^{\mathsf{J}}_{n+1}\widehat{v}_{n+1},$ (67c)
$\displaystyle m_{n+1}$
$\displaystyle=\widehat{m}_{n+1}+K_{n}({y}^{\dagger}_{n+1}-\mathbb{E}^{\mathsf{J}}_{n+1}\widehat{h}_{n+1}),$
(67d) $\displaystyle v_{n+1}^{(j)}$
$\displaystyle=m_{n+1}+\sum_{i=1}^{J}\left(\widehat{v}_{n+1}^{(i)}-\widehat{m}_{n+1}\right)(Z_{n})_{ij}.$

It follows from (67d) that

(68) $V_{n+1}=\widehat{V}_{n+1}Z_{n}.$

If we define

(69)
$Z_{n}=\left(I+\widehat{H}_{n+1}^{\top}\Gamma^{-1}\widehat{H}_{n+1}\right)^{-1/2}\in\mathbb{R}^{J\times
J}.$

then, as desired, $V_{n+1}V_{n+1}^{\top}=C_{n+1}$ as defined by (64), by
virtue of (65). The calculations in (65) can also be utilized to verify that
the empirical Kalman gain matrix defined by (62) satisfies

$K_{n}=\widehat{V}_{n+1}Z_{n}^{2}\widehat{H}_{n+1}^{\top}\Gamma^{-1}.$

Using this formula for $K_{n}$ in (63) leads to an algorithm which matches
first and second order statistics of the Gaussian projected filter, at $n+1$,
requiring only matrix inversions in space of dimension defined by the number
of particles $J$.

We finally note that (67c) and (67d) can be combined into a single
transformation step of the form

(70) $v_{n+1}^{(j)}=\sum_{i=1}^{J}\widehat{v}_{n+1}^{(i)}(S_{n})_{ij}$

where $S_{n}\in\mathbb{R}^{J\times J}$ replaces the matrix $Z_{n}$ in (67d)
such that

$m_{n+1}=\sum_{j=1}^{J}v_{n+1}^{(j)}=\sum_{i,j=1}^{J}\widehat{v}_{n+1}^{(i)}(S_{n})_{ij}$

holds in addition to (68) with $Z_{n}$ replaced by $S_{n}$. Formulation (70)
has a number of attractive features. First, it clearly reveals that the
analysis $\\{{v}^{(j)}_{n+1}\\}_{j\in\mathsf{J}}$ lies in the span of the
space spanned by the predictions
$\\{\widehat{v}^{(j)}_{n+1}\\}_{j\in\mathsf{J}}$, which is relevant whenever
$J<d_{v}$. Second, all particle implementations of the ensemble Kalman filter
and many of its extensions can be put into the framework (70) with the
(possibly random) matrix $S_{n}$ chosen appropriately. Third, it encodes a
coupling between the prediction
$\\{\widehat{v}^{(j)}_{n+1}\\}_{j\in\mathsf{J}}$ and the analysis
$\\{{v}^{(j)}_{n+1}\\}_{j\in\mathsf{J}}$ at the level of their associated
empirical measures $\widehat{\mu}_{n}^{\mathsf{J}}$ and
$\mu_{n}^{\mathsf{J}}$, respectively. See the bibliographic Section 2.7 below
for more details.

###  2.7 Bibliographical Notes

For the control theoretic approach to the state estimation problem see
Luenberger (1964) and Luenberger (1971). Recent analysis of the closely
related 3DVAR method rests heavily on ideas introduced in Olson & Titi (2003),
Hayden et al. (2011) and Foias et al. (2016); the papers Law & Stuart (2012),
Law et al. (2012), Law et al. (2016a) and Sanz-Alonso & Stuart (2015)
essentially establish the stability of these deterministic results to small
noise perturbations. In the context of using observations to control the
instability of chaotic systems, all of this work may be seen as building on
the study of synchronization, reviewed in Ashwin (2003).

The Kalman filter (Kalman, 1960) is perhaps the first systematic analysis of
an algorithm for incorporation of discrete time data into estimation of a
discrete time stochastic dynamical system; it applies only to linear Gaussian
systems. The monograph Jazwinski (2007) provides an introduction to nonlinear
filtering both in discrete and continuous time.

Particle-based extensions of the classical Kalman filter to nonlinear
filtering problems include the unscented Kalman filter and the ensemble Kalman
filter. While this review focuses primarily on ensemble Kalman filter
techniques, the unscented Kalman filter is an approach based on application of
quadrature to the Gaussian projected filter; see Julier et al. (2000) as well
as Särkkä (2013). A discussion and evaluation in the context of ensemble
square root filters may be found in Wang et al. (2004).

The books Reich & Cotter (2015), Asch et al. (2016), Law et al. (2015), Harlim
& Majda (2010), Abarbanel (2013), and Evensen et al. (2022) provide overviews
of a variety of filtering methods, and ensemble Kalman methods in particular.
The numerical implementation of ensemble Kalman methods is a very rich subject
on its own with many important details which are beyond the scope of this
mathematical survey paper. We refer the reader to the excellent monographs by
Asch et al. (2016) and Evensen et al. (2022) for texts with emphasis on these
implementation details. These include techniques such as inflation and
localization, which we do not cover in this mathematical survey. Here we wish
to only point to two particular aspects. First, the matrix $Z_{n}$ in (68) is
not uniquely defined by the requirement $V_{n+1}V_{n+1}^{\top}=C_{n+1}$.
Formula (69) constitutes one possible choice, which leads to a symmetric
$Z_{n}$. See Livings et al. (2008) for more details. Second, finite particle
implementations of the stochastic EnKF entail that the random realizations
$\widehat{y}^{(j)}_{n+1}$ appear both in the Kalman gain $K_{n}$ as well as in
the innovation term in (60c). As first observed by Houtekamer & Mitchell
(2005), this leads to a systematic underestimation of the ensemble spread,
which vanishes in the $J\to\infty$ limit; but can affect the performance of
the EnKF for small particle sizes. In Subsection 4.6 we will highlight the
same effect when discussing finite particle implementations (Nüsken & Reich,
2019; Garbuno-Inigo et al., 2020a) of the ensemble Kalman sampler for Bayesian
inversion, based on the mean field model proposed in Garbuno-Inigo et al.
(2020b).

Much of the development of ensemble Kalman methods reflects the historical
roots of the subject in the geophysical sciences, the atmosphere-ocean
sciences in particular, including Lagrangian data assimilation, and in the
modeling of subsurface flow (Burgers et al., 1998; Houtekamer & Mitchell,
1998; Anderson, 2001; Bishop et al., 2001; Whitaker & Hamill, 2002; Tippett et
al., 2003; Hunt et al., 2007; Li & Reynolds, 2009; Sakov et al., 2012; Bocquet
& Sakov, 2014; Evensen, 2019; Bocquet & Sakov, 2012; Bocquet et al., 2017;
Gurumoorthy et al., 2017; Sampson et al., 2021; Kuznetsov et al., 2003; Salman
et al., 2006). We also mention the randomized maximum likelihood (RML)
approach to Bayesian inference which is closely related to the analysis step
of a stochastic EnKF and has also been developed primarily through application
in the geophysical sciences Kitanidis (1995), Oliver et al. (1997), and Oliver
et al. (2008).

The following papers concern analysis and development of ensemble methods with
an emphasis on applications in complex and turbulent flows: Grooms et al.
(2014), Grooms et al. (2015), Robinson et al. (2018), Lee et al. (2017),
Gottwald & Majda (2013), Kelly et al. (2015), Tong et al. (2016), Tong et al.
(2015), Kelly et al. (2015), Harlim et al. (2014), Majda & Tong (2018), Fertig
et al. (2007), Harlim & Hunt (2007b), Harlim & Hunt (2007a). The conceptual
fluid dynamics models of Lorenz (Lorenz, 1996) (often referred to,
collectively, as Lorenz ’96 models) have been particularly influential in
germinating this body of work and we use them exclusively in our illustrative
Examples 2.1, 2.2, 2.3, 2.13 and 4.33. Furthermore we will make use of the
relationship between the multiscale and singlescale version of the model as
developed in (Fatkullin & Vanden-Eijnden, 2004).

Recently ensemble Kalman methods have been developed for potential use in
machine learning (Haber et al., 2018; Kovachki & Stuart, 2019; Guth et al.,
2020; Grooms, 2021; Gottwald & Reich, 2021; Yang & Grooms, 2021; Pidstrigach &
Reich, 2021); see also Bocquet et al. (2017), Chen et al. (2022) for research
at the intersection of machine learning with ensemble Kalman methodology.

While we have started from mean field equations in this survey and then
discretized by $J$ particles, it remains to be demonstrated that the arising
discrete formulations converge to the mean field equations in the $J\to\infty$
limit. This has indeed been establish for the standard ensemble Kalman filter
formulations by Le Gland et al. (2011), Mandel et al. (2011), and Kwiatkowski
& Mandel (2015); see Law et al. (2016b) for related work studying ensemble
Kalman methods in the context of non-Gaussian problems. The papers Ding et al.
(2020), Ding & Li (2021a, b) study particle approximation of mean-field limits
beyond the Gaussian setting, primarily in the context of the solution of
inverse problems; see the discussion in Subsection 4.6. The papers Hoel et al.
(2016); Chernov et al. (2021) study the use of multilevel approximation of the
mean field limit, coupling ensemble approximations at different levels of
space or time discretization.

The non-uniqueness of second order transport maps is studied in continuous
time, for linear Gaussian stochastic differential equations, in Taghvaei &
Mehta (2020); this work is closely related to our analysis in the first two
subsections of Section 8. Non-Gaussian extensions are discussed in Taghvaei &
Hosseini (2022). We also highlight that it is possible to construct second
order transport maps $\widetilde{T}$, which satisfy conditions different from
(43b) and (48b), respectively. For example, one could request that

$G\bigl{(}(\widetilde{T}^{S}_{n})^{\sharp}\nu\bigr{)}=GB_{n}(\nu).$

This approximation has been utilized by Lei & Bickel (2011). The analog
deterministic approach has been put forward by Tödter & Ahrens (2015) and has
been explored further, for example, in Acevedo et al. (2017). Alternatively,
one can also replace the definition (2.4.1) of the Gaussian projection
operator $G$. To this end, recall that the Kullback–Leibler divergence between
probability measures $\mu$ and $\pi$ on $\mathbb{R}^{d}$ is defined as

$d_{\rm{KL}}(\mu||\pi)=\int\pi(du)\log\frac{d\pi}{d\mu}(u);$

it is not symmetric in its two arguments. Gaussian variational inference
Bishop (2011) is for example based of the definition

(71) $G\mu={\rm argmin}_{\pi\in\mathfrak{G}}d_{\rm{KL}}(\pi||\mu)$

in place of (2.4.1); note, however, that the minimizer of (71) may not be
unique, whilst the minimizer of (2.4.1) is always unique.

While we follow the moment matching perspective on the derivation of ensemble
Kalman filter methods in this survey, we mention in passing that there is an
alternative perspective based on linear minimum variance estimators. See van
Leeuwen (2020) in the context of the stochastic ensemble Kalman filter and the
Appendix 8.3 as well as Lei & Bickel (2011) for nonlinear extensions.

Even in the mean field limit $J\to\infty$, the ensemble Kalman filter provides
approximations only to approximate transport based filters. As mentioned in
Subsection 1.1, sequential Monte Carlo methods can be designed to be
consistent with the underlying nonlinear filtering problem as defined, for
example, by perfect transport based filters. Foundational analysis of these
particle methods is undertaken in Crisan et al. (1998) and del Moral (2004);
but we reiterate that, in contrast to ensemble Kalman based methods, they do
not scale well to high dimensions. We also point to Del Moral et al. (2006)
for an application of the sequential Monte Carlo method to Bayesian inference
problems; the approach therein is closely related to iterative implementations
of the EnKF that were subsequently developed in the papers Li & Reynolds
(2009), Gu & Oliver (2007), Sakov et al. (2012).

Despite only providing approximations to the exact filtering distribution,
accuracy and stability results for the ensemble Kalman filter, viewed as a
state estimator, have been derived, for example, in González-Tokman & Hunt
(2013), Kelly et al. (2014), Tong et al. (2015), and Tong et al. (2016).
Mechanisms for finite time filter divergence have also been identified
(Gottwald & Majda, 2013; Kelly et al., 2015).

Extending the ensemble Kalman filter to strongly nonlinear and high
dimensional state estimation problems constitute an area of active ongoing
research. The current state of the art has been summarized in van Leeuwen et
al. (2019) in the context of high dimensional geophysical applications.
Extensions of the transport framework (39), which build on approximating the
perfect transport maps $T^{D}$ in (39b) in an asymptotically consistent
manner, include the work by Reich (2013), Cheng & Reich (2015), Spantini et
al. (2019b), and Zech & Marzouk (2022). In an alternative line of research
there have been several proposals to construct hybrid methods, which aim to
adaptively bridge between the ensemble Kalman and particle filters; including
the work by Stordal et al. (2011), Frei & Künsch (2013), and Chustagulprom et
al. (2016).

In this context, the transformation formula (70) proves to be rather useful
since most existing particle based methods can be covered by appropriate
choices of $S_{n}$, where $S_{n}$ is typically the realization of a random
matrix. See Reich & Cotter (2015) for more details. For example, a resampling
step in a sequential Monte Carlo method gives rise to a matrix $S_{n}$ with a
single non-zero entry equal to one in each of its columns. More generally, it
holds that

$\sum_{i\in\mathsf{J}}(S_{n})_{ij}=1.$

The subject of optimal transport is given a comprehensive treatment in Villani
(2008); see also Villani (2021). Computational aspects of the subject are
discussed in Peyré et al. (2019). Ideas from the field are impacting filtering
and Bayesian inference (Corenflos et al., 2021), building on the computational
tractability of entropy-regularized optimal transport (Cuturi, 2013). This
regularization is linked to the Schrödinger bridge problem, and connections
with data assimilation are developed in Reich (2019). See also Subsection 1.1
for discussion of transport-based methodologies within the context of ensemble
Kalman methods. See Example 8.9 and Reich & Cotter (2015) for the connection
between the map (56) and optimal transport. For a derivation of the standard
formulae for mean and covariance of conditioned Gaussians used in the proof of
Lemma 2.5 see Eaton (2007).

##  3 Continuous Time

This section is devoted to deriving, and studying properties of, continuous
time analogues of concepts introduced in the previous section. We start in
Subsection 3.1 defining the set-up. Thereafter the subsections mirror those
from the previous section, describing the relevant continuous time analogues;
in particular we conclude with Subsection 3.7 containing bibliographic notes.

All problems arising in practice are implemented as algorithms in discrete
time, so it is important to establish motivation for the continuous time
formulations. There are two primary reasons for introducing them. The first is
that continuous time limits of the discrete algorithms provide a way to
understand and interpret the behavior of the discrete algorithms; results
about accuracy, stability and uncertainty quantification, which shed light on
the relative merits of different algorithmic approaches, are often cleanest in
the continuous time setting. The second is that many problems arising in
practice involve physical processes which evolve in continuous time; when the
data informing these models is discrete in time, but very high frequency, it
is insightful to consider the idealization of continuous time data. Both of
these motivations underlie the developments in this section.

We introduce a small increment in time, denoted by $\Delta t$. Recalling the
map $\Psi(\cdot)$ defining the systematic component of the state dynamics, we
now define an infinitesimal analogue $f(\cdot)$; we also introduce the
rescaled observation operators $\mathsf{h}(\cdot)$ from the original nonlinear
observation operator $h(\cdot)$; and we introduce state/observational
covariances $(\mathsf{\Gamma},\mathsf{\Sigma})$ by rescaling
$(\Gamma,\Sigma)$:

(72a) $\displaystyle\Psi(v)$ $\displaystyle=v+\Delta tf(v),\quad h(v)=\Delta
t\mathsf{h}(v),$ (72b) $\displaystyle\Sigma$ $\displaystyle=\Delta
t\mathsf{\Sigma},\quad\Gamma=\Delta t\mathsf{\Gamma}.$

In the linear setting $\Psi(\cdot)=M\cdot,h(\cdot)=H\cdot$ we will also
introduce an infinitesimal vector field $f(\cdot)=F\cdot$ for matrix $F$, and
rescaled linear observation operator $\mathsf{H}$ 101010Note the difference
between the discrete time objects $(h(\cdot),H,\Gamma,\Sigma)$ and the related
continuous time objects
$(\mathsf{h}(\cdot),\mathsf{H},\mathsf{\Gamma},\mathsf{\Sigma})$.

$\displaystyle M=\text{Id}+\Delta tF,\quad H=\Delta t\mathsf{H}.$

The observation $\\{y_{n}\\}$ is best thought of, in the scalings we
introduce, as capturing increments of a process $\\{z_{n}\\}$. To capture
this, and extend it to the specific realization of the data appearing in the
algorithms, and the artificial data used in some algorithms, we introduce the
variables $z_{n},{z}^{\dagger}_{n},\widehat{z}_{n}$ by assuming that

(73a) $\displaystyle y_{n+1}$ $\displaystyle:=z_{n+1}-z_{n}=\Delta z_{n+1},$
(73b) $\displaystyle{y}^{\dagger}_{n+1}$
$\displaystyle:={z}^{\dagger}_{n+1}-{z}^{\dagger}_{n}=\Delta{z}^{\dagger}_{n+1},$
(73c) $\displaystyle\widehat{y}_{n+1}$
$\displaystyle:=\widehat{z}_{n+1}-\widehat{z}_{n}=\Delta\widehat{z}_{n+1}.$

Note that $z_{n},{z}^{\dagger}_{n},\widehat{z}_{n}$ have dimension
$d_{z}=d_{y}$. We assume that $z_{0}={z}^{\dagger}_{0}=\widehat{z}_{0}=0$.
Then $z_{n},{z}^{\dagger}_{n},\widehat{z}_{n}$ are uniquely defined from
$y_{n},{y}^{\dagger}_{n},\widehat{y}_{n}$, respectively.

In the following we set $t_{n}=n\Delta t.$ With this definition, and the
scaling above in hand, we may view the state $v_{n}$ and observation $y_{n}$
as relating to approximations of continuous time processes $v(\cdot)$ and
$z(\cdot)$: $v_{n}\approx v(t_{n})$, $z_{n}\approx z(t_{n})$. We also
introduce continuous time process $\widehat{v}(\cdot)$, which we will be used
in the prediction part of algorithms developed in what follows, and
$\bigl{(}{z}^{\dagger}(\cdot),\widehat{z}(\cdot)\bigr{)}$, denoting the
continuous time observed data, which we are conditioning on, and predicted
data, used in the algorithm. We assume that
$z(0)={z}^{\dagger}(0)=\widehat{z}(0)=0.$ Under the rescalings above, and in
the limit $\Delta t\to 0$, the data assimilation problem may be reformulated
in terms of SDEs; and the related mappings on measures, and discrete-time
algorithms that stem from them, may be reformulated in terms of SDEs and
SPDEs; we now go on to identify these continuous time stochastic processes.

###  3.1 Set-Up

We turn our attention to deriving the continuous-time analogue of the
discrete-time set-up (1) for state-observation coevolution: for all
$n\in\mathbb{Z}^{+}$ we have

$\displaystyle v_{n+1}$ $\displaystyle=\Psi(v_{n})+\xi_{n},$ $\displaystyle
y_{n+1}$ $\displaystyle=h(v_{n+1})+\eta_{n+1};$

recall that we assume that $v_{0},\\{\xi_{n}\\}_{n\in\mathbb{Z}^{+}}$ and
$\\{\eta_{n}\\}_{n\in\mathbb{N}}$ are mutually independent Gaussians defined
by

$v_{0}\sim\mathsf{N}(m_{0},C_{0}),\quad\xi_{n}\sim\mathsf{N}(0,\Sigma)\,\,\text{i.i.d.},\quad\eta_{n}\sim\mathsf{N}(0,\Gamma)\,\,\text{i.i.d.}\,\,.$

Applying the rescalings in (72) and the reparametrization of $y_{n+1}$ in
(73), we obtain the system

(74a) $\displaystyle v_{n+1}$ $\displaystyle=v_{n}+\Delta tf(v_{n})+\xi_{n},$
(74b) $\displaystyle z_{n+1}$ $\displaystyle=z_{n}+\Delta
t\mathsf{h}(v_{n+1})+\eta_{n+1},$

for all $n\in\mathbb{Z}^{+}$, where we assume
$v_{0},\\{\xi_{n}\\}_{n\in\mathbb{Z}^{+}}$ and
$\\{\eta_{n}\\}_{n\in\mathbb{N}}$ are mutually independent Gaussians defined
by

$v_{0}\sim\mathsf{N}(m_{0},C_{0}),\quad\xi_{n}\sim\mathsf{N}(0,\Delta
t\mathsf{\Sigma})\,\,\text{i.i.d.},\quad\eta_{n}\sim\mathsf{N}(0,\Delta
t\mathsf{\Gamma})\,\,\text{i.i.d.}\,\,.$

Note that (74) is a variant on the Euler-Maruyama discretization of an SDE.
Indeed, by taking the $\Delta t\rightarrow 0$ limit it is clear that the
natural continuous time analogue of equations (1) is the pair of coupled SDEs

(75a) $\displaystyle dv$ $\displaystyle=f(v)dt+\sqrt{\mathsf{\Sigma}}dW,\quad
v_{0}\sim\mathsf{N}(m_{0},C_{0})$ (75b) $\displaystyle dz$
$\displaystyle=\mathsf{h}(v)dt+\sqrt{\mathsf{\Gamma}}dB,\quad z(0)=0,$

taken to hold for all $t\in\mathbb{R}^{+}.$ The vector fields $f(\cdot)$ and
$\mathsf{h}(\cdot)$ describe the systematic, deterministic components of the
dynamics and observation processes and are assumed known. The systematic
components of the model are subjected to unit Brownian motions $W$ and $B$, in
$\mathbb{R}^{d_{v}}$ and $\mathbb{R}^{d_{z}}$ respectively, and the initial
condition for $v$ is Gaussian. Analogous to the discrete time setting, we
assume that

(76) $C_{0}\succeq 0,\quad\mathsf{\Sigma}\succeq 0,\quad\mathsf{\Gamma}\succ
0,$

where $v_{0},W$ and $B$ are mutually independent random variables. Note that,
for each fixed $t\in\mathbb{R}^{+}$, the state $v(t)\in\mathbb{R}^{d_{v}}$ and
the observations $z(t)\in\mathbb{R}^{d_{z}}.$

Throughout we use $\dagger$ again to denote a specific realization of a
process, as in the discrete time setting. We assume that we have available to