$C_{R}^{-1}=L_{R}^{\top}\mathsf{\Gamma}_{R}^{-1}L_{R}$ and the equilibrium
mean $m_{R}$ solves the equation

$L_{R}^{\top}\mathsf{\Gamma}_{R}^{-1}L_{R}m_{\infty}=L_{R}^{\top}\mathsf{\Gamma}_{R}^{-1}w_{R};$

the discussion following (120) shows that these are the desired expressions
for posterior mean and covariance given by (121) and (122).

######  Remark 4.29.

As for the classical inversion methodology introduced in Subsection 4.4.1 the
convergence to equilibrium is exponential; see Subsection 4.6 for details.
$\blacksquare$

Motivated by this property, we may use Kalman transport algorithms applied to
(173) to solve the Bayesian inverse problem defined by (117) in the general
nonlinear setting. Of course, because Kalman transport is only exact up to
second order moments, this does not produce the exact filtering distribution.
Furthermore, to make resulting algorithms tractable we will predict in (173a)
using the covariance $C_{n}$ of the Kalman transport algorithm, rather than
the covariance under the true filtering distribution. The previous lemma shows
that in the linear case, though, these do coincide. With these considerations
in hand, application of the Kalman transport mean field model (46) to (173)
yields

(174a) $\displaystyle\widehat{u}_{n+1}$ $\displaystyle=u_{n}+\xi_{n},$ (174b)
$\displaystyle\widehat{w}_{n+1}$ $\displaystyle=G_{R}(u_{n+1})+\eta_{n+1},$
(174c) $\displaystyle u_{n+1}$
$\displaystyle=\widehat{u}_{n+1}+\widehat{C}^{uG}_{R,n+1}\bigl{(}\widehat{C}^{GG}_{R,n+1}+(1+\alpha^{-1})\mathsf{\Gamma}_{R}\bigr{)}^{-1}\bigl{(}w_{R}-\widehat{w}_{n+1}\bigr{)}.$

Here $\alpha>0$, $\xi_{n}\sim\mathsf{N}(0,\alpha C_{n}),$
$\eta_{n+1}\sim\mathsf{N}\bigl{(}0,(1+\alpha^{-1})\mathsf{\Gamma}_{R})$ and

$\displaystyle C_{n}$
$\displaystyle=\mathbb{E}\bigl{(}(u_{n}-\mathbb{E}u_{n})\otimes(u_{n}-\mathbb{E}u_{n})\Bigr{)},$
$\displaystyle\widehat{C}^{uG}_{R,n+1}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}\hat{u}_{n+1}-\mathbb{E}\hat{u}_{n+1}\bigr{)}\otimes\bigl{(}G_{R}(\hat{u}_{n+1})-\mathbb{E}G_{R}(\hat{u}_{n+1})\bigr{)}\Bigr{)},$
$\displaystyle\widehat{C}^{GG}_{R,n+1}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}G_{R}(\hat{u}_{n+1})-\mathbb{E}G_{R}(\hat{u}_{n+1})\bigr{)}\otimes\bigl{(}G_{R}(\hat{u}_{n+1})-\mathbb{E}G_{R}(\hat{u}_{n+1})\bigr{)}\Bigr{)}.$

We note that the stochastic dynamics in (174a) can be replaced by the
deterministic update

(175)
$\widehat{u}_{n+1}=\mathbb{E}u_{n}+\sqrt{1+\alpha}\,(u_{n}-\mathbb{E}u_{n})$

leading to an equivalent dynamical behavior, in law, in the case of linear
Gaussian Bayesian inverse problems. The form of (175) is closely related to
multiplicative ensemble inflation as widely used in ensemble Kalman filtering.
Similarly, the stochastic Kalman update step (174c) could be replaced by a
corresponding deterministic one. In this context, it is noteworthy that
$\sigma^{\prime}>0$ in (168) is related to the concept of additive ensemble
inflation as also widely used in ensemble Kalman filtering. See the
Bibliography 4.6 for appropriate references.

####  4.4.4 Steady State Filtering Distribution and Bayesian Inversion:
Continuous Time

We now study continuous time limit of the original state-observation system,
and the second-order Kalman transport algorithm applied to it. We set
$\alpha=\Delta t$ in (173) to obtain the continuous time limit

(176a) $\displaystyle du$ $\displaystyle=\sqrt{C}dW,$ (176b) $\displaystyle
dz$ $\displaystyle=G_{R}(u)dt+\sqrt{\mathsf{\Gamma}_{R}}dB$

where $W,B$ are independent unit Brownian motions of appropriate dimensions
and $C$ is covariance computed under the filtering distribution of $u(t)|Z(t)$
where $Z(t)=\\{z(s)\\}_{0\leq s\leq t}.$ We are interested in the filtering
distribution $u(t)|{Z}^{\dagger}(t)$ in the case where
$\frac{d{z}^{\dagger}}{dt}\equiv w_{R}.$ Now taking the continuous time limit
in equation (174), as before setting $\alpha=\Delta t$, we obtain the
following mean field Kalman transport model for the inverse problem:

(177a) $\displaystyle du$
$\displaystyle=C^{uG}_{R}(\mathsf{\Gamma}_{R})^{-1}(w_{R}\,dt-d\widehat{z})+\sqrt{C}dW,$
(177b) $\displaystyle d\widehat{z}$
$\displaystyle=G_{R}(u)dt+\sqrt{\mathsf{\Gamma}_{R}}dB,$ (177c) $\displaystyle
C$
$\displaystyle=\mathbb{E}\bigl{(}(u-\mathbb{E}u)\otimes(u-\mathbb{E}u)\bigr{)},$
(177d) $\displaystyle C^{uG}_{R}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}u-\mathbb{E}u\bigr{)}\otimes\bigl{(}G_{R}(u)-\mathbb{E}G_{R}(u)\bigr{)}\Bigr{)}.$

Again $W,B$ are independent unit Brownian motions of appropriate dimensions,
and expectation $\mathbb{E}$ is under the law of $u$ itself.

######  Remark 4.30.

We discuss the continuous time mean-field model (177) in the context of Remark
4.19. It again has the character of an approximation of mean-field gradient
descent, using statistical linearization; but now the regularized objective
function $\Phi_{R}$ is used in place of $\Phi.$ Furthermore, in contrast to
the continuous time model (172), it is affine invariant. And finally, in the
linear, Gaussian setting where statistical linearization is exact, we note
that the underlying dynamics differs from (167), but nonetheless has the same
steady state, equal to the posterior distribution. $\blacksquare$

We now study the steady state of (177) in the linear Gaussian setting;
analogous calculations for (167) may be found in the papers cited in
Subsection 4.6.

######  Lemma 4.31.

Assume that $C_{0}\succ 0$ and $\mathsf{\Gamma}\succ 0$. Consider the
filtering distribution $u(t)|Z^{\dagger}(t)$ for the mean-field dynamical
system (176) in the setting where $z^{\dagger}(t)=tw_{R}$ for all $t\geq 0.$
and where $u(0)$ is initialized at a Gaussian.212121 Again, for simplicity of
exposition, we denote this Gaussian as $\mathsf{N}(m_{0},C_{0})$ but emphasize
that here it need not be the prior distribution. Assume that $G(\cdot)=L\cdot$
so that $u(t)|Z^{\dagger}(t)$ is Gaussian $\mathsf{N}(m(t),C(t)).$ Then the
mean-field model (177) has a steady state mean $m_{R}$ and precision matrix
$C^{-1}_{R}$ given by those of the posterior distribution (118) arising from
the Bayesian inverse problem defined by (117).

######  Proof 4.32.

In this linear setting the law of $u$ given by (177) solves the filtering
problem and $u$ itself evolves according to the SDE

$du=CL_{R}^{\top}\Gamma_{R}^{-1}\bigl{(}w_{R}-L_{R}u-\sqrt{\Gamma_{R}}dB\bigr{)}+\sqrt{C}dW,$

where $C$ is the covariance of $u$. Solution $u$ is Gaussian
$\mathsf{N}\bigl{(}m(t),C(t)\bigr{)}$ and use of the Itô formula reveals that

(178a) $\displaystyle\frac{dm}{dt}$
$\displaystyle=CL_{R}^{\top}\mathsf{\Gamma}_{R}^{-1}(w_{R}-L_{R}m)$ (178b)
$\displaystyle\frac{dC}{dt}$
$\displaystyle=-CL_{R}^{\top}\mathsf{\Gamma}_{R}^{-1}L_{R}C+C.$

Setting the right hand sides to zero, and seeking a solution with non-zero
covariance, gives the desired result that $(m_{R},C_{R})$ as defined by (122)
and (121), the posterior mean and covariance in this linear setting.

###  4.5 Ensemble Kalman Methods For Inversion

In this section we make a brief foray into ensemble approximations of the mean
field models introduced in our discussion of inverse problems. We illustrate,
in a succinct fashion, both the transport and iterative approaches of
Subsections 4.2 and 4.4 respectively. By working in specific parametric
regimes we illustrate them both by application of a single algorithm, and we
demonstrate the use of that algorithm in an example.

We note that particle approximations of the mean field models (137) and (170)
yield implementable schemes for inversion. These are described in Algorithms 3
and 4, respectively. Clearly, these two algorithmic frameworks reduce to one
another in the cases $\Delta t=1$ in Algorithm 3 and $\alpha=1$,
$\sigma^{\prime}=0$ and $\gamma^{\prime}=1$ in Algorithm 4, and are then both
derived from ensemble approximations of the mean field model

(179a) $\displaystyle u_{n+1}$
$\displaystyle=u_{n}+C^{uG}_{n}(C^{GG}_{n}+\mathsf{\Gamma})^{-1}\Bigl{(}w+\xi_{n}-G(u_{n})\Bigr{)},$
(179b) $\displaystyle C^{uG}_{n}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}u_{n}-\mathbb{E}u_{n}\bigr{)}\otimes\bigl{(}G(u_{n})-\mathbb{E}G(u_{n})\bigr{)}\Bigr{)},$
(179c) $\displaystyle C^{GG}_{n}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}G(u_{n})-\mathbb{E}G(u_{n})\bigr{)}\otimes\bigl{(}G(u_{n})-\mathbb{E}G(u_{n})\bigr{)}\Bigr{)},$

where $\xi_{n}\sim\text{N}(0,\Gamma)$ and the expectation $\mathbb{E}$ is
computed under the law of $u_{n}$ itself. Algorithm 5 describes a particle
approximation of (179).

In general we recommend using methodology based on (170) with
$\sigma^{\prime}>0$. Furthermore for methodology based on (137) we recommend
using $\Delta t<1.$ However for the specific example we consider below the
algorithm performs well with $\sigma^{\prime}=0$ (for (170)) and with $\Delta
t=1$ (for (137)). Thus by considering ensemble Kalman approximation of (179)
we are able to illustrate both transport and iterative approaches to inversion
within a single numerical experiment.

Regarding the connection to transport, we note that the mean field model (179)
can be derived by setting $\Delta t=1$ in (137), initialized at the prior and
iterated to $N$ such that $N\Delta t=1$ and so $N=1$. Indeed, iterating (179)
initialized at the prior once then yields an approximate sample from the
posterior $\mu^{w}$. The mean field model (179) can also be derived from the
iterative optimization approach (170) by setting $\alpha=1$,
$\sigma^{\prime}=0$ and $\gamma^{\prime}=1$. In this setting (179) can be
initialized anywhere and needs to be iterated to convergence. These
observations explain that implementing a particle approximation of (179),
initialized at the prior, will yield an approximate ensemble of samples from
$\mu^{w}$ after one iteration; and by iterating beyond one step it will
provide approximate solution to classical optimization based framing of the
inverse problem. In the optimization context particle collapse (consensus) may
be taken as a proxy for convergence and used as a stopping criterion, as
demonstrated in the continuous time setting, for linear inverse problems, in
Lemma 4.20.

The mean field model (179) thus constitutes the basis for an implementable
algorithm that can be employed to solve the inverse problem (117) by either
Bayesian or optimzation approaches. Indeed, a particle approximation of (179)
leads to the specific ensemble Kalman inversion method described in Algorithm
5. We have chosen to highlight ensemble Kalman methods based on (179) because
it allows illustration of both the transport and iterative approaches to
solving inverse problems within a single example. Since this coincidence of
the two methods requires the choice $\alpha=1$ in (170) it is natural (see
Remark 4.25) to consider an over-determined inverse problem. The following
illustrative example has precisely this property. In it we study application
of Algorithm 5, and then compare with application of Algorithm 3 with a choice
of $\Delta t<1.$

######  Example 4.33.

We again consider the Lorenz ’96 (singlescale) model for $v\in
C(\mathbb{R}^{+},\mathbb{R}^{L})$ satisfying the equations

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

As before we set $L=9,h_{v}=-0.8$ and $u=10$. We recall that function $m$ is
shown in Figure 2. We have used the notation $u$ instead of $F$, for the
forcing parameter, to align with the notation for the unknown parameter used
through the section concerning inverse problems. Our objective is to recover
parameter $u$ from time-averaged data. We consider an observation
$w\in\mathbb{R}^{2}$ arising from the model

(181) $w=G_{T}(u)+\gamma,$

where $\gamma\sim\text{N}(0,\Gamma)$ represents the noise polluting the data,
$G_{T}:\mathbb{R}\rightarrow R^{2}$ is the forward model that outputs, given
forcing parameter $u$, a two-dimensional vector consisting of the averages
over the $L$ variables, of the mean and variance over time $T$ of the state of
system (180). In particular, considering the state of the system
${v}^{\dagger}$ to evolve according to

(182) ${v}^{\dagger}_{n+1}=\Psi({v}^{\dagger}_{n}),$

where $\Psi$ is the solution operator for (180) over the observation time
interval $\tau$, the action of forward operator $G_{T}$ on $u$ is defined as

$G_{T}(u)=\begin{pmatrix}w_{1}\\\ w_{2}\end{pmatrix},$

with

$\displaystyle w_{1}=\frac{1}{L}\sum_{l=1}^{L}\bar{v}^{\dagger}_{l}\quad
w_{2}=\frac{1}{L\cdot{T/\tau}}\sum_{n=1}^{T/\tau}\sum_{l=1}^{L}(v^{\dagger}_{n;l}-\bar{v}^{\dagger}_{l})^{2},\quad\bar{v}^{\dagger}_{l}=\frac{1}{T/\tau}\sum_{n=1}^{T/\tau}v^{\dagger}_{n;l},$

where we have used $v^{\dagger}_{n;l}$ to denote the $l^{th}$ variable in
vector ${v}^{\dagger}_{n}$, and where the states ${v}^{\dagger}_{n}$ evolve
according to (180) with true forcing parameter $u$, in discrete time form
(182). We note that $G_{T}$ also depends on the initialization of the system
but, by ergodicity, this effect will be small if $T$ is large enough compared
with the Lyapunov time of the system; for this reason we neglect it,
notationally.

To solve the parameter estimation problem consisting of inferring unknown
forcing parameter $u$ from data $w$ given by (181) we use ensemble Kalman
methods as formulated in Algorithms 3 and 5. The data $w$ is generated with
$T=10;$ the algorithms are implemented with $T=20.$ In each instance of the
evaluation of $G_{T}$ the initial condition is chosen at random from a
Gaussian distribution with mean $0$ and variance $40$. We initialize an
ensemble of size $30$ from a prior Gaussian of mean $0$ and variance $10$.

In Figure 9 we show application of Algorithm 5, running for $15$ steps. In
Figure 10 we show application of Algorithm 3 with $\Delta t$ set to $5\cdot
10^{-2}$, running for $40$ steps. Notice that, in the first figure, step $1$
of the algorithm outputs an ensemble Kalman approximation of the posterior and
in the second figure it is at step $20$ that the algorithm does so; the output
of both algorithms at these specific steps corresponds to approximation of the
posterior distribution, given by formula (128) at time $t=1.$ In both cases we
also iterate beyond $t=1$ which takes the algorithms into the realm of
approximate minimization of $\Phi$, as formula (128) shows as $t\to\infty.$
The left panel of both figures shows evolution of the ensemble mean as the
iteration progresses; furthermore the ensemble spread is summarized in the
shading, which denotes a one standard deviation band around the mean. The
right panel shows the distribution of the particles after $1$ and $20$ steps
respectively, and a Gaussian fit with the mean and variance of the particles;
in principle this is an approximation of the posterior distribution. Clearly
both algorithms behave well as an iterative method to solve the inverse
problems from an optimization perspective: the true value of the parameter is
captured well by the ensemble mean; and ensemble collapse, which is a form of
consensus in the interacting particle system, is shown by the shrinking of the
ensemble standard deviation as the iteration progresses (left-hand panels).
The right hand panels illustrate that the putative Bayesian samples from the
posterior should be interpreted with some caution: the distribution is very
different for the two different values of $\Delta t$ that are used and cannot
both be good approximations of the true posterior.

(a)

(b)

Figure 9: The plot in (a) displays the results obtained applying Ensemble
Kalman Inversion (Transport and Iteration to Infinity), i.e. Algorithm 5, to
the forcing parameter estimation problem. The value of the of the true forcing
parameter is displayed in black. The green line represents the parameter
estimates given by the ensemble averages, while the shaded light green region
quantifies one standard deviation around the ensemble mean. The plot in (b)
displays the values of the particles in the ensemble after one iteration step
of Algorithm 5 and a Gaussian fit with the mean and variance of the particles.
This ensemble represents an ensemble Kalman approximate set of samples from
the posterior.

(a)

(b)

Figure 10: The plot in (a) displays the results obtained applying Ensemble
Kalman Inversion (Transport), i.e. Algorithm 3, with $\Delta t=5\cdot 10^{-2}$
to the forcing parameter estimation problem. The value of the of the true
forcing parameter is displayed in black. The green line represents the
parameter estimates given by the ensemble averages, while the shaded light
green region quantifies one standard deviation around the ensemble mean. The
plot in (b) displays the values of the particles in the ensemble after $20$
steps and a Gaussian fit with the mean and variance of the particles. This
ensemble represents an ensemble Kalman approximate set of samples from the
posterior; note that it differs substantially from that obtained with $\Delta
t=1$ in Figure 9.

###  4.6 Bibliographical Notes

Sequential Monte Carlo methods may be used to approximately morph one
probability distribution (source) into another (target), using empirical
approximation and a discrete-time homotopy (Del Moral et al., 2006; Chopin &
Papaspiliopoulos, 2020). In general the methodology does not scale well to
high dimensional problems (Beskos et al., 2014). However some success has been
achieved in this direction (Kantas et al., 2014), and a basic underlying
theory is described in Beskos et al. (2015). Our presentation in this paper is
confined to the setting of ensemble Kalman methods because of their empirical
success and scalability to high dimensions.

The development of ensemble Kalman methods for inverse problems was pioneered
in the study of reservoir simulation, in the context of learning subsurface
properties from localized flow measurements (Chen & Oliver, 2002; Gu & Oliver,
2007; Li & Reynolds, 2009; Emerick & Reynolds, 2013a, b; Evensen, 2018).
Subsequent work has studied parameter estimation in chaotic dynamical systems,
such as those arising in weather forecasting using ensemble methods for joint
state and parameter estimation (Pulido et al., 2018; Bocquet et al., 2020),
and by matching to time-averaged statistics (Schneider et al., 2017; Cleary et
al., 2021; Dunbar et al., 2020), motivated by climate modeling. The chaotic
dynamics that underlie weather and climate models lead to complicated energy
landscapes for minimization and sampling (Lea et al., 2000). The papers Huang
et al. (2022b); Duncan et al. (2021) demonstrate the benefits of using
ensemble methods, rather than computing exact derivatives, for such problems:
the ensemble approach smoothes the energy landscape.

The idea of transporting prior to posterior, as developed in Subsection 4.2,
is widely used in the statistics literature; see Del Moral et al. (2006) for a
unified perspective and for citations to earlier works which characterize the
deformation of one measure to the other either through incrementally building
up the available data, or through a temperature like annealing parameter in
the likelihood. In the context of data assimilation, and ensemble Kalman
methods in particular, these ideas were developed by Li & Reynolds (2009), Gu
& Oliver (2007), and Sakov et al. (2012). Their continuous time limits were
first derived in Bergemann & Reich (2010a) and Bergemann & Reich (2010b) and
further explored in the context of continuous time transport in Reich (2011).
A connection between the non-stochastic Kalman transport equations (142) and
preconditioned gradient descent was first identified in Bergemann & Reich
(2010b) for finite ensemble sizes and in Reich & Cotter (2015) for the mean
field limit. See also Yang et al. (2014) for related formulations based on the
feedback particle filter approach to continuous time filtering.

The papers Schillings & Stuart (2017) and Schillings & Stuart (2018) studied
the use of ensemble Kalman methods for optimization problems, taking a
continuous time limit and making a connection to preconditioned gradient
descent. This led to work on approximate sampling from the preconditioned
Langevin equation in Garbuno-Inigo et al. (2020b), Garbuno-Inigo et al.
(2020a), Nüsken & Reich (2019) and Liu et al. (2022); in particular the papers
Nüsken & Reich (2019); Garbuno-Inigo et al. (2020a) demonstrated an important
finite ensemble size correction to the mean field limit introduced in Garbuno-
Inigo et al. (2020b). The paper Garbuno-Inigo et al. (2020b) used the non-
standard Kalman-Wasserstein metric, first introduced in Reich & Cotter (2015),
to provide a framework to analyze the preconditioned Langevin equation. It
remains open to fully develop the mathematical foundations of gradient flows
using this metric. These papers demonstrate the role that the ensemble plays
in preconditioning the dynamics. This makes a link to the important paper
Goodman & Weare (2010) which introduced the concept of affine invariant
ensemble samplers, an idea developed further in Leimkuhler et al. (2018).

Filtering using dynamical systems with equilibrium filtering distribution
which solves the inverse problem, the approach adopted in Subsection 4.4, is
an idea developed for optimization in Huang et al. (2022b) and for sampling in
Huang et al. (2022a). See those papers for details concerning uniqueness of,
and exponential convergence to, steady state solutions in the linear Gaussian
setting, as discussed in Remarks 4.25, 4.29. See also Pidstrigach & Reich
(2021) for related formulations in the time continuous setting. The stochastic
perturbations utilized in (168a) and (173a) are closely related to additive
and multiplicative, respectively, ensemble inflation methods as widely used in
ensemble Kalman filter implementations (Asch et al., 2016; Evensen et al.,
2022).

The optimization and sampling approaches for inverse problems, in both
discrete and continuous time, can in principle be combined with ideas from
stochastic annealing (Kushner & Yin, 2003), and stochastic gradient descent
(Goodfellow et al., 2016). The papers by Haber et al. (2018); Kovachki &
Stuart (2019); Pidstrigach & Reich (2021) demonstrate the use of ensemble
Kalman methods for inversion, when combined with stochastic gradient descent,
and mini-batching in particular, as well as the application of ensemble Kalman
methods beyond the setting of the $L_{2}$-loss functions $\Phi$ and $\Phi_{R}$
considered here.

However, despite the growing use of ensemble Kalman methods to solve inverse
problems, it is important to appreciate that they all invoke approximations
which amount to matching only first and second order statistics, at some point
in the algorithmic development. For this reason the methods are intuitively
only useful as samplers for problems with posterior distribution close to a
Gaussian. This idea is carefully explained in Ernst et al. (2015). As is the
case for state estimation, analysis is required to justify use of ensemble
Kalman methods beyond the linear and Gaussian regime.

The idea of using ensemble methods for performing the optimization step within
variational data assimilation was introduced in Zupanski (2005). The
connection between iterative applications of the ensemble Kalman filter and
optimization were first investigated in Iglesias et al. (2013) and developed
to include constraints in Albers et al. (2019), Chada et al. (2019) and
Tikhonov regularization Chada et al. (2020b). Recall from Remark 4.22 the
algebraic rates of convergence inherent in mean-field gradient descent. This
undesirable feature of optimization methods based on statistical linearization
of mean-field gradient descent can be ameliorated to some extent by the use of
the use adaptive time-steps, connections to the Levenberg-Marquadt algorithm,
and the use of stopping criteria; see Iglesias (2015), Iglesias (2016), and
Iglesias & Yang (2021). The recent paper Chada et al. (2020a) overviews the
optimization perspective and provides a unifying framework, going beyond
gradient descent-based methods. In particular, framing ensemble Kalman based
optimization methods in terms of statistical linearization, as we do in
Subsection 4.3, originates in that paper.

Other interacting particle system approaches to optimization have been
proposed, including feedback particle (Zhang et al., 2017) and unscented
Kalman approaches (Huang et al., 2022b, a). Furthermore the papers Reich &
Weissmann (2021) and Pavliotis et al. (2021) show how to construct a
derivative free Langevin sampler using localized ensembles. An alternative
interacting particle system approach to solving inverse problems and
optimization tasks is the use of consensus based methods (Tsianos et al.,
2012; Ha et al., 2021; Fornasier et al., 2021; Carrillo et al., 2018, 2021;
Chen et al., 2020; Pinnau et al., 2017; Fornasier et al., 2020).

The papers Ding et al. (2020), Ding & Li (2021a, b) undertake a systematic
analysis of the link between interacting particle systems and mean-field
systems, mostly focused on the solution of inverse problems; however the
methods developed are more widely applicable.

Finally we note that Kalman methods have been related to approximate Bayesian
computation (ABC) methodologies, utilizing a linear regression ansatz (Sisson
et al., 2018; Nott et al., 2012). Such methods are in turn closely related to
Bayes linear and best linear unbiased estimators (BLUE) as discussed in
Goldstein & Rougier (2006), Goldstein & Wooff (2007), Lei & Bickel (2011),
Nott et al. (2012), Snyder (2014), Goldstein (2014), Reich & Cotter (2015),
and Latz (2016). BLUE is discussed in more detail in Appendix 8.3.

##  5 Conclusions and Open Problems

This paper presents a unifying perspective on the derivation, interpretation
and analysis of ensemble Kalman methods through use of the ideas of mean field
models, second order approximate transport and particle approximation. Both
state estimation and parameter estimation (inverse problems) are studied, and
similar ideas may be developed for joint parameter-state estimation problems.
Furthermore, ideas have been presented in discrete time and, through specific
parametric scalings, continuous time limits have been identified. Our unifying
approach constitutes a novel presentation of the subject, and creates a
framework for the mathematical development of the subject area.

Ensemble Kalman methods have been enormously impactful in the geosciences,
where they originated, and are starting to be used in numerous other
application domains. However, if they are to realize their potential for
widespread adoption and application, many research challenges remain. These
challenges are both in mathematical analysis and in the development of
methodology. One of the biggest challenges is the following: some theory, and
abundant numerical evidence, show that ensemble Kalman methods perform well at
state estimation and at parameter estimation; however, there is very little
theory, or empirical evidence, which identifies situations in which the
statistical information in the ensemble constitutes valid approximate Bayesian
inference.

In mathematical analysis a number of substantial challenges are presented by
ensemble Kalman methods, which we list here.

  * •

Determine conditions, related to small noise scenarios, under which the true
state is well-approximated by the mean or sample path of mean field models
based on second order transport; find sharp error estimates.

  * •

More generally, determine conditions under which the filtering distribution is
well-approximated by mean field models based on second order transport; find
sharp error estimates and appropriate metrics for the analysis; furthermore,
identify which problems satisfy these conditions.

  * •

For inverse problems, determine conditions under which the optimizer of a (to-
be-identified) loss function is well-approximated by the mean or sample path
of mean field models based on second order transport; find sharp error
estimates.

  * •

Again, more generally, for inverse problems, determine conditions under which
the Bayesian posterior is well-approximated by mean field models based on
second order transport, in both the transport and iterative approaches to
inversion; find sharp error estimates and appropriate metrics for the
analysis; furthermore, identify which inverse problems satisfy these
conditions.

  * •

For all of the preceding four scenarios derive error bounds for particle
approximations of the mean-field models; when low-rank structure is present in
covariances prove that ensemble Kalman methods can correctly identify it, and
exploit the low-rank structure in the analysis of particle approximations.

  * •

In all of the particle methods arising above, compare the cost/error trade-off
with that arising for other competitive methods, to determine when ensemble
Kalman methods are competitive.

  * •

All of the algorithms in the paper are studied in idealized scenarios, in the
absence of widely employed techniques such as covariance inflation and
localization; developing analyses which account for covariance inflation and
localization will be highly desirable.

On the methodology side there are also a number of significant challenges,
which we also list here.

  * •

Given ability to compute an ensemble of evaluations of the combined state-
observation dynamical system, what is the optimal way to combine them to
either estimate the state given an observation sequence, or the filtering
distribution.

  * •

Given ability to compute an ensemble of evaluations of the forward model, what
is the optimal way to combine them to either estimate the parameter given an
observation, or the posterior distribution, for the corresponding inverse
problem.

  * •

What role might be played by machine learning in addressing the design of
algorithms, and in particular in addressing the preceding questions.

  * •

Develop an overarching interacting particle and mean field framework that
subsumes ensemble Kalman and alternative particle filters as well as
derivative-free and consensus based optimization methods and use the framework
to create new methods.

  * •

Develop principles for the deployment of covariance inflation, and
generalizations of localization, so that the resulting methodology is widely
applicable and does not need application-specific principles to be applied.

  * •

Expand the preceding scenarios beyond the additive error models discussed in
the survey, to include classification and other machine learning tasks.

Acknowledgments EC is grateful to the Kortschak Scholars Program within the
CMS Department at Caltech for financial support. The work of SR was supported
by Deutsche Forschungsgemeinschaft (DFG) - Project-ID 318763901 - SFB1294..
The work of AMS is supported by NSF award AGS1835860, the Office of Naval
Research award N00014-17-1-2079 and by a Department of Defense Vannevar Bush
Faculty Fellowship. The authors are grateful to Dmitry Burov for helpful
advice regarding the numerical experiments; they are also grateful to Ricardo
Baptista and Daniel Sanz-Alonso for helpful discussions which improved the
paper.

## References

  * Abarbanel (2013) Abarbanel, H. (2013)  Predicting the future: Completing models of observed complex systems.  Springer.
  * Abdulle et al. (2012) Abdulle, A., Weinan, E., Engquist, B. & Vanden-Eijnden, E. (2012)  The heterogeneous multiscale method.  Acta Numerica, 21, 1–87.
  * Acevedo et al. (2017) Acevedo, W., de Wiljes, J. & Reich, S. (2017)  Second-order accurate ensemble transform particle filters.  SIAM J. Sci. Comput., 39, A1834–A1850.
  * Agapiou et al. (2017) Agapiou, S., Papaspiliopoulos, O., Sanz-Alonso, D. & Stuart, A. (2017)  Importance sampling: Intrinsic dimension and computational cost.  Statistical Science, 405–431.
  * Albers et al. (2019) Albers, D. J., Blancquart, P.-A., Levine, M. E., Seylabi, E. E. & Stuart, A. (2019)  Ensemble Kalman methods with constraints.  Inverse Problems, 35, 095007.
  * Amezcua et al. (2014) Amezcua, J., Kalnay, E., Ide, K. & Reich, S. (2014)  Ensemble transform Kalman-Bucy filters.  Q.J.R. Meteor. Soc., 140, 995–1004.
  * Anderson & Moore (2012) Anderson, B. D. & Moore, J. B. (2012)  Optimal filtering.  Courier Corporation.
  * Anderson (2001) Anderson, J. L. (2001)  An ensemble adjustment Kalman filter for data assimilation.  Monthly weather review, 129, 2884–2903.
  * Asch et al. (2016) Asch, M., Bocquet, M. & Nodet, M. (2016)  Data assimilation: Methods, algorithms, and applications.  SIAM.
  * Ashwin (2003) Ashwin, P. (2003)  Synchronization from chaos.  Nature, 422, 384–385.
  * Azouani et al. (2014) Azouani, A., Olson, E. & Titi, E. S. (2014)  Continuous data assimilation using general interpolant observables.  Journal of Nonlinear Science, 24, 277–304.
  * Bain & Crisan (2008) Bain, A. & Crisan, D. (2008)  Fundamentals of stochastic filtering, vol. 60.  Springer Science & Business Media.
  * Bergemann & Reich (2010a) Bergemann, K. & Reich, S. (2010a)  A localization technique for ensemble Kalman filters.  Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography, 136, 701–707.
  * Bergemann & Reich (2010b) Bergemann, K. & Reich, S. (2010b)  A mollified ensemble Kalman filter.  Quarterly Journal of the Royal Meteorological Society, 136, 1636–1643.
  * Bergemann & Reich (2012) Bergemann, K. & Reich, S. (2012)  An ensemble Kalman-Bucy filter for continuous data assimilation.  Meteorologische Zeitschrift, 21, 213.
  * Beskos et al. (2014) Beskos, A., Crisan, D. & Jasra, A. (2014)  On the stability of sequential monte carlo methods in high dimensions.  The Annals of Applied Probability, 24, 1396–1445.
  * Beskos et al. (2015) Beskos, A., Jasra, A., Muzaffer, E. A. & Stuart, A. M. (2015)  Sequential monte carlo methods for bayesian elliptic inverse problems.  Statistics and Computing, 25, 727–737.
  * Bickel et al. (2008) Bickel, P., Li, B. & Bengtsson, T. (2008)  Sharp failure rates for the bootstrap particle filter in high dimensions.  Pushing the limits of contemporary statistics: Contributions in honor of Jayanta K. Ghosh.  Institute of Mathematical Statistics, pp. 318–329.
  * Bishop et al. (2001) Bishop, C., Etherton, B. & Majumdar, S. (2001)  Adaptive sampling withe ensemble transform Kalman filter. Part I: Theoretical aspects.  Monthly Weather Review, 129, 420–436.
  * Bishop (2011) Bishop, C. (2011)  Pattern Recognition and Machine Learning, 2nd edn.  New York: Springer-Verlag.
  * Blömker et al. (2013) Blömker, D., Law, K., Stuart, A. M. & Zygalakis, K. C. (2013)  Accuracy and stability of the continuous-time 3DVAR filter for the Navier–Stokes equation.  Nonlinearity, 26, 2193.
  * Blömker et al. (2018) Blömker, D., Schillings, C. & Wacker, P. (2018)  A strongly convergent numerical scheme from ensemble Kalman inversion.  SIAM Journal on Numerical Analysis, 56, 2537–2562.
  * Blömker et al. (2021) Blömker, D., Schillings, C., Wacker, P. & Weissmann, S. (2021)  Continuous time limit of the stochastic ensemble Kalman inversion: Strong convergence analysis.  Technical Report.  arXiv:2107.14508.
  * Bocquet et al. (2017) Bocquet, M., Gurumoorthy, K. S., Apte, A., Carrassi, A., Grudzien, C. & Jones, C. K. (2017)  Degenerate Kalman filter error covariances and their convergence onto the unstable subspace.  SIAM/ASA Journal on Uncertainty Quantification, 5, 304–333.
  * Bocquet et al. (2020) Bocquet, M., Brajard, J., Carrassi, A. & Bertino, L. (2020)  Bayesian inference of chaotic dynamics by merging data assimilation, machine learning and expectation-maximization.  Foundations of Data Science, 2, 55–80.
  * Bocquet & Sakov (2012) Bocquet, M. & Sakov, P. (2012)  Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems.  Nonlinear Processes in Geophysics, 19, 383–399.
  * Bocquet & Sakov (2014) Bocquet, M. & Sakov, P. (2014)  An iterative ensemble Kalman smoother.  Q. J. R. Meteorol. Soc., 140, 1521–1535.
  * Burgers et al. (1998) Burgers, G., van Leeuwen, P. & Evensen, G. (1998)  Analysis scheme in the ensemble Kalman filter.  Monthly Weather Review, 126, 1719–1724.
  * Burov et al. (2021) Burov, D., Giannakis, D., Manohar, K. & Stuart, A. (2021)  Kernel analog forecasting: Multiscale test problems.  Multiscale Modeling & Simulation, 19, 1011–1040.
  * Carrillo et al. (2021) Carrillo, J., Hoffmann, F., Stuart, A. & Vaes, U. (2021)  Consensus based sampling.  arXiv preprint arXiv:2106.02519.
  * Carrillo et al. (2018) Carrillo, J. A., Choi, Y.-P., Totzeck, C. & Tse, O. (2018)  An analytical framework for consensus-based global optimization method.  Mathematical Models and Methods in Applied Sciences, 28, 1037–1066.
  * Chada et al. (2019) Chada, N. K., Schillings, C. & Weissmann, S. (2019)  On the incorporation of box-constraints for ensemble Kalman inversion.  arXiv preprint arXiv:1908.00696.
  * Chada et al. (2020a) Chada, N. K., Chen, Y. & Sanz-Alonso, D. (2020a)  Iterative ensemble Kalman methods: a unified perspective with some new variants.  arXiv preprint arXiv:2010.13299.
  * Chada et al. (2020b) Chada, N. K., Stuart, A. M. & Tong, X. T. (2020b)  Tikhonov regularization within ensemble Kalman inversion.  SIAM Journal on Numerical Analysis, 58, 1263–1294.
  * Chen et al. (2020) Chen, J., Jin, S. & Lyu, L. (2020)  A consensus-based global optimization method with adaptive momentum estimation.  arXiv preprint arXiv:2012.04827.
  * Chen et al. (2022) Chen, Y., Sanz-Alonso, D. & Willett, R. (2022)  Autodifferentiable Ensemble Kalman filters.  SIAM Journal on Mathematics of Data Science, 4, 801–833.
  * Chen & Oliver (2002) Chen, Y. & Oliver, D. (2002)  Ensemble randomized maximum likelihood method as an iterative ensemble smoother.  Mathematical Geosciences, 44, 1–26.
  * Cheng & Reich (2015) Cheng, Y. & Reich, S. (2015)  Data assimilation: A coupling of measure perspective.  Frontiers in Applied Dynamical Systems, vol. 2.  New York: Springer-Verlag.
  * Chernov et al. (2021) Chernov, A., Hoel, H., Law, K. J., Nobile, F. & Tempone, R. (2021)  Multilevel ensemble kalman filtering for spatio-temporal processes.  Numerische Mathematik, 147, 71–125.
  * Chopin & Papaspiliopoulos (2020) Chopin, N. & Papaspiliopoulos, O. (2020)  An Introduction to Sequential Monte Carlo.  Cham, Switzerland: Springer Nature Switzerland AG.
  * Chustagulprom et al. (2016) Chustagulprom, N., Reich, S. & Reinhardt, M. (2016)  A hybrid ensemble transform filter for nonlinear and spatially extended dynamical systems.  SIAM/ASA J. Uncertainty Quantification, 4, 592–608.
  * Clark & Crisan (2005) Clark, J. & Crisan, D. (2005)  On a robust version of the integral representation formula of nonlinear filtering.  Probab. Theory Relat. Fields, 133, 43–56.
  * Cleary et al. (2021) Cleary, E., Garbuno-Inigo, A., Lan, S., Schneider, T. & Stuart, A. M. (2021)  Calibrate, emulate, sample.  Journal of Computational Physics, 424, 109716.
  * Coghi et al. (2021) Coghi, M., Nilssen, T., Nüsken, N. & Reich, S. (2021)  Rough McKean–Vlasov dynamics for robust ensemble Kalman filtering.  Technical Report.  arXiv:2107.06621.
  * Corenflos et al. (2021) Corenflos, A., Thornton, J., Deligiannidis, G. & Doucet, A. (2021)  Differentiable particle filtering via entropy-regularized optimal transport.  International Conference on Machine Learning.  PMLR, pp. 2100–2111.
  * Cotter & Reich (2013) Cotter, C. & Reich, S. (2013)  Ensemble filter techniques for intermittent data assimilation.  Radon Ser. Comput. Appl. Math., 13, 91–134.
  * Crisan et al. (1998) Crisan, D., Del Moral, P. & Lyons, T. (1998)  Discrete filtering using branching and interacting particle systems.  Citeseer.
  * Crisan & Lyons (1999) Crisan, D. & Lyons, T. (1999)  A particle approximation of the solution of the kushner–stratonovitch equation.  Probability Theory and Related Fields, 115, 549–578.
  * Crisan & Xiong (2010) Crisan, D. & Xiong, J. (2010)  Approximate McKean-Vlasov representations for a class of SPDEs.  Stochastics, 82, 53–68.
  * Cuturi (2013) Cuturi, M. (2013)  Sinkhorn distances: Lightspeed computation of optimal transport.  Advances in neural information processing systems, 26.
  * Daum et al. (2010) Daum, F., Huang, J. & Noushin, A. (2010)  Exact particle flow for nonlinear filters.  Signal Processing, Sensor Fusion, and Target Recognition XIX (I. Kadar ed.), vol. 7697.  International Society for Optics and Photonics, SPIE, pp. 92 – 110.
  * de Wiljes et al. (2018) de Wiljes, J., Reich, S. & Stannat, W. (2018)  Long-time stability and accuracy of the ensemble Kalman–Bucy filter for fully observed processes and small measurement noise.  SIAM Journal on Applied Dynamical Systems, 17, 1152–1181.
  * de Wiljes & Tong (2020) de Wiljes, J. & Tong, X. (2020)  Analysis of a localised nonlinear ensemble Kalman Bucy filter with complete and accurate observations.  Nonlinearity, 33, 4752–4782.
  * del Moral (2004) del Moral, P. (2004)  Feynman–Kac Formulae: Genealogical and Interacting Particle Systems with Applications.  New York: Springer-Verlag.
  * Del Moral et al. (2006) Del Moral, P., Doucet, A. & Jasra, A. (2006)  Sequential Monte Carlo samplers.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68, 411–436.
  * Del Moral & Tugaut (2018) Del Moral, P. & Tugaut, J. (2018)  On the stability and the uniform propagation of chaos properties of ensemble Kalman–Bucy filters.  The Annals of Applied Probability, 28, 790–850.
  * Ding et al. (2020) Ding, Z., Li, Q. & Lu, J. (2020)  Ensemble Kalman inversion for nonlinear problems: weights, consistency, and variance bounds.  arXiv preprint arXiv:2003.02316.
  * Ding & Li (2021a) Ding, Z. & Li, Q. (2021a)  Ensemble Kalman inversion: mean-field limit and convergence analysis.  Statistics and Computing, 31, 1–21.
  * Ding & Li (2021b) Ding, Z. & Li, Q. (2021b)  Ensemble Kalman sampler: Mean-field limit and convergence analysis.  SIAM Journal on Mathematical Analysis, 53, 1546–1578.
  * Doucet et al. (2001) Doucet, A., De Freitas, N. & Gordon, N. (2001)  An introduction to sequential Monte Carlo methods.  Sequential Monte Carlo Methods in Practice.  Springer, pp. 3–14.
  * Dunbar et al. (2020) Dunbar, O. R., Garbuno-Inigo, A., Schneider, T. & Stuart, A. M. (2020)  Calibration and uncertainty quantification of convective parameters in an idealized GCM.  arXiv preprint arXiv:2012.13262.
  * Duncan et al. (2021) Duncan, A., Stuart, A. & Wolfram, M.-T. (2021)  Ensemble inference methods for models with noisy and expensive likelihoods.  arXiv preprint arXiv:2104.03384.
  * Eaton (2007) Eaton, M. L. (2007)  Multivariate statistics: A vector space approach.  Lecture Notes-Monograph Series, 53, i–512.
  * Eberle (2013) Eberle, A. (2013)  Stochastic analysis.  Technical Report.  Lecture Notes, University of Bonn.
  * El Moselhy & Marzouk (2012) El Moselhy, T. A. & Marzouk, Y. M. (2012)  Bayesian inference with optimal maps.  Journal of Computational Physics, 231, 7815–7850.
  * Emerick & Reynolds (2013a) Emerick, A. A. & Reynolds, A. C. (2013a)  Ensemble smoother with multiple data assimilation.  Computers & Geosciences, 55, 3–15.
  * Emerick & Reynolds (2013b) Emerick, A. & Reynolds, A. (2013b)  Investigation of the sampling performance of ensemble-based methods with a simple reservoir model.  Computational Geosciences, 17, 325–350.
  * Ernst et al. (2015) Ernst, O. G., Sprungk, B. & Starkloff, H.-J. (2015)  Analysis of the ensemble and polynomial chaos Kalman filters in bayesian inverse problems.  SIAM/ASA Journal on Uncertainty Quantification, 3, 823–851.
  * Evans (2012) Evans, L. C. (2012)  An introduction to stochastic differential equations, vol. 82.  American Mathematical Soc.
  * Evensen (1994) Evensen, G. (1994)  Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics.  Journal of Geophysical Research: Oceans, 99, 10143–10162.
  * Evensen (2018) Evensen, G. (2018)  Analysis of iterative ensemble smoothers for solving inverse problems.  Computational Geosciences, 22, 885–908.
  * Evensen (2019) Evensen, G. (2019)  Accounting for model errors in iterative ensemble smoothers.  Computational Geosciences, 23, 761–775.
  * Evensen et al. (2022) Evensen, G., Vossepoel, F. & van Leeuwen, P. (2022)  Data Assimilation Fundamentals: A unified Formulation of the State and Parameter Estimation Problem.  Cham, Switzerland: Springer Nature Switzerland AG.
  * Fatkullin & Vanden-Eijnden (2004) Fatkullin, I. & Vanden-Eijnden, E. (2004)  A computational strategy for multiscale systems with applications to Lorenz 96 model.  Journal of Computational Physics, 200, 605–638.
  * Fertig et al. (2007) Fertig, E. J., Harlim, J. & Hunt, B. R. (2007)  A comparative study of 4D-VAR and a 4D ensemble Kalman filter: Perfect model simulations with Lorenz-96.  Tellus A: Dynamic Meteorology and Oceanography, 59, 96–100.
  * Foias et al. (2016) Foias, C., Mondaini, C. F. & Titi, E. S. (2016)  A discrete data assimilation scheme for the solutions of the two-dimensional Navier–Stokes equations and their statistics.  SIAM Journal on Applied Dynamical Systems, 15, 2109–2142.
  * Fornasier et al. (2020) Fornasier, M., Huang, H., Pareschi, L. & Sünnen, P. (2020)  Consensus-based optimization on hypersurfaces: Well-posedness and mean-field limit.  Mathematical Models and Methods in Applied Sciences, 30, 2725–2751.
  * Fornasier et al. (2021) Fornasier, M., Klock, T. & Riedl, K. (2021)  Consensus-based optimization methods converge globally in mean-field law.  arXiv preprint arXiv:2103.15130.
  * Frei & Künsch (2013) Frei, M. & Künsch, H. R. (2013)  Bridging the ensemble Kalman and particle filters.  Biometrika, 100, 781–800.
  * Garbuno-Inigo et al. (2020a) Garbuno-Inigo, A., Nüsken, N. & Reich, S. (2020a)  Affine invariant interacting Langevin dynamics for bayesian inference.  SIAM Journal on Applied Dynamical Systems, 19, 1633–1658.
  * Garbuno-Inigo et al. (2020b) Garbuno-Inigo, A., Hoffmann, F., Li, W. & Stuart, A. M. (2020b)  Interacting Langevin diffusions: Gradient structure and ensemble Kalman sampler.  SIAM Journal on Applied Dynamical Systems, 19, 412–441.
  * Gesho et al. (2016) Gesho, M., Olson, E. & Titi, E. S. (2016)  A computational study of a data assimilation algorithm for the two-dimensional Navier–Stokes equations.  Communications in Computational Physics, 19, 1094–1110.
  * Ghil et al. (1981) Ghil, M., Cohn, S., Tavantzis, J., Bube, K. & Isaacson, E. (1981)  Applications of estimation theory to numerical weather prediction.  Dynamic meteorology: Data assimilation methods.  Springer, pp. 139–224.
  * Goldstein (2014) Goldstein, M. (2014)  Bayes linear analysis.  Wiley StatsRef: Statistics Reference Online, 1–7.
  * Goldstein & Rougier (2006) Goldstein, M. & Rougier, J. (2006)  Bayes linear calibrated prediction for complex systems.  Journal of the American Statistical Association, 101, 1132–1143.
  * Goldstein & Wooff (2007) Goldstein, M. & Wooff, D. (2007)  Bayes linear statistics: Theory and methods, vol. 716.  John Wiley & Sons.
  * Gonzalez & Stuart (2008) Gonzalez, O. & Stuart, A. M. (2008)  A first course in continuum mechanics.  Cambridge University Press.
  * González-Tokman & Hunt (2013) González-Tokman, C. & Hunt, B. R. (2013)  Ensemble data assimilation for hyperbolic systems.  Physica D: Nonlinear Phenomena, 243, 128–142.
  * Goodfellow et al. (2016) Goodfellow, I., Bengio, Y., Courville, A. & Bengio, Y. (2016)  Deep learning, vol. 1.  MIT press Cambridge.
  * Goodman & Weare (2010) Goodman, J. & Weare, J. (2010)  Ensemble samplers with affine invariance.  Communications in applied mathematics and computational science, 5, 65–80.
  * Gottwald & Majda (2013) Gottwald, G. A. & Majda, A. J. (2013)  A mechanism for catastrophic filter divergence in data assimilation for sparse observation networks.  Nonlinear Processes in Geophysics, 20, 705–712.
  * Gottwald & Reich (2021) Gottwald, G. A. & Reich, S. (2021)  Supervised learning from noisy observations: Combining machine-learning techniques with data assimilation.  Physica D: Nonlinear Phenomena, 423, 132911.
  * Grooms et al. (2014) Grooms, I., Lee, Y. & Majda, A. J. (2014)  Ensemble Kalman filters for dynamical systems with unresolved turbulence.  Journal of Computational Physics, 273, 435–452.
  * Grooms et al. (2015) Grooms, I., Lee, Y. & Majda, A. J. (2015)  Ensemble filtering and low-resolution model error: Covariance inflation, stochastic parameterization, and model numerics.  Monthly Weather Review, 143, 3912–3924.
  * Grooms (2021) Grooms, I. (2021)  Analog ensemble data assimilation and a method for constructing analogs with variational autoencoders.  Quarterly Journal of the Royal Meteorological Society, 147, 139–149.
  * Gu & Oliver (2007) Gu, Y. & Oliver, D. S. (2007)  An iterative ensemble Kalman filter for multiphase fluid flow data assimilation.  SPE Journal, 12, 438–446.
  * Gurtin (1982) Gurtin, M. E. (1982)  An introduction to continuum mechanics.  Academic press.
  * Gurumoorthy et al. (2017) Gurumoorthy, K. S., Grudzien, C., Apte, A., Carrassi, A. & Jones, C. K. (2017)  Rank deficiency of Kalman error covariance matrices in linear time-varying system with deterministic evolution.  SIAM Journal on Control and Optimization, 55, 741–759.
  * Guth et al. (2020) Guth, P. A., Schillings, C. & Weissmann, S. (2020)  Ensemble Kalman filter for neural network based one-shot inversion.  arXiv preprint arXiv:2005.02039.
  * Ha et al. (2021) Ha, S.-Y., Jin, S. & Kim, D. (2021)  Convergence and error estimates for time-discrete consensus-based optimization algorithms.  Numerische Mathematik, 147, 255–282.
  * Haber et al. (2018) Haber, E., Lucka, F. & Ruthotto, L. (2018)  Never look back - The EnKF method and its application to the training of neural networks without back propagation.  arXiv preprint arXiv:1805.08034.
  * Harlim et al. (2014) Harlim, J., Mahdi, A. & Majda, A. J. (2014)  An ensemble Kalman filter for statistical estimation of physics constrained nonlinear regression models.  Journal of Computational Physics, 257, 782–812.
  * Harlim & Hunt (2007a) Harlim, J. & Hunt, B. R. (2007a)  Four-dimensional local ensemble transform Kalman filter: Numerical experiments with a global circulation model.  Tellus A: Dynamic Meteorology and Oceanography, 59, 731–748.
  * Harlim & Hunt (2007b) Harlim, J. & Hunt, B. R. (2007b)  A non-Gaussian ensemble filter for assimilating infrequent noisy observations.  Tellus A: Dynamic Meteorology and Oceanography, 59, 225–237.
  * Harlim & Majda (2010) Harlim, J. & Majda, A. J. (2010)  Filtering turbulent sparsely observed geophysical flows.  Monthly Weather Review, 138, 1050–1083.
  * Hayden et al. (2011) Hayden, K., Olson, E. & Titi, E. S. (2011)  Discrete data assimilation in the Lorenz and 2D Navier–Stokes equations.  Physica D: Nonlinear Phenomena, 240, 1416–1425.
  * Higham (2001) Higham, D. J. (2001)  An algorithmic introduction to numerical simulation of stochastic differential equations.  SIAM review, 43, 525–546.
  * Hoel et al. (2016) Hoel, H., Law, K. J. & Tempone, R. (2016)  Multilevel ensemble kalman filtering.  SIAM Journal on Numerical Analysis, 54, 1813–1839.
  * Houtekamer & Mitchell (2005) Houtekamer, P. & Mitchell, H. (2005)  Ensemble Kalman filtering.  Q. J. Royal Meteorological Soc., 131, 3269–3289.
  * Houtekamer & Mitchell (1998) Houtekamer, P. & Mitchell, M. (1998)  Data assimilation using an ensemble Kalman filter techniques.  Monthly Weather Review, 126, 796–811.
  * Hu et al. (2002) Hu, Y., Kallianpur, G. & Xiong, J. (2002)  An approximation for the Zakai equation.  Applied Mathematics and Optimization, 45, 23–44.
  * Huang et al. (2022a) Huang, D. Z., Huang, J., Reich, S. & Stuart, A. M. (2022a)  Efficient derivative-free bayesian inference for large-scale inverse problems.  arXiv preprint arXiv:2204.04386.
  * Huang et al. (2022b) Huang, D. Z., Schneider, T. & Stuart, A. M. (2022b)  Iterated Kalman methodology for inverse problems.  Journal of Computational Physics, 463, 111262.
  * Hunt et al. (2007) Hunt, B., Kostelich, E. & Szunyogh, I. (2007)  Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter.  Physica D, 230, 112–126.
  * Iglesias et al. (2013) Iglesias, M. A., Law, K. J. & Stuart, A. M. (2013)  Ensemble Kalman methods for inverse problems.  Inverse Problems, 29, 045001.
  * Iglesias (2015) Iglesias, M. A. (2015)  Iterative regularization for ensemble data assimilation in reservoir models.  Computational Geosciences, 19, 177–212.
  * Iglesias (2016) Iglesias, M. A. (2016)  A regularizing iterative ensemble Kalman method for pde-constrained inverse problems.  Inverse Problems, 32, 025002.
  * Iglesias & Yang (2021) Iglesias, M. & Yang, Y. (2021)  Adaptive regularisation for ensemble Kalman inversion.  Inverse Problems, 37, 025008.
  * Jazwinski (2007) Jazwinski, A. H. (2007)  Stochastic processes and filtering theory.  Courier Corporation.
  * Julier et al. (2000) Julier, S., Uhlmann, J. & Durrant-Whyte, H. F. (2000)  A new method for the nonlinear transformation of means and covariances in filters and estimators.  IEEE Transactions on Automatic Control, 45, 477–482.
  * Kaipio & Somersalo (2006) Kaipio, J. & Somersalo, E. (2006)  Statistical and computational inverse problems, vol. 160.  Springer Science & Business Media.
  * Kalman (1960) Kalman, R. E. (1960)  A New Approach to Linear Filtering and Prediction Problems.  Journal of Basic Engineering, 82, 35–45.
  * Kalman & Bucy (1961) Kalman, R. E. & Bucy, R. S. (1961)  New Results in Linear Filtering and Prediction Theory.  Journal of Basic Engineering, 83, 95–108.
  * Kantas et al. (2014) Kantas, N., Beskos, A. & Jasra, A. (2014)  Sequential monte carlo methods for high-dimensional inverse problems: A case study for the navier–stokes equations.  SIAM/ASA Journal on Uncertainty Quantification, 2, 464–489.
  * Kelly et al. (2015) Kelly, D., Majda, A. J. & Tong, X. T. (2015)  Concrete ensemble Kalman filters with rigorous catastrophic filter divergence.  Proceedings of the National Academy of Sciences, 112, 10589–10594.
  * Kelly et al. (2014) Kelly, D. T., Law, K. J. & Stuart, A. M. (2014)  Well-posedness and accuracy of the ensemble Kalman filter in discrete and continuous time.  Nonlinearity, 27, 2579.
  * Kitanidis (1995) Kitanidis, P. K. (1995)  Quasi-linear geostatistical theory for inversion.  Water Resour. Res., 31, 2411–2419.
  * Kloeden & Platen (1991) Kloeden, P. & Platen, E. (1991)  Numerical methods for stochastic differential equations.  New York: Springer.
  * Kovachki & Stuart (2019) Kovachki, N. B. & Stuart, A. M. (2019)  Ensemble Kalman inversion: A derivative-free technique for machine learning tasks.  Inverse Problems, 35, 095005.
  * Kushner & Yin (2003) Kushner, H. & Yin, G. G. (2003)  Stochastic approximation and recursive algorithms and applications, vol. 35.  Springer Science & Business Media.
  * Kuznetsov et al. (2003) Kuznetsov, L., Ide, K. & Jones, C. K. (2003)  A method for assimilation of Lagrangian data.  Monthly Weather Review, 131, 2247–2260.
  * Kwiatkowski & Mandel (2015) Kwiatkowski, E. & Mandel, J. (2015)  Convergence of the square root ensemble Kalman filter in the large ensemble limit.  SIAM/ASA Journal on Uncertainty Quantification, 3, 1–17.
  * Lange & Stannat (2019) Lange, T. & Stannat, W. (2019)  On the continuous time limit of Ensemble Square Root Filters.  arXiv preprint arXiv:1910.12493.
  * Lange & Stannat (2021) Lange, T. & Stannat, W. (2021)  On the continuous time limit of the ensemble Kalman filter.  Mathematics of Computation, 90, 233–265.
  * Larios & Pei (2017) Larios, A. & Pei, Y. (2017)  Nonlinear continuous data assimilation.  arXiv preprint arXiv:1703.03546.
  * Latz (2016) Latz, J. (2016)  Bayes linear methods for inverse problems.  Ph.D. thesis, Master’s thesis, University of Warwick.
  * Law et al. (2015) Law, K., Stuart, A. & Zygalakis, K. (2015)  Data assimilation.  Cham, Switzerland: Springer, 214.
  * Law et al. (2016a) Law, K., Sanz-Alonso, D., Shukla, A. & Stuart, A. (2016a)  Filter accuracy for the Lorenz 96 model: Fixed versus adaptive observation operators.  Physica D: Nonlinear Phenomena, 325, 1–13.
  * Law et al. (2012) Law, K. J., Shukla, A. & Stuart, A. M. (2012)  Analysis of the 3DVAR filter for the partially observed Lorenz’63 model.  arXiv preprint arXiv:1212.4923.
  * Law et al. (2016b) Law, K. J., Tembine, H. & Tempone, R. (2016b)  Deterministic mean-field ensemble kalman filtering.  SIAM Journal on Scientific Computing, 38, A1251–A1279.
  * Law & Stuart (2012) Law, K. J. & Stuart, A. M. (2012)  Evaluating data assimilation algorithms.  Monthly Weather Review, 140, 3757–3782.
  * Le Gland et al. (2011) Le Gland, F., Monbet, V. & Tran, V.-D. (2011)  Large sample asymptotics for the ensemble Kalman filter.  Oxford: Oxford Univ. Press, pp. 598–631.
  * Lea et al. (2000) Lea, D. J., Allen, M. R. & Haine, T. W. (2000)  Sensitivity analysis of the climate of a chaotic system.  Tellus A: Dynamic Meteorology and Oceanography, 52, 523–532.
  * Lee et al. (2017) Lee, Y., Majda, A. J. & Qi, D. (2017)  Preventing catastrophic filter divergence using adaptive additive inflation for baroclinic turbulence.  Monthly Weather Review, 145, 669–682.
  * Lei & Bickel (2011) Lei, J. & Bickel, P. (2011)  A moment matching ensemble filter for nonlinear non-Gaussian data assimilation.  Monthly Weather Review, 139, 3964–3973.
  * Leimkuhler et al. (2018) Leimkuhler, B., Matthews, C. & Weare, J. (2018)  Ensemble preconditioning for markov chain monte carlo simulation.  Statistics and Computing, 28, 277–290.
  * Li & Reynolds (2009) Li, G. & Reynolds, A. C. (2009)  An iterative ensemble Kalman filter for data assimilation.  Society of Petroleum Engineers, pp. 496–505.
  * Liu et al. (2022) Liu, Z., Stuart, A. & Wang, Y. (2022)  Second order ensemble Langevin method for sampling and inverse problems.  arXiv:2208.04506.
  * Livings et al. (2008) Livings, D., Dance, S. & Nichols, N. (2008)  Unbiased ensemble square root filters.  Physica D, 237, 1021–1028.
  * Lorenz (1996) Lorenz, E. N. (1996)  Predictability: A problem partly solved.  Proc. Seminar on Predictability, vol. 1.  ECMWF.
  * Luenberger (1971) Luenberger, D. (1971)  An introduction to observers.  IEEE Transactions on automatic control, 16, 596–602.
  * Luenberger (1964) Luenberger, D. G. (1964)  Observing the state of a linear system.  IEEE transactions on military electronics, 8, 74–80.
  * Majda & Tong (2018) Majda, A. J. & Tong, X. T. (2018)  Performance of ensemble Kalman filters in large dimensions.  Communications on Pure and Applied Mathematics, 71, 892–937.
  * Mandel et al. (2011) Mandel, J., Cobb, L. & Beezley, J. (2011)  On the convergence of the ensemble Kalman filter.  Appl Math, 56, 533–541.
  * McLachlan & Quispel (2002) McLachlan, R. I. & Quispel, G. R. W. (2002)  Splitting methods.  Acta Numerica, 11, 341–434.
  * Mondaini & Titi (2018) Mondaini, C. F. & Titi, E. S. (2018)  Uniform-in-time error estimates for the postprocessing Galerkin method applied to a data assimilation algorithm.  SIAM Journal on Numerical Analysis, 56, 78–110.
  * Nott et al. (2012) Nott, D. J., Marshall, L. & Ngoc, T. M. (2012)  The ensemble Kalman filter is an ABC algorithm.  Statistics and Computing, 22, 1273–1276.
  * Nüsken & Reich (2019) Nüsken, N. & Reich, S. (2019)  Note on interacting Langevin diffusions: Gradient structure and ensemble Kalman sampler by Garbuno-Inigo, Hoffmann, Li and Stuart.  arXiv preprint arXiv:1908.10890.
  * Oksendal (2013) Oksendal, B. (2013)  Stochastic differential equations: an introduction with applications.  Springer Science & Business Media.
  * Oliver et al. (2008) Oliver, D., Reynolds, A. & Liu, N. (2008)  Inverse Theory for Petroleum Reservoir Characterization and History Matching.  Cambridge: Cambridge University Press.
  * Oliver et al. (1997) Oliver, D. S., Cunha, L. B. & Reynolds, A. C. (1997)  Markov chain Monte Carlo methods for conditioning a permeability field to pressure data.  Mathematical Geology, 29, 61–91.
  * Olson & Titi (2003) Olson, E. & Titi, E. S. (2003)  Determining modes for continuous data assimilation in 2D turbulence.  Journal of statistical physics, 113, 799–840.
  * Pathiraja et al. (2021) Pathiraja, S., Reich, S. & Stannat, W. (2021)  McKean-Vlasov SDEs in nonlinear filtering.  SIAM J Control and Optimization, 59, 4188–4212.
  * Pavliotis et al. (2021) Pavliotis, G., Stuart, A. & Vaes, U. (2021)  Derivative-free Bayesian inversion using multiscale dynamics.  arXiv preprint arXiv:2102.00540.
  * Pavliotis & Stuart (2008) Pavliotis, G. & Stuart, A. (2008)  Multiscale methods: Averaging and homogenization.  Springer Science & Business Media.
  * Peyré et al. (2019) Peyré, G., Cuturi, M. et al. (2019)  Computational optimal transport: With applications to data science.  Foundations and Trends® in Machine Learning, 11, 355–607.
  * Pidstrigach & Reich (2021) Pidstrigach, J. & Reich, S. (2021)  Affine-invariant ensemble transform methods for logistic regression.  arXiv preprint arXiv:2104.08061.
  * Pinnau et al. (2017) Pinnau, R., Totzeck, C., Tse, O. & Martin, S. (2017)  A consensus-based model for global optimization and its mean-field limit.  Mathematical Models and Methods in Applied Sciences, 27, 183–204.
  * Pulido et al. (2018) Pulido, M., Tandeo, P., Bocquet, M., Carrassi, A. & Lucini, M. (2018)  Stochastic parameterization identification using ensemble Kalman filtering combined with maximum likelihood methods.  Tellus A, 70, 1442099.
  * Rebeschini & Van Handel (2015) Rebeschini, P. & Van Handel, R. (2015)  Can local particle filters beat the curse of dimensionality?  The Annals of Applied Probability, 25, 2809–2866.
  * Reich (2011) Reich, S. (2011)  A dynamical systems framework for intermittent data assimilation.  BIT Numerical Mathematics, 51, 235–249.
  * Reich (2013) Reich, S. (2013)  A nonparametric ensemble transform method for Bayesian inference.  SIAM J. Sci. Comput., 35, A2013–A2024.
  * Reich (2019) Reich, S. (2019)  Data assimilation: the schrödinger perspective.  Acta Numerica, 28, 635–711.
  * Reich & Cotter (2015) Reich, S. & Cotter, C. (2015)  Probabilistic forecasting and Bayesian data assimilation.  Cambridge University Press.
  * Reich & Weissmann (2021) Reich, S. & Weissmann, S. (2021)  Fokker-Planck particle systems for Bayesian inference: Computational approaches.  SIAM/ASA J. Uncertain. Quantification, 9, 446–482.
  * Robinson et al. (2018) Robinson, G., Grooms, I. & Kleiber, W. (2018)  Improving particle filter performance by smoothing observations.  Monthly Weather Review, 146, 2433–2446.
  * Sakov et al. (2012) Sakov, P., Oliver, D. S. & Bertino, L. (2012)  An iterative EnKF for strongly nonlinear systems.  Mon. Wea. Rev., 140, 1988–2004.
  * Sakov & Oke (2008) Sakov, P. & Oke, P. (2008)  A deterministic formulation of the ensemble Kalman filter: An alternative to ensemble square root filters.  Tellus A, 60, 361–371.
  * Salman et al. (2006) Salman, H., Kuznetsov, L., Jones, C. & Ide, K. (2006)  A method for assimilating Lagrangian data into a shallow-water-equation ocean model.  Monthly Weather Review, 134, 1081–1101.
  * Sampson et al. (2021) Sampson, C., Carrassi, A., Aydoğdu, A. & Jones, C. K. (2021)  Ensemble Kalman filter for nonconservative moving mesh solvers with a joint physics and mesh location update.  Quarterly Journal of the Royal Meteorological Society, 147, 1539–1561.
  * Sanz-Alonso & Stuart (2015) Sanz-Alonso, D. & Stuart, A. M. (2015)  Long-time asymptotics of the filtering distribution for partially observed chaotic dynamical systems.  SIAM/ASA Journal on Uncertainty Quantification, 3, 1200–1220.
  * Särkkä (2013) Särkkä, S. (2013)  Bayesian Filtering and Smoothing.  Cambridge: Cambridge University Press.
  * Schillings & Stuart (2017) Schillings, C. & Stuart, A. M. (2017)  Analysis of the ensemble Kalman filter for inverse problems.  SIAM Journal on Numerical Analysis, 55, 1264–1290.
  * Schillings & Stuart (2018) Schillings, C. & Stuart, A. M. (2018)  Convergence analysis of ensemble Kalman inversion: The linear, noisy case.  Applicable Analysis, 97, 107–123.
  * Schneider et al. (2017) Schneider, T., Lan, S., Stuart, A. M. & Teixeira, J. (2017)  Earth system modeling 2.0: A blueprint for models that learn from observations and targeted high-resolution simulations.  Geophysical Research Letters, 44, 12–396.
  * Sisson et al. (2018) Sisson, S. A., Fan, Y. & Beaumont, M. (2018)  Handbook of approximate Bayesian computation.  CRC Press.
  * Snyder et al. (2008) Snyder, C., Bengtsson, T., Bickel, P. & Anderson, J. (2008)  Obstacles to high-dimensional particle filtering.  Monthly Weather Review, 136, 4629–4640.
  * Snyder (2014) Snyder, C. (2014)  Introduction to the Kalman filter.  Advanced Data Assimilation for Geosciences, Lecture Notes of the Les Houches School of Physics, June 2012.  Oxford University Press.
  * Spantini et al. (2019a) Spantini, A., Baptista, R. & Marzouk, Y. (2019a)  Coupling techniques for nonlinear ensemble filtering.  arXiv preprint arXiv:1907.00389.
  * Spantini et al. (2019b) Spantini, A., Baptista, R. & Marzouk, Y. (2019b)  Coupling techniques for nonlinear ensemble filtering.  arXiv preprint arXiv:1907.00389.
  * Stordal et al. (2011) Stordal, A., Karlsen, H., Nævdal, G., Skaug, H. & Vallés, B. (2011)  Bridging the ensemble Kalman filter and particle filters: the adaptive Gaussian mixture filter.  Comput. Geosci., 15, 293–305.
  * Strang (1968) Strang, G. (1968)  On the construction and comparison of difference schemes.  SIAM journal on numerical analysis, 5, 506–517.
  * Stuart (2010) Stuart, A. M. (2010)  Inverse problems: a bayesian perspective.  Acta numerica, 19, 451–559.
  * Taghvaei et al. (2017) Taghvaei, A., de Wiljes, J., Mehta, P. & Reich, S. (2017)  Kalman filter and its modern extensions for the continuous-time nonlinear filtering problem.  ASME. J. Dyn. Sys., Meas., Control., 140, 030904–030904–11.
  * Taghvaei et al. (2020) Taghvaei, A., Mehta, P. & Meyn, S. (2020)  Diffusion map-based algorithm for gain function approximation in the feedback particle filter.  SIAM/ASA J. Uncertain. Quantif., 8, 1090–1117.
  * Taghvaei & Hosseini (2022) Taghvaei, A. & Hosseini, B. (2022)  An optimal transport formulation of bayes’ law for nonlinear filtering algorithms.  arXiv preprint arXiv:2203.11869.
  * Taghvaei & Mehta (2020) Taghvaei, A. & Mehta, P. G. (2020)  An optimal transport formulation of the ensemble kalman filter.  IEEE Transactions on Automatic Control, 66, 3052–3067.
  * Tippett et al. (2003) Tippett, M., Anderson, J., Bishop, C., Hamill, T. & Whitaker, J. (2003)  Ensemble square root filters.  Mon. Weather Rev., 1485–1490.
  * Tödter & Ahrens (2015) Tödter, J. & Ahrens, B. (2015)  A second-order exact ensemble square root filter for nonlinear data assimilation.  Mon. Wea. Rev., 143, 1347–1367.
  * Tong et al. (2015) Tong, X. T., Majda, A. J. & Kelly, D. (2015)  Nonlinear stability of the ensemble Kalman filter with adaptive covariance inflation.  arXiv preprint arXiv:1507.08319.
  * Tong et al. (2016) Tong, X. T., Majda, A. J. & Kelly, D. (2016)  Nonlinear stability and ergodicity of ensemble based Kalman filters.  Nonlinearity, 29, 657.
  * Tsianos et al. (2012) Tsianos, K. I., Lawlor, S. & Rabbat, M. G. (2012)  Consensus-based distributed optimization: Practical issues and applications in large-scale machine learning.  2012 50th annual allerton conference on communication, control, and computing (allerton).  IEEE, pp. 1543–1550.
  * van Leeuwen et al. (2019) van Leeuwen, P., Künsch, H., Nerger, L., Potthast, R. & Reich, S. (2019)  Particle filters for high-dimensional geoscience applications: A review.  Q. J. Royal Meteorol. Soc., 145, 2335–2365.
  * van Leeuwen (2020) van Leeuwen, P. J. (2020)  A consistent interpretation of the stochastic version of the Ensemble Kalman filter.  Quarterly Journal of the Royal Meteorological Society, 146, 2815–2825.
  * Van Leeuwen & Evensen (1996) Van Leeuwen, P. J. & Evensen, G. (1996)  Data assimilation and inverse methods in terms of a probabilistic formulation.  Monthly Weather Review, 124, 2898–2913.
  * Vanden-Eijnden et al. (2003) Vanden-Eijnden, E. et al. (2003)  Fast communications: Numerical techniques for multi-scale dynamical systems with stochastic effects.  Communications in Mathematical Sciences, 1, 385–391.
  * Villani (2008) Villani, C. (2008)  Optimal transport: Old and new, vol. 338.  Springer Science & Business Media.
  * Villani (2021) Villani, C. (2021)  Topics in optimal transportation, vol. 58.  American Mathematical Soc.
  * Wang et al. (2004) Wang, X., Bishop, C. & Julier, S. (2004)  Which is better, an ensemble of positive-nagative pairs or a centered spherical simplex ensemble?  Monthly Weather Review, 132, 1590–1605.
  * Welch et al. (1995) Welch, G., Bishop, G. et al. (1995)  An introduction to the Kalman filter.
  * Whitaker & Hamill (2002) Whitaker, J. & Hamill, T. (2002)  Ensemble data assimilation without perturbed observations.  Monthly Weather Review, 130, 1913–1924.
  * Yang & Grooms (2021) Yang, L. M. & Grooms, I. (2021)  Machine learning techniques to construct patched analog ensembles for data assimilation.  arXiv preprint arXiv:2103.00318.
  * Yang et al. (2013) Yang, T., Mehta, P. G. & Meyn, S. P. (2013)  Feedback particle filter.  IEEE Trans. Automat. Control, 58, 2465–2480.
  * Yang et al. (2014) Yang, T., Blom, H. A. P. & Mehta, P. G. (2014)  The continuous-discrete time feedback particle filter.  IEEE, pp. 648–653.
  * Zech & Marzouk (2022) Zech, J. & Marzouk, Y. (2022)  Sparse approximation of triangular transports, part I: The finite-dimensional case.  Constructive Approximation, 55, 919–986.
  * Zhang et al. (2017) Zhang, C., Taghvaei, A. & Mehta, P. G. (2017)  A controlled particle filter for global optimization.  Technical Report.  arXiv:1701.02413.
  * Zupanski (2005) Zupanski, M. (2005)  Maximum likelihood ensemble filter: Theoretical aspects.  Mon. Weath. Rev., 133, 1710–1726.

##  6 Appendix A (Pseudo-Code)

In this appendix we provide pseudo-code describing several of the algorithms
that we present and deploy in this paper. Algorithms 1 and 2, 3DVAR and the
ensemble Kalman filter (EnKF) respectively, are presented 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. We refer to Algorithm
3 as Ensemble Kalman Inversion (Transport), as it arises from the approach to
inversion described in Subsection 4.2; we refer to Algorithm 4 as Ensemble
Kalman Inversion (Iteration to Infinity), as it arises from the approach to
inversion described in Subsection 4.4; Algorithm 5 corresponds to particular
cases of Algorithms 3 and 4 where they coincide as explained in Subsection
4.5.

Algorithm 1 3DVAR

Input: Initial $v_{0}\in\mathbb{R}^{d_{v}}$ and fixed gain matrix
$K\in\mathbb{R}^{d_{v}\times d_{y}}$.

for $n=0$ to $N-1$ do

Prediction: $\widehat{v}_{n+1}=\Psi(v_{n}).$

Analysis:
$v_{n+1}=\widehat{v}_{n+1}+K\Bigl{(}{y}^{\dagger}_{n+1}-h\bigr{(}\widehat{v}_{n+1}\bigr{)}\Bigr{)}.$

end for

Output: Estimates $\\{v_{n}\\}_{n=0}^{N}$.

Algorithm 2 EnKF

Input: Ensemble size $J$, initial ensemble $\\{v_{0}^{(j)}\\}_{j=1}^{J}$.

for $n=0$ to $N-1$ do

Prediction: for $j=1,\dots,J$ do  $\displaystyle\xi_{n}^{(j)}$
$\displaystyle\sim\text{N}(0,\Sigma),\quad\eta_{n}^{(j)}\sim\text{N}(0,\Gamma)$
$\displaystyle\widehat{v}_{n+1}^{(j)}$
$\displaystyle=\Psi(v_{n}^{(j)})+\xi_{n}^{(j)},$
$\displaystyle\widehat{y}_{n+1}^{(j)}$
$\displaystyle=h(\widehat{v}_{n+1}^{(j)})+\eta_{n}^{(j)}.$

Compute  $\displaystyle\widehat{m}_{n+1}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\widehat{v}_{n+1}^{(j)},\quad\widehat{o}_{n+1}=\frac{1}{J}\sum_{j=1}^{J}h(\widehat{v}_{n+1}^{(j)}),$
$\displaystyle\widehat{C}^{vh}_{n+1}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}(\widehat{v}_{n+1}^{(j)}-\widehat{m}_{n+1})\otimes(\widehat{v}_{n+1}^{(j)}-\widehat{o}_{n+1}),$
$\displaystyle\widehat{C}^{hh}_{n+1}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}(\widehat{v}_{n+1}^{(j)}-\widehat{o}_{n+1})\otimes(\widehat{v}_{n+1}^{(j)}-\widehat{o}_{n+1}),$
$\displaystyle K_{n+1}$
$\displaystyle=\widehat{C}_{n+1}^{vh}\left(\widehat{C}_{n+1}^{hh}+\Gamma\right)^{-1}.$

Analysis: for $j=1,\dots,J$ do
$v_{n+1}^{(j)}=\widehat{v}_{n+1}^{(j)}+K_{n+1}\Bigl{(}y_{n+1}^{(j)}-h\bigr{(}\widehat{v}_{n+1}^{(j)}\bigr{)}\Bigr{)}.$

end for

Output: Ensembles $\\{v_{n}^{(j)}\\}_{j=1}^{J}$ for $n=0,\dots,N$.

Algorithm 3 Ensemble Kalman Inversion (Transport)

Input: Data $w$, $N$ and $\Delta t$ such that $N\Delta t=1$, ensemble size
$J$, initial ensemble $\\{u_{0}^{(j)}\\}_{j=1}^{J}$.

for $n=0$ to $N-1$ do

Prediction: for $j=1,\dots,J$ do  $\displaystyle\eta_{n}^{(j)}$
$\displaystyle\sim\text{N}\Bigl{(}0,\frac{\Gamma}{\Delta t}\Bigr{)}$
$\displaystyle\widehat{u}_{n+1}^{(j)}$ $\displaystyle=u_{n}^{(j)}$
$\displaystyle\widehat{w}_{n+1}^{(j)}$ $\displaystyle=w+\eta_{n}^{(j)}.$

Compute  $\displaystyle\widehat{m}_{n+1}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\widehat{u}_{n+1}^{(j)},\quad\widehat{o}_{n+1}=\frac{1}{J}\sum_{j=1}^{J}G(\widehat{u}_{n+1}^{(j)}),$
$\displaystyle\widehat{C}^{uG}_{n+1}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\bigl{(}\widehat{u}_{n+1}^{(j)}-\widehat{m}_{n+1}\bigr{)}\otimes\bigl{(}G(\widehat{u}_{n+1}^{(j)})-\widehat{o}_{n+1}\bigr{)},$
$\displaystyle\widehat{C}^{GG}_{n+1}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\bigl{(}G(\widehat{u}_{n+1}^{(j)})-\widehat{o}_{n+1}\bigr{)}\otimes\bigl{(}G(\widehat{u}_{n+1}^{(j)})-\widehat{o}_{n+1}\bigr{)}.$

Analysis: for $j=1,\dots,J$ do  $u_{n+1}^{(j)}=\widehat{u}_{n+1}^{(j)}+\Delta
t\widehat{C}_{n+1}^{uG}\left(\Delta
t\widehat{C}_{n+1}^{GG}+\Gamma\right)^{-1}\bigl{(}\widehat{w}_{n+1}^{(j)}-G(\widehat{u}_{n+1}^{(j)})\bigr{)}.$

end for

Output: Ensembles $\\{u_{n}^{(j)}\\}_{j=1}^{J}$ for $n=0,\dots,N$.

Algorithm 4 Ensemble Kalman Inversion (Iteration to Infinity)

Input: Data $w$, $N$, $\alpha\in(0,1)$, $r_{0}$,
$\sigma^{\prime},\gamma^{\prime}$ ensemble size $J$, initial ensemble
$\\{u_{0}^{(j)}\\}_{j=1}^{J}$.

for $n=0$ to $N-1$ do

Prediction: for $j=1,\dots,J$ do  $\displaystyle\xi_{n}^{(j)}$
$\displaystyle\sim\text{N}(0,\sigma^{\prime}\mathsf{\Sigma}),\quad\eta_{n}^{(j)}\sim\text{N}(0,\gamma^{\prime}\mathsf{\Gamma})$
$\displaystyle\widehat{u}_{n+1}^{(j)}$ $\displaystyle=\alpha
u_{n}^{(j)}+(1-\alpha)r_{0}+\xi_{n}^{(j)}$
$\displaystyle\widehat{w}_{n+1}^{(j)}$ $\displaystyle=w+\eta_{n}^{(j)}.$

Compute  $\displaystyle\widehat{m}_{n+1}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\widehat{u}_{n+1}^{(j)},\quad\widehat{o}_{n+1}=\frac{1}{J}\sum_{j=1}^{J}G(\widehat{u}_{n+1}^{(j)}),$
$\displaystyle\widehat{C}^{uG}_{n+1}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\bigl{(}\widehat{u}_{n+1}^{(j)}-\widehat{m}_{n+1}\bigr{)}\otimes\bigl{(}G(\widehat{u}_{n+1}^{(j)})-\widehat{o}_{n+1}\bigr{)},$
$\displaystyle\widehat{C}^{GG}_{n+1}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\bigl{(}G(\widehat{u}_{n+1}^{(j)})-\widehat{o}_{n+1}\bigr{)}\otimes\bigl{(}G(\widehat{u}_{n+1}^{(j)})-\widehat{o}_{n+1}\bigr{)}.$

Analysis: for $j=1,\dots,J$ do
$u_{n+1}^{(j)}=\widehat{u}_{n+1}^{(j)}+\widehat{C}_{n+1}^{uG}\left(\widehat{C}_{n+1}^{GG}+\gamma^{\prime}\mathsf{\Gamma}\right)^{-1}\bigl{(}\widehat{w}_{n+1}^{(j)}-G(\widehat{u}_{n+1}^{(j)})\bigr{)}.$

end for

Output: Ensembles $\\{u_{n}^{(j)}\\}_{j=1}^{J}$ for $n=0,\dots,N$.

Algorithm 5 Ensemble Kalman Inversion (Transport and Iteration to Infinity)

Input: Data $w$, ensemble size $J$, initial ensemble
$\\{u_{0}^{(j)}\\}_{j=1}^{J}$.

for $n=0$ to $N-1$ do

Prediction: for $j=1,\dots,J$ do  $\displaystyle\eta_{n}^{(j)}$
$\displaystyle\sim\text{N}(0,\Gamma)$ $\displaystyle\widehat{u}_{n+1}^{(j)}$
$\displaystyle=u_{n}^{(j)}$ $\displaystyle\widehat{w}_{n+1}^{(j)}$
$\displaystyle=w+\eta_{n}^{(j)}.$

Compute  $\displaystyle\widehat{m}_{n+1}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\widehat{u}_{n+1}^{(j)},\quad\widehat{o}_{n+1}=\frac{1}{J}\sum_{j=1}^{J}G(\widehat{u}_{n+1}^{(j)}),$
$\displaystyle\widehat{C}^{uG}_{n+1}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\bigl{(}\widehat{u}_{n+1}^{(j)}-\widehat{m}_{n+1}\bigr{)}\otimes\bigl{(}G(\widehat{u}_{n+1}^{(j)})-\widehat{o}_{n+1}\bigr{)},$
$\displaystyle\widehat{C}^{GG}_{n+1}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\bigl{(}G(\widehat{u}_{n+1}^{(j)})-\widehat{o}_{n+1}\bigr{)}\otimes\bigl{(}G(\widehat{u}_{n+1}^{(j)})-\widehat{o}_{n+1}\bigr{)}.$

Analysis: for $j=1,\dots,J$ do
$u_{n+1}^{(j)}=\widehat{u}_{n+1}^{(j)}+\widehat{C}_{n+1}^{uG}\left(\widehat{C}_{n+1}^{GG}+\Gamma\right)^{-1}\bigl{(}\widehat{w}_{n+1}^{(j)}-G(\widehat{u}_{n+1}^{(j)})\bigr{)}.$

end for

Output: Ensembles $\\{u_{n}^{(j)}\\}_{j=1}^{J}$ for $n=0,\dots,N$.

##  7 Appendix B (Lorenz ’96 Multiscale)

To illustrate the problems of both state estimation and parameter estimation
we use, throughout this paper, variants on the Lorenz ’96 model. In particular
we use both the Lorenz ’96 multiscale system and a singlescale closure derived
from it in the scale-separated case; see Subsection 2.7, and the developments
in this section, for bibliographic details. If we generate data with the
singlescale model, and assimilate using this model, we are able to test
algorithms in their basic (perfect model) form. If we generate data with the
multiscale model, and assimilate using the singlescale model, this will also
allow us to study the effect of model misspecification on data assimilation.

###  7.1 Lorenz ’96 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$, the multiscale model is defined by equations (11) and (12); the
boundary conditions (13) complete definition of the equations and provide
further couplings beyond those in the differential equations themselves. 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.

###  7.2 Lorenz ’96 Singlescale Model

If $\epsilon\ll 1$ then the dynamics for the $w=\\{w_{\ell,j}\\}$ governed
(11b) evolve on a much faster timescale than the dynamics for the
$v=\\{v_{\ell}\\}$ governed by (11a). Thus it is a reasonable approximation to
think of $v$ as frozen in (11b). In the following let
$\overline{w}=\\{\overline{w}_{\ell}\\}.$ If we assume that the dynamics of
$w$ with $v$ frozen is ergodic with invariant measure $\mu^{v}(dw)$ (a measure
in $w$, parameterized by $v$) then the averaging principle Vanden-Eijnden et
al. (2003); Abdulle et al. (2012); Pavliotis & Stuart (2008) suggests that we
may make the approximation

$\overline{w}\approx M(v):=\int\overline{w}\,\mu^{v}(dw).$

If we also invoke the approximation $M_{\ell}(v)\approx m(v_{\ell})$, which is
shown to be valid for large $J$ in Fatkullin & Vanden-Eijnden (2004), then we
arrive at the singlescale Lorenz ’96 model (6). This program of analysis for
the Lorenz ’96 model, and studies of the validity of the resulting singlescale
approximation, was established in the paper Fatkullin & Vanden-Eijnden (2004).
The function $m(\cdot)$ is not given explicitly, but may be fit to data in
various different ways, as explained in Fatkullin & Vanden-Eijnden (2004).
Figure 2 shows such an $m$, fit using Gaussian process regression methodology
as detailed in Subsection 4.3 of Burov et al. (2021).

##  8 Appendix C (Mean Field Maps)

In Subsection 8.1 we discuss the existence, form and properties of mean field
maps which carry out the program of approximate transport described in
Subsection 2.5.3; these maps require access to simulated data and are hence
referred as as stochastic second order transport maps. In Subsection 8.2 we
discuss the existence, form and properties of mean field maps which carry out
the program of approximate transport described in Subsection 2.5.4; these maps
do not require access to simulated data and are hence referred as
deterministic second order transport maps. The two approaches provide
fundamental underpinnings of ensemble Kalman methods as we develop them here,
but were not adopted in its historical development. Subsection 8.3 is devoted
to the minimum variance approximation, a way of deriving the Kalman transport
map (46), which is part of the historical development of the subject, but
which does not play a central role in our presentation and analysis of the
subject.

We note that the maps of interest in Subsection 8.1 take
$\textit{Law}\,(\widehat{v}_{n+1},\widehat{y}_{n+1})$ into
$\textit{Law}\,(v_{n+1})$ and hence we may drop the subscript $n+1$ throughout
the analysis. Similar considerations apply in Subsection 8.2 and in Subsection
8.3.

###  8.1 Mean Field Maps – Simulated Data

Let $(\widehat{v},\widehat{y})\sim\nu$ where

$G\nu=\mathsf{N}\Bigl{(}\begin{bmatrix}\widehat{m}\\\
\widehat{o}\end{bmatrix},\begin{bmatrix}\widehat{C}&\widehat{C}^{vy}\\\
(\widehat{C}^{vy})^{\top}&\widehat{C}^{yy},\end{bmatrix}\Bigr{)}.$

Define

(183a) $\displaystyle m$
$\displaystyle=\widehat{m}+\widehat{C}^{vy}(\widehat{C}^{yy})^{-1}\bigl{(}{y}^{\dagger}-\widehat{o}\bigr{)},$
(183b) $\displaystyle C$
$\displaystyle=\widehat{C}-\widehat{C}^{vy}(\widehat{C}^{yy})^{-1}(\widehat{C}^{vy})^{\top}.$

Note that these are the mean and covariance of the Gaussian random variable
$\widehat{v}$ conditioned on $\widehat{y}={y}^{\dagger}$ under Gaussian
measure $G\nu;$ see (30). The quantities in the above identities are as
defined in Subsection 2.3.2, with the exception that the subscripts $n+1$ have
been dropped for clarity of exposition.

Our goal is to identify maps of the form

(184) $v=A\widehat{v}+B\widehat{y}+a$

so that, if $(\widehat{v},\widehat{y})\sim\nu$, then $v$ has mean $m$ and
covariance $C$ given by (183). We make the following assumptions on the
covariance under $G\nu$ and on the matrices $A,B$ and vector $a:$

######  Assumptions 8.1.

The covariance under $G\nu$ is invertible. The matrices $A,B$ and vector $a$
may depend on ${y}^{\dagger}$ and measure $\nu$ but not on the random
variables $(\widehat{v},\widehat{y}).\quad\blacksquare$

Recall the discussion of pushforward of measures in the introduction to
Subsection 2.5. Under Assumptions 8.1 pushforward under the map (184), when
chosen to match the desired first and second order statistics, defines a
nonlinear map on the space of measures, and in particular on $\nu$ itself. In
what follows all expectations are computed under $\nu.$ Since the covariance
under $G\nu$ is strictly positive-definite so are the marginal covariances
$\widehat{C}$ and $\widehat{C}^{yy}$ (see Stuart (2010)[Lemma 6.21]). Thus we
may define the conditional mean and covariance by (183), as well as Kalman
gain $K$, and conditional covariance $\widetilde{C}$, given by

$\displaystyle K$ $\displaystyle=\widehat{C}^{vy}(\widehat{C}^{yy})^{-1},$
$\displaystyle\widetilde{C}$
$\displaystyle=\widehat{C}^{yy}-(\widehat{C}^{vy})^{\top}(\widehat{C})^{-1}\widehat{C}^{vy};$

we note that the conditional covariances $C$ and $\widetilde{C}$ are also
strictly positive-definite (see Stuart (2010)[Lemma 6.21]). From equations
(183a) and (184) the following is immediate:

######  Lemma 8.2.

Let Assumptions 8.1 hold and let $(\widehat{v},\widehat{y})\sim\nu.$ If $v$
given by (184) has mean given by (183a) then

$a=(I-A)\mathbb{E}\widehat{v}+K{y}^{\dagger}-(B+K)\mathbb{E}\widehat{y}.$

As a consequence it follows that

(185)
$v=\mathbb{E}\widehat{v}+A(\widehat{v}-\mathbb{E}\widehat{v})+B(\widehat{y}-\mathbb{E}\widehat{y})+K({y}^{\dagger}-\mathbb{E}\widehat{y}).$

and that

$\mathbb{E}\Bigl{(}(v-\mathbb{E}v)\otimes(v-\mathbb{E}v)\Bigr{)}=A\widehat{C}A^{\top}+B\widehat{C}^{yy}B^{\top}+A\widehat{C}^{vy}B^{\top}+B(\widehat{C}^{vy})^{\top}A.$

Thus, to match the covariance of the conditioned random variable, we obtain

(186)
$C=A\widehat{C}A^{\top}+B\widehat{C}^{yy}B^{\top}+A\widehat{C}^{vy}B^{\top}+B(\widehat{C}^{vy})^{\top}A.$

######  Theorem 8.3.

Let Assumptions 8.1 hold and let $(\widehat{v},\widehat{y})\sim\nu.$ If $a$ is
given by Lemma 8.2 then $v$ defined by (184) has covariance (183b) if and only
if real-valued matrices $A$ and $B$ are related by the identity

(187) $F\widehat{C}^{-1}F^{\top}=C-B\widetilde{C}B^{\top},$

where

(188) $F=A\widehat{C}+B(\widehat{C}^{vy})^{\top}.$

######  Proof 8.4.

We complete the square on the right hand side of (186) to obtain

$(A\widehat{C}+B(\widehat{C}^{vy})^{\top})\widehat{C}^{-1}(A\widehat{C}+B(\widehat{C}^{vy})^{\top})^{\top}=C^{\prime}.$

where

$C^{\prime}=C-B\widetilde{C}B^{\top}.$

Rearranging gives the desired result.

Define

$\mathcal{B}=\\{B\in\mathbb{R}^{d_{v}\times d_{y}}:C^{\prime}\succ 0\\},$

noting that this set is non-empty, and contains an open (and hence
uncountable) set of $B$, since $C\succ 0$. For $B\in\mathcal{B}$ consider the
eigenvalue problem

$\displaystyle C^{\prime}\varphi^{(i)}$
$\displaystyle=(s^{(i)})^{2}\varphi^{(i)},$
$\displaystyle\langle\varphi^{(i)},\varphi^{(j)}\rangle$
$\displaystyle=\delta_{ij}.$

Note that this has $d_{v}$ real solutions, up to sign changes in the
eigenvectors and assuming the $s^{(i)}$ to be non-negative. We now seek to
express $F$ in terms of $B\in\mathcal{B}.$ Writing the SVD
$F\widehat{C}^{-\frac{1}{2}}=U\Sigma V^{\top}$, where
$U,V\in\mathbb{R}^{d_{v}\times d_{v}}$ are orthogonal matrices and
$\Sigma\in\mathbb{R}^{d_{v}\times d_{v}}$ is a diagonal matrix, we see from
(187) that

$U\Sigma^{2}U^{\top}=C^{\prime}$

so that $U$ has columns given by the $\\{\varphi^{(i)}\\}_{i=1}^{d_{v}}$ and
corresponding diagonal entries of $\Sigma$, $\pm s^{(i)}$. We define

(189)
$U=\bigl{(}\varphi^{(1)},\cdots,\varphi^{(d_{v})}\bigr{)},\quad\Sigma={\rm
diag}\bigl{(}\pm s^{(1)},\cdots,\pm s^{(d_{v})}\bigr{)}.$

######  Corollary 8.5.

For every $B\in\mathcal{B}$ the choices of $A$ such that the pair $(A,B)$
satisfies the criterion of Theorem 8.3 are defined as follows. For $U,\Sigma$
as given in (189), set

$F=U\Sigma V^{\top}\widehat{C}^{\frac{1}{2}}$

where $V$ is an arbitrary orthogonal matrix in $\mathbb{R}^{d_{v}\times
d_{v}}.$ Then

$A=\Bigl{(}F-B(\widehat{C}^{vy})^{\top}\Bigr{)}\widehat{C}^{-1}.$

######  Example 8.6.

Among the uncountably many possible solutions for pairs $(A,B)$ we highlight
two. The choice $B=0$ is interesting because it does not require the data
variable $\widehat{y}$ in the definition of $v$. The choice $A=I$ is
interesting because it does not require action of an operator acting on
$\widehat{v}.$

The first, with $B=0$, allows the choice
$F=C^{\frac{1}{2}}\widehat{C}^{\frac{1}{2}}$ and then
$A=C^{\frac{1}{2}}\widehat{C}^{-\frac{1}{2}}.$ Thus the map (185) becomes

$\displaystyle v$
$\displaystyle=\mathbb{E}\widehat{v}+C^{\frac{1}{2}}\widehat{C}^{-\frac{1}{2}}(\widehat{v}-\mathbb{E}\widehat{v})+\widehat{C}^{vy}(\widehat{C}^{yy})^{-1}({y}^{\dagger}-\mathbb{E}\widehat{y})$
$\displaystyle=m+C^{\frac{1}{2}}\widehat{C}^{-\frac{1}{2}}(\widehat{v}-\mathbb{E}\widehat{v}).$

The second comes from setting $B=-K$ which leads to the possible choice $F=C$
and $A=I$ under which the map (185) becomes

$v=\widehat{v}+\widehat{C}^{vy}(\widehat{C}^{yy})^{-1}({y}^{\dagger}-\widehat{y}).$

We refer to this as the Kalman transport solution. $\blacksquare$

Given the plethora of solutions for matrices $(A,B)$, all of which effect the
desired measure transport from $\nu$ into the conditional, it is natural to
ask how to choose a specific pair $(A,B)$. One possibility is to use optimal
transport. To this end we define, for positive definite
$W\in\mathbb{R}^{d_{v}\times d_{v}}$,

$I_{W}(A,B)=\frac{1}{2}\mathbb{E}\langle(v-\widehat{v}),W(v-\widehat{v})\rangle.$

######  Theorem 8.7.

Let $v$ be given by (185). Consider the problem of finding minimizers of
$I_{W}(A,B)$ over pairs $(A,B)$ satisfying (187), (188); we refer to the
resulting map evaluated at such an $(A,B)$ as an optimal transport in the
$W-$weighted Euclidean distance. For any positive definite $W$, such
minimizers satisfy $B=0.$ In particular the Kalman transport solution is not
an optimal transport solution.

######  Proof 8.8.

We formulate the optimization problem over the pair $(F,B)$ since, because
$\widehat{C}$ is invertible, there is a bijection between $(A,B)$ and $(F,B).$
Let $:$ denotes the Frobenius inner-product. Then, under constraint (188),

(190a) $\displaystyle I_{W}(A,B)$ $\displaystyle=J_{W}(F)+{\rm const},$ (190b)
$\displaystyle J_{W}(F)$ $\displaystyle=-F:W,$

where ${\rm const}$ denotes a matrix independent of $A,B,F$ (with exact value
changing from instance to instance.) To see this note that

$v-\widehat{v}=-(\widehat{v}-\mathbb{E}\widehat{v})+A(\widehat{v}-\mathbb{E}\widehat{v})+B(\widehat{y}-\mathbb{E}\widehat{y})+K({y}^{\dagger}-\mathbb{E}\widehat{y}).$

Now, using (186),

$\mathbb{E}\Bigl{(}(v-\widehat{v})\otimes(v-\widehat{v})\Bigr{)}=\widehat{C}-\widehat{C}A^{\top}-A\widehat{C}-\widehat{C}^{vy}B^{\top}-B(\widehat{C}^{vy})^{\top}+C+{\rm
const.}$

Noting that

$I_{W}(A,B)=\frac{1}{2}\mathbb{E}\Bigl{(}(v-\widehat{v})\otimes(v-\widehat{v})\Bigr{)}:W,$

and that $D:W=D^{\top}:W$ for all $D$ since $W$ is symmetric, identity (190)
follows from (188). It then also follows that minimization of $I_{W}(A,B)$
subject to the constraints given by (187), (188) is equivalent to minimization
of $J_{W}(F)$ subject to the constraint (187).

To effect this latter minimization we introduce the Lagrange multiplier
$L\in\mathbb{R}^{d_{v}\times d_{v}}$, symmetric because the constraint is
symmetric, and define

$\tilde{\mathsf{J}}_{W}(F,B,L)=-F:W+L:\bigl{(}F\widehat{C}^{-1}F^{\top}+B\widetilde{C}B^{\top}-C\bigr{)}.$

Differentiating with respect to $F,B$ and $L$ respectively gives

(191a) $\displaystyle-W+2LF\widehat{C}^{-1}$ $\displaystyle=0,$ (191b)
$\displaystyle 2LB\widetilde{C}$ $\displaystyle=0,$ (191c) $\displaystyle
F\widehat{C}^{-1}F^{\top}+B\widetilde{C}B^{\top}$ $\displaystyle=C.$

Since $F$, $W$, and $\widehat{C}$ in (191a) are all invertible, the Lagrange
multiplier $L$ is necessarily invertible. Furthermore, since $\widetilde{C}$
is invertible because $\Sigma$ is, it follows from (191b) that $B=0$ as
required.

######  Example 8.9.

Define matrix $S$, and from it $A$, by

$S=\Bigl{(}C^{\frac{1}{2}}\widehat{C}C^{\frac{1}{2}}\Bigr{)}^{-\frac{1}{2}},\quad
A=C^{\frac{1}{2}}SC^{\frac{1}{2}}.$

We notice that if $B=0$, as required for an optimal transport solution, then
equation (191c) has as a solution $F=A\widehat{C}$ and then (191a) delivers
the Lagrange multiplier $L$. Note further that the solution is independent of
the specific choice of positive-definite $W$. $\blacksquare$

###  8.2 Mean Field Maps – No Simulated Data

Let $\widehat{v}\sim\widehat{\mu}$ and assume that
$\widehat{y}=h(\widehat{v})+\eta$ where $\eta\sim\mathsf{N}(0,\Gamma)$ is
independent of $\widehat{v}.$ Implicitly we have defined the joint
distribution $\nu$ of $(\widehat{v},\widehat{y}).$ We note that then,
expressed in terms of $\widehat{h}:=h(\widehat{v})$ and quantities defined in
(8.1),

$\displaystyle\widehat{C}^{vh}:=$
$\displaystyle\mathbb{E}(\widehat{v}-\mathbb{E}\widehat{v})\otimes(\widehat{h}-\mathbb{E}\widehat{h})=\widehat{C}^{vy},$
$\displaystyle\widehat{C}^{hh}:=$
$\displaystyle\mathbb{E}(\widehat{h}-\mathbb{E}\widehat{h})\otimes(\widehat{h}-\mathbb{E}\widehat{h})=\widehat{C}^{yy}-\Gamma.$

Thus we may rewrite (183) as

(192a) $\displaystyle m$
$\displaystyle=\widehat{m}+\widehat{C}^{vh}(\widehat{C}^{hh}+\Gamma)^{-1}\bigl{(}{y}^{\dagger}-\mathbb{E}\widehat{h}\bigr{)},$
(192b) $\displaystyle C$
$\displaystyle=\widehat{C}-\widehat{C}^{vh}(\widehat{C}^{hh}+\Gamma)^{-1}(\widehat{C}^{vh})^{\top},$

We have now eliminated reference to $\widehat{y}$ and our goal is to identify
maps of the form

(193) $v=R\widehat{v}+S\widehat{h}+r$

so that, if $\widehat{v}\sim\widehat{\mu}$, then $v$ has mean $m$ and
covariance $C$ given by (192). We make the following assumptions on the
covariance under $G\nu$ and on the matrices $R,S$ and vector $r:$

######  Assumptions 8.10.

The covariance under $G\nu$ is invertible. The matrices $R,S$ and vector $r$
may depend on ${y}^{\dagger}$ and measure $\widehat{\mu}$ but not on the
random variable $(\widehat{v},\widehat{h}).\quad\blacksquare$

With these assumptions the pushforward under map (184), when constrained to
match the desired first and second order statistics, defines a nonlinear map
on the space of measures, and in particular on measure $\widehat{\mu}$. In
what follows all expectations are computed under $\widehat{\mu}.$ We note that
matrix $\widehat{C}$ is invertible and that $K$ in (8.1) may be rewritten as

(194) $K=\widehat{C}^{vh}(\widehat{C}^{hh}+\Gamma)^{-1}.\\\ $

From equations (192a) and (193) the following is immediate:

######  Lemma 8.11.

Let Assumptions 8.10 hold and let $\widehat{v}\sim\widehat{\mu}.$ If $v$ given
by (193) has mean given by (192a) then

$r=(I-R)\mathbb{E}\widehat{v}+K{y}^{\dagger}-(S+K)\mathbb{E}\widehat{h}.$

As a consequence it follows that

(195)
$v=\mathbb{E}\widehat{v}+R(\widehat{v}-\mathbb{E}\widehat{v})+S(\widehat{h}-\mathbb{E}\widehat{h})+K({y}^{\dagger}-\mathbb{E}\widehat{h}).$

and that

(196)
$\mathbb{E}\Bigl{(}(v-\mathbb{E}v)\otimes(v-\mathbb{E}v)\Bigr{)}=R\widehat{C}R^{\top}+S\widehat{C}^{hh}S^{\top}+R\widehat{C}^{vh}S^{\top}+S(\widehat{C}^{vh})^{\top}R.$

Thus, to match the covariance of the conditioned random variable, we obtain

(197)
$C=R\widehat{C}R^{\top}+S\widehat{C}^{hh}S^{\top}+R\widehat{C}^{vh}S^{\top}+S(\widehat{C}^{vh})^{\top}R.$

Define

$\check{C}=\widehat{C}^{hh}-(\widehat{C}^{vh})^{\top}(\widehat{C})^{-1}\widehat{C}^{vh}.$

######  Theorem 8.12.

Let Assumptions 8.10 hold and let $\widehat{v}\sim\widehat{\mu}.$ If $s$ is
given by Lemma 8.11 then $v$ defined by (193) has covariance (192b) if and
only if real-valued matrices $R$ and $S$ are related by the identity

$E\widehat{C}^{-1}E^{\top}=C-S\check{C}S^{\top},$

where

$E=R\widehat{C}+S(\widehat{C}^{vh})^{\top}.$

######  Proof 8.13.

We complete the square on the right hand side of (197) to obtain

$(R\widehat{C}+S(\widehat{C}^{vh})^{\top})\widehat{C}^{-1}(R\widehat{C}+S(\widehat{C}^{vh})^{\top})^{\top}+S\check{C}S^{\top}=C.$

Rearranging gives the desired result.

Define

$\mathcal{S}=\\{S\in\mathbb{R}^{d_{v}\times d_{y}}:C-S\check{C}S^{\top}\succ
0\\}$

and, for $S\in\mathcal{S}$ consider the eigenvalue problem

$\displaystyle(C-S\check{C}S^{\top})\psi^{(i)}$
$\displaystyle=(o^{(i)})^{2}\psi^{(i)},$
$\displaystyle\langle\psi^{(i)},\psi^{(j)}\rangle$
$\displaystyle=\delta_{ij}.$

noting that this has $d_{v}$ real solutions, upto sign changes in the
eigenvectors and assuming the $o^{(i)}$ to be non-negative. We now seek to
express $E$ in terms of $S\in\mathcal{S}.$ Writing the SVD
$E\widehat{C}^{-\frac{1}{2}}=W\Omega Z^{\top}$, where
$W,Z\in\mathbb{R}^{d_{v}\times d_{v}}$ are orthogonal matrices and
$\Omega\in\mathbb{R}^{d_{v}\times d_{v}}$ is a diagonal matrix, we see from
(187) that

$W\Omega^{2}W^{\top}=C-S\check{C}S^{\top}$

so that $W$ has columns given by the $\\{\psi^{(i)}\\}_{i=1}^{d_{v}}$ and
corresponding diagonal entries of $\Omega$, $\pm o^{(i)}$. We define

(199) $W=\bigl{(}\psi^{(1)},\cdots,\psi^{(d_{v})}\bigr{)},\quad\Omega={\rm
diag}(\pm o^{(1)},\dots,\pm o^{(d_{v})}.$

######  Corollary 8.14.

For every $S\in\mathcal{S}$ the choices of $R$ such that the pair $(R,S)$
satisfies the criterion of Theorem 8.12 are defined as follows. For $W,\Omega$
as given in (199), set

$E=W\Omega Z^{\top}\widehat{C}^{\frac{1}{2}}$

where $Z$ is an arbitrary orthogonal matrix in $\mathbb{R}^{d_{v}\times
d_{v}}.$ Then

$R=\Bigl{(}E-S(\widehat{C}^{vh})^{\top}\Bigr{)}\widehat{C}^{-1}.$

######  Example 8.15.

As in Example 8.6 we highlight two examples, here corresponding to $S=0$ and
to $R=I.$ The first, with $S=0$, allows the choice
$E=C^{\frac{1}{2}}\widehat{C}^{\frac{1}{2}}$ and hence
$R=C^{\frac{1}{2}}\widehat{C}^{-\frac{1}{2}}.$ We thus obtain (8.6) again.

The second comes from setting $R=I.$ This leads from (196) to the following
equation for $S:$

$S\widehat{C}^{hh}S^{\top}+\widehat{C}^{vh}S^{\top}+S(\widehat{C}^{vh})^{\top}=-\widehat{C}^{vh}(\widehat{C}^{hh}+\Gamma)^{-1}(\widehat{C}^{vh})^{\top}.$

We seek a solution for $S$ in the form

$S=-\widehat{C}^{vh}Y^{-1}$

and determine $Y$. We obtain the equation

$Y^{-1}\widehat{C}^{hh}Y^{-T}-Y^{-T}-Y^{-1}=-(\widehat{C}^{hh}+\Gamma)^{-1}.$

Thus, premultiplying by $Y$ and post multiplying by $Y^{\top}$, we obtain

$Y(\widehat{C}^{hh}+\Gamma)^{-1}Y^{\top}-Y-Y^{\top}+\widehat{C}^{hh}=0$

which factorizes to give

$\Bigl{(}Y(\widehat{C}^{hh}+\Gamma)^{-1}-I\Bigr{)}(\widehat{C}^{hh}+\Gamma)\Bigl{(}Y(\widehat{C}^{hh}+\Gamma)^{-1}-I\Bigr{)}^{\top}=\Gamma.$

Thus, taking the symmetric square root, we have

$\Bigl{(}Y(\widehat{C}^{hh}+\Gamma)^{-1}-I\Bigr{)}(\widehat{C}^{hh}+\Gamma)^{\frac{1}{2}}=\Gamma^{\frac{1}{2}}.$

Hence

$\Bigl{(}Y(\widehat{C}^{hh}+\Gamma)^{-1}-I\Bigr{)}=\Gamma^{\frac{1}{2}}(\widehat{C}^{hh}+\Gamma)^{-\frac{1}{2}}.$

Rearranging gives

$Y(\widehat{C}^{hh}+\Gamma)^{-1}=\Gamma^{\frac{1}{2}}(\widehat{C}^{hh}+\Gamma)^{-\frac{1}{2}}+(\widehat{C}^{hh}+\Gamma)^{\frac{1}{2}}(\widehat{C}^{hh}+\Gamma)^{-\frac{1}{2}}.$

Thus

$Y=\Bigl{(}\Gamma^{\frac{1}{2}}+(\widehat{C}^{hh}+\Gamma)^{\frac{1}{2}}\Bigr{)}(\widehat{C}^{hh}+\Gamma)^{\frac{1}{2}}.$

We obtain $S=-\widetilde{K}$ where

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

The map (195) becomes

$\displaystyle v$
$\displaystyle=\widehat{v}-\widetilde{K}(\widehat{h}-\mathbb{E}\widehat{h})+\widehat{C}^{vh}(\widehat{C}^{hh}+\Gamma)^{-1}({y}^{\dagger}-\mathbb{E}\widehat{h}),$
$\displaystyle=\widehat{v}-\widetilde{K}(\widehat{h}-\mathbb{E}\widehat{h})+K({y}^{\dagger}-\mathbb{E}\widehat{h}),$
$\displaystyle=m+(\widehat{v}-\widehat{m})-\widetilde{K}(\widehat{h}-\mathbb{E}\widehat{h}).$

###  8.3 Minimum Variance Approximation

In the two preceding subsections we have identified an uncountable set of
transport maps that match the second order statistics of true transport. Among
all these, the _Kalman transport_ (46) has a particular appeal because it is
constructed around the _Kalman gain_ familiar from filtering in the linear
Gaussian setting. In this subsection we show how the principle of minimizing
the variance within a class of _linear_ estimators of the state, given
observation, leads to this choice of transport map. We believe that the second
order transport approach, which we highlight in the main text, provides a more
fundamental viewpoint on the mean-field models which underpin ensemble Kalman
methods; however the minimum variance perspective is widely adopted in the
statistics and geophysics communities (see Subsection 2.7) and hence has an
important place in the subject of ensemble Kalman methods.

Recall that the Kalman transport approach leads back to the map (14c),
motivated by control theoretic considerations, and identifies a specific
choice of Kalman gain $K_{n}$, a choice which depends on the law of the
predicted state and data. In the Gaussian case the Kalman transport map
exactly generates the desired transport, a fact which we establish in Theorem
2.11. The general form of approximate stochastic second order transport maps
which we study using the second order transport approach in the main text is
(42). Our goal here is to motivate a specific choice of $\widetilde{T}^{S}$ in
(42c). In this subsection we achieve this by first defining, and identifying,
the _best linear unbiased estimator_ , BLUE for short. From this we will
derive Kalman transport.

######  Lemma 8.16.

Assume that $\Gamma\succ 0.$ Let all expectations be computed under
$\textit{Law}\,(\widehat{v}_{n+1},\widehat{y}_{n+1})$ and consider
$m^{\tiny{\rm BL}}_{n+1}$ in the form

$m^{\tiny{\rm BL}}_{n+1}=a^{\prime}+B\widehat{y}_{n+1}.$

Define $\widehat{C}_{n+1},\widehat{C}^{vy}_{n+1}$ and $\widehat{C}^{yy}_{n+1}$
by (23), (24) and

$\mathsf{I}(a^{\prime},B):=\mathbb{E}|\widehat{v}_{n+1}-m^{\tiny{\rm
BL}}_{n+1}|^{2}.$

Then $m^{\tiny{\rm BL}}_{n+1}$ is said to be the BLUE of $\widehat{v}_{n+1}$
given $\widehat{y}_{n+1}$, if vector $a^{\prime}$ and matrix $B$ are
independent of $(\widehat{v}_{n+1},\widehat{y}_{n+1})$, but may depend on
$\textit{Law}\,(\widehat{v}_{n+1},\widehat{y}_{n+1})$, and are chosen to
minimize $\mathsf{I}(a^{\prime},B).$ Then

(200) $m^{\tiny{\rm
BL}}_{n+1}=\mathbb{E}\widehat{v}_{n+1}+\widehat{C}^{vy}_{n+1}(\widehat{C}^{yy}_{n+1})^{-1}(\widehat{y}_{n+1}-\mathbb{E}\widehat{y}_{n+1}).$

Furthermore, the estimator (200) is unbiased, that is,

$\mathbb{E}(\widehat{v}_{n+1}-m^{\tiny{\rm BL}}_{n+1})=0,$

and its covariance with respect to
$\textit{Law}\,(\widehat{v}_{n+1},\widehat{y}_{n+1})$,

$C^{\tiny{\rm BL}}_{n+1}:=\mathbb{E}\bigl{(}(\widehat{v}_{n+1}-m^{\tiny{\rm
BL}}_{n+1})(\widehat{v}_{n+1}-m^{\tiny{\rm BL}}_{n+1})^{\top}\bigr{)},$

is given by

(201) $\displaystyle C^{\tiny{\rm BL}}_{n+1}$
$\displaystyle=\widehat{C}_{n+1}-\widehat{C}_{n+1}^{vy}(\widehat{C}_{n+1}^{yy})^{-1}({\widehat{C}_{n+1}^{vy}})^{\top}.$

######  Proof 8.17.

Noting that we may, without loss of generality, reparametrize the proposed
form of $m^{\tiny{\rm BL}}_{n+1}$ as

$m^{\tiny{\rm
BL}}_{n+1}=\mathbb{E}\widehat{v}_{n+1}+a+B(\widehat{y}_{n+1}-\mathbb{E}\widehat{y}_{n+1}),$

we see that the desired objective to be minimized over vector-matrix pair
$(a,B)$ is

$\displaystyle\mathsf{J}(a,B)$
$\displaystyle:=\frac{1}{2}\mathbb{E}|\widehat{v}_{n+1}-\mathbb{E}\widehat{v}_{n+1}-a-B(\widehat{y}_{n+1}-\mathbb{E}\widehat{y}_{n+1})|^{2}$
$\displaystyle=\frac{1}{2}\mathbb{E}|\widehat{v}_{n+1}-\mathbb{E}\widehat{v}_{n+1}|^{2}+\frac{1}{2}|a|^{2}$
$\displaystyle\quad\quad\quad+\frac{1}{2}\mathbb{E}|B(\widehat{y}_{n+1}-\mathbb{E}\widehat{y}_{n+1})|^{2}-\mathbb{E}\langle\widehat{v}_{n+1}-\mathbb{E}\widehat{v}_{n+1},B(\widehat{y}_{n+1}-\mathbb{E}\widehat{y}_{n+1})\rangle$
$\displaystyle=\frac{1}{2}\mathbb{E}|\widehat{v}_{n+1}-\mathbb{E}\widehat{v}_{n+1}|^{2}+\frac{1}{2}|a|^{2}+\frac{1}{2}(BB^{\top}):\widehat{C}_{n+1}^{yy}-B:\widehat{C}_{n+1}^{vy}.$

Clearly the minimizer with respect to $a$ is achieved by setting $a=0$.
Differentiating with respect to $B$, noting that $\widehat{C}_{n+1}^{yy}$ is
symmetric, shows that $B=K_{n}$ given by (25). That the resulting pair is
indeed a minimizer, and not another critical point, follows from the fact that
$\widehat{C}_{n+1}^{yy}$ is positive-definite; this is implied by the
assumption $\Gamma\succ 0.$ Thus we obtain the estimator (200). A
straightforward calculation shows that the estimator is unbiased and that its
covariance is given by (201).

We now connect the BLUE to an approximate (second order) transport map. To
this end, note that (200) may be viewed as mapping $\widehat{y}_{n+1}$ into
$m^{\tiny{\rm BL}}_{n+1}=m^{\tiny{\rm BL}}(\widehat{y}_{n+1}).$ With this
notation we define $m_{n+1}^{\dagger}$ by

(202a) $\displaystyle m_{n+1}^{\dagger}$ $\displaystyle=m^{\tiny{\rm
BL}}({y}^{\dagger}_{n+1})$ (202b)
$\displaystyle=\mathbb{E}\widehat{v}_{n+1}+\widehat{C}^{vy}_{n+1}(\widehat{C}^{yy}_{n+1})^{-1}({y}^{\dagger}_{n+1}-\mathbb{E}\widehat{y}_{n+1}).$

We may now make the following connection between BLUE and Kalman transport:

######  Theorem 8.18.

Let $(\widehat{v}_{n+1},\widehat{y}_{n+1})$ be distributed according to
measure $\nu_{n+1}$. Then the transport map (46c) has the properties:

$\displaystyle\mathbb{E}v_{n+1}$ $\displaystyle=m_{n+1}^{\dagger}$
$\displaystyle\mathbb{E}\bigl{(}(v_{n+1}-m_{n+1}^{\dagger})(v_{n+1}-m_{n+1}^{\dagger})^{\top}\bigr{)}$
$\displaystyle=C^{\tiny{\rm BL}}_{n+1},$

where $m_{n+1}^{\dagger}$ and $C^{\tiny{\rm BL}}_{n+1}$ are given by (202) and
(201) respectively.

######  Proof 8.19.

Note that the transport map (2.5.3) can now be reformulated as

$v_{n+1}=\widehat{v}_{n+1}+(m_{n+1}^{\dagger}-m^{\tiny{\rm BL}}_{n+1}).$

Thus

$v_{n+1}=m_{n+1}^{\dagger}+\widehat{v}_{n+1}-m^{\tiny{\rm BL}}_{n+1}.$

The desired properties of the mean and covariance follow from Lemma 8.16.

######  Remark 8.20.

We now have two derivations of the Kalman gain, the approximate transport
derivation from Subsection 2.5.3, and the minimum variance derivation given
here. We include a third, dimensional, argument that motivates its form. Let
state denote the physical units associated with the state variable and data
those associated with the observation. Then the physical units of the gain
matrix $K$ should equal state/data. We note that $\widehat{C}^{yy}$ has units
of data squared whilst the units of $\widehat{C}^{vy}$ are the product of
state and data. If we then constrain the gain to be determined by covariance
matrices involving the state and the data it is natural to choose it to be
formed as $\widehat{C}^{vy}(\widehat{C}^{yy})^{-1}$ where $\widehat{C}^{yy}$
is an estimate of covariance in the data, and $\widehat{C}^{vy}$ an estimate
of covariance between state and data. Making a choice of this type leads to
the right units for the gain, since the innovation has units data, and relies
only on use of first and second order statistics. This dimensional argument
motivates the form (25) as derived in both Subsections 8.3 and 2.4.
$\blacksquare$

##  9 Appendix D (Stochastic Calculus Considerations)

###  9.1 Derivation of the Kushner-Stratonovich Equation

In Section 3, where we derived continuous time limits from discrete models, we
used the two small parameters $\delta$ and $\Delta t$. The first characterized
the data frequency and the second a time-increment. In this appendix we make
the choice $\delta=\Delta t$ and study the limit $\Delta t\to 0$ to derive
continuum models. Studying this limit enables us to convert between different
modes of stochastic integration in an explicit fashion; this may be helpful to
some readers as it avoids the need to invoke abstract results on covariation.

######  Lemma 9.1.

The Itô and Stratonovich interpretations of the nonlocal nonlinear stochastic
evolution equation for density $r(v,t)$ in (87) are related through

(203a) $\displaystyle dr$
$\displaystyle=r\left\langle\mathsf{h}-\mathbb{E}\mathsf{h},\circ
d{z}^{\dagger}\right\rangle_{\mathsf{\Gamma}}-\frac{r}{2}\left\\{\left|\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}-\mathbb{E}\left|\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}\right\\}dt$
(203b)
$\displaystyle=r\left\langle\mathsf{h}-\mathbb{E}\mathsf{h},d{z}^{\dagger}-\mathbb{E}\mathsf{h}dt\right\rangle_{\mathsf{\Gamma}}$

######  Proof 9.2.

Consider (203a) and the term $r\langle\mathsf{h}-\mathbb{E}\mathsf{h},\circ
d{z}^{\dagger}\rangle_{\mathsf{\Gamma}}$ in particular. Define the
corresponding Itô and Stratonovich integrals

$I:=\int_{0}^{T}r\langle\mathsf{h}-\mathbb{E}\mathsf{h},d{z}^{\dagger}\rangle_{\mathsf{\Gamma}},\quad
S:=\int_{0}^{T}r\langle\mathsf{h}-\mathbb{E}\mathsf{h},\circ
d{z}^{\dagger}\rangle_{\mathsf{\Gamma}}.$

We aim to write $S$ as the sum of $I$ and an added correction. Consider the
space

$\mathcal{X}:=\\{{\varrho}\in
L^{1}(\mathbb{R}^{d_{v}};\mathbb{R}^{+}):\|\varrho\|_{L^{1}}=1\\}$

of probability density functions. We first choose an increasing sequence
$(t_{j})_{j=1,\dots,N}$ with $\Delta t:=t_{j+1}-t_{j}$ such that $N\Delta
t=T$. We next define

$\Delta
r_{t_{j}}:=r_{t_{j+1}}-r_{t_{j}}\quad\Delta\mathsf{h}_{t_{j}}:=\mathsf{h}_{t_{j+1}}-\mathsf{h}_{t_{j}}\quad\Delta{z}^{\dagger}_{t_{j}}:={z}^{\dagger}_{t_{j+1}}-{z}^{\dagger}_{t_{j}}.$

We recall the driving evolution equation for ${z}^{\dagger}$

$d{z}^{\dagger}=\mathsf{h}({v}^{\dagger})dt+\sqrt{\mathsf{\Gamma}}dW,$

from which, recalling the properties of Brownian motion $W$, we deduce the
discretization

(204) $\Delta{z}^{\dagger}_{t_{j}}=\mathsf{h}({v}^{\dagger}_{t_{j}})\Delta
t+\sqrt{\mathsf{\Gamma}\Delta t}\xi_{t_{j}}+\mathcal{O}(\Delta t^{3/2}),$

where $\xi_{t_{j}}$ are independent mean-zero normal random variables with
variance $I_{d_{y}}$ for each $j$. The fact that this discretization is
accurate up to terms of $\mathcal{O}(\Delta t^{3/2})$ follows from Itô-Taylor
expansion Kloeden & Platen (1991). Recalling the definition of the
Stratonovich stochastic integral as a limit, we consider the following finite
sum approximation of $S$:

(205) $S_{\Delta
t}:=\sum_{j=1}^{N}\frac{r_{t_{j+1}}+r_{t_{j}}}{2}\left\langle\mathsf{h}-\int\frac{r_{t_{j+1}}+r_{t_{j}}}{2}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}};$

this converges in the
$L^{2}_{\mathbb{P}}\left(\Omega;\mathcal{C}([0,T];\mathcal{X})\right)$ sense
to $S$ as $\Delta t\rightarrow 0$ (and thus as $N\rightarrow\infty$).
Expanding now the sum in (205), we obtain

(206a) $\displaystyle S_{\Delta t}$
$\displaystyle=\sum_{j=1}^{N}\left(r_{t_{j}}+\frac{\Delta
r_{t_{j}}}{2}\right)\left\langle\mathsf{h}-\int\left(r_{t_{j}}+\frac{\Delta
r_{t_{j}}}{2}\right)\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}$
(206b) $\displaystyle=\sum_{j=1}^{N}r_{t_{j}}\left\langle\mathsf{h}-\int
r_{t_{j}}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}+\frac{1}{2}\sum_{j=1}^{N}\Delta
r_{t_{j}}\left\langle\mathsf{h}-\int
r_{t_{j}}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}$
(206c)
$\displaystyle\quad\quad\quad-\frac{1}{2}\sum_{j=1}^{N}r_{t_{j}}\left\langle\int\Delta
r_{t_{j}}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}-\frac{1}{4}\sum_{j=1}^{N}\Delta
r_{t_{j}}\left\langle\int\Delta
r_{t_{j}}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}.$

Now, we note that by discretizing (203a) we can write

$\Delta r_{t_{j}}=r_{t_{j}}\left\langle\mathsf{h}-\int
r_{t_{j}}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}+\mathcal{O}(\Delta
t).$

By substituting this expression in the expanded sum in (206b) we notice that
the last term in (206c) will be of order $\mathcal{O}(N\Delta t^{3/2})$. We
can thus write (206) as

(207) $S_{\Delta t}=I_{\Delta t}+J_{1,\Delta t}+J_{2,\Delta
t}+\mathcal{O}(N\Delta t^{3/2}),$

where we have defined

$\displaystyle I_{\Delta t}:=$
$\displaystyle\sum_{j=1}^{N}r_{t_{j}}\left\langle\mathsf{h}-\int
r_{t_{j}}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}},\quad\quad
J_{1,\Delta
t}:=\frac{1}{2}\sum_{j=1}^{N}r_{t_{j}}\left|\left\langle\mathsf{h}-\int
r_{t_{j}}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}\right|^{2},$
$\displaystyle J_{2,\Delta
t}:=-\frac{1}{2}\sum_{j=1}^{N}r_{t_{j}}\left\langle\int
r_{t_{j}}\left\langle\mathsf{h}-\int
r_{t_{j}}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}.$

The first term in the sum (207), $I_{\Delta t}$, converges in the
$L^{2}_{\mathbb{P}}$ sense to $I$ as $\Delta t\rightarrow 0$. Considering now
the term $J_{1,\Delta t}$ and using the discretization (204) we have that

$\displaystyle J_{1,\Delta t}$
$\displaystyle=\frac{1}{2}\sum_{j=1}^{N}r_{t_{j}}\left|\left\langle\mathsf{h}-\int
r_{t_{j}}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}\right|^{2}$
$\displaystyle=\frac{1}{2}\sum_{j=1}^{N}r_{t_{j}}\Delta t\left(\mathsf{h}-\int
r_{t_{j}}\mathsf{h}dv\right)^{\top}\mathsf{\Gamma}^{-1/2}\xi_{t_{j}}\xi_{t_{j}}^{\top}\mathsf{\Gamma}^{-1/2}\left(\mathsf{h}-\int
r_{t_{j}}\mathsf{h}dv\right)+\sum_{j=1}^{N}\mathcal{O}(\Delta t^{3/2}).$

Since $\sum_{j=1}^{N}\mathcal{O}(\Delta t^{3/2})=\mathcal{O}(\Delta t^{1/2})$,
taking the $\Delta t\rightarrow 0$ limit and using the independence of the
$\xi_{t_{j}}$ for each $j$, we see that

$J_{1,\Delta
t}\rightarrow\;\frac{1}{2}\int_{0}^{T}\left[r\left(\mathsf{h}-\mathbb{E}\mathsf{h}\right)^{\top}\mathsf{\Gamma}^{-1/2}\mathbb{E}\xi\xi^{\top}\mathsf{\Gamma}^{-1/2}\left(\mathsf{h}-\mathbb{E}\mathsf{h}\right)\right]dt=\int_{0}^{\top}\frac{r}{2}\left|\mathsf{h}-\mathbb{E}\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}dt,$

in the $L^{2}_{\mathbb{P}}$ sense. Similarly for $J_{2,\Delta t}$ we have
that,

$\displaystyle J_{2,\Delta t}$
$\displaystyle=-\frac{1}{2}\sum_{j=1}^{N}r_{t_{j}}\left\langle\int
r_{t_{j}}\left\langle\mathsf{h}-\int
r_{t_{j}}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}$
$\displaystyle=-\frac{1}{2}\sum_{j=1}^{N}r_{t_{j}}\int
r_{t_{j}}\left\langle\left\langle\mathsf{h}-\int
r_{t_{j}}\mathsf{h}dv,\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}\mathsf{h},\Delta{z}^{\dagger}_{t_{j}}\right\rangle_{\mathsf{\Gamma}}dv$
$\displaystyle=-\frac{1}{2}\sum_{j=1}^{N}r_{t_{j}}\Delta t\int
r_{t_{j}}\mathsf{h}^{\top}\mathsf{\Gamma}^{-1/2}\xi_{t_{j}}\xi_{t_{j}}^{\top}\mathsf{\Gamma}^{-1/2}\left(\mathsf{h}-\int
r_{t_{j}}\mathsf{h}dv\right)dv+\sum_{j=1}^{N}\mathcal{O}(\Delta t^{3/2}).$

Since $\sum_{j=1}^{N}\mathcal{O}(\Delta t^{3/2})=\mathcal{O}(\Delta t^{1/2})$,
taking the $\Delta t\rightarrow 0$ limit and using the independence of the
$\xi_{t_{j}}$ for each $j$, we see that, in the $L^{2}_{\mathbb{P}}$ sense,

$J_{2,\Delta t}\rightarrow\;-\frac{1}{2}\int_{0}^{T}r\int
r\mathsf{h}^{\top}\mathsf{\Gamma}^{-1/2}\mathbb{E}\left(\xi\xi^{\top}\right)\mathsf{\Gamma}^{-1/2}\left(\mathsf{h}-\mathbb{E}\mathsf{h}\right)dvdt=-\int_{0}^{T}\frac{r}{2}\mathbb{E}\left|\mathsf{h}-\mathbb{E}\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}dt.$

Finally, by taking the $\Delta t\rightarrow 0$ limit on both sides of (207),
we conclude

$S=I+\int_{0}^{T}\frac{r}{2}\left|\mathsf{h}-\mathbb{E}\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}dt-\int_{0}^{T}\frac{r}{2}\mathbb{E}\left|\mathsf{h}-\mathbb{E}\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}dt,$

in the $L^{2}_{\mathbb{P}}$ sense. Hence by using the above Stratonovich-to-
Îto correction in (203a), we obtain

$\displaystyle dr$
$\displaystyle=r\left\langle\mathsf{h}-\mathbb{E}\mathsf{h},\circ
d{z}^{\dagger}\right\rangle_{\mathsf{\Gamma}}-\frac{r}{2}\left\\{\left|\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}-\mathbb{E}\left|\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}\right\\}dt$
$\displaystyle=r\left\langle\mathsf{h}-\mathbb{E}\mathsf{h},d{z}^{\dagger}\right\rangle_{\mathsf{\Gamma}}+\frac{r}{2}\left|\mathsf{h}-\mathbb{E}\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}dt-\frac{r}{2}\mathbb{E}\left|\mathsf{h}-\mathbb{E}\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}dt-\frac{r}{2}\left\\{\left|\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}-\mathbb{E}\left|\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}\right\\}dt$
$\displaystyle=r\left\langle\mathsf{h}-\mathbb{E}\mathsf{h},d{z}^{\dagger}-\mathbb{E}\mathsf{h}dt\right\rangle_{\mathsf{\Gamma}}.$

######  Lemma 9.3.

The probability density for the solution $v$ of (100), with respect to
randomness induced by the law of $\widehat{z}$ and the initial condition
$v(0)$, with ${z}^{\dagger}$ a fixed data sample path, is given by (101).

######  Proof 9.4.

Recall that in equation (100) the evolution of ${z}^{\dagger}$ is given by
(79), with ${B}^{\dagger}$ and $B$ independent draws from unit Brownian motion
in $\mathbb{R}^{d_{y}}.$ Thus

$\displaystyle dv$
$\displaystyle=a(v;\rho)dt+K(v;\rho)\Bigl{(}\bigl{(}\mathsf{h}({v}^{\dagger})-\mathsf{h}(v)\bigr{)}dt+\sqrt{\mathsf{\Gamma}}d{B}^{\dagger}-\sqrt{\mathsf{\Gamma}}dB\Bigr{)}.$

We wish to find the probability density function $\rho(v,t)$ for $v$, with
respect to randomness induced by the law of $B$ and the initial condition
$v(0)$, but with ${v}^{\dagger},{B}^{\dagger}$ fixed signal and observational
noise sample paths. Let $\phi:\mathbb{R}^{d_{v}}\to\mathbb{R}$ be smooth and
note that the Itô formula shows that

$d\phi\bigl{(}v(t)\bigr{)}=\bigl{(}\mathcal{L}\phi\bigr{)}\bigl{(}v(t)\bigr{)}dt+\big{\langle}K\bigl{(}v(t);\rho(\cdot,t)\bigr{)}\sqrt{\mathsf{\Gamma}}(d{B}^{\dagger}-dB),\nabla\phi\bigl{(}v(t)\bigr{)}\big{\rangle},$

where

$\mathcal{L}\psi(v)=\Big{\langle}\Bigl{(}a(v;\rho)+K(v;\rho)\bigl{(}\mathsf{h}({v}^{\dagger})-\mathsf{h}(v)\bigr{)}\Bigr{)},\nabla\psi\Big{\rangle}+K(v;\rho)\mathsf{\Gamma}K(v;\rho)^{\top}:\nabla\nabla\psi,$

and the integrals are to be interpreted in the Itô sense. Note the factor $1$
in the second order term, arising because of the independent quadratic
variation contributions from both $B$ and ${B}^{\dagger}.$ Taking expectation
$\mathbb{E}$ with respect to $B$ and initial condition $v(0)$, with
${B}^{\dagger}$ fixed, for $a$ and $K$ functions of $v$ and $\rho$, yields,
using again that ${z}^{\dagger}$ is given by (79),

$\displaystyle d\mathbb{E}\phi(v)$
$\displaystyle=\mathbb{E}\Big{\langle}\Bigl{(}a+K\bigl{(}\mathsf{h}({v}^{\dagger})-\mathsf{h}\bigr{)}\Bigr{)},\nabla\phi\Bigr{\rangle}dt+\mathbb{E}\langle
K\sqrt{\mathsf{\Gamma}}d{B}^{\dagger},\nabla\phi\rangle\rangle+\mathbb{E}K\Gamma
K^{\top}:\nabla\nabla\phi dt$
$\displaystyle=\mathbb{E}\big{\langle}\bigl{(}a-K\mathsf{h}\bigr{)},\nabla\phi\bigr{\rangle}dt+\mathbb{E}\langle
Kd{z}^{\dagger},\nabla\phi\rangle\rangle+\mathbb{E}K\mathsf{\Gamma}K^{\top}:\nabla\nabla\phi
dt.$

Noting that the expectation operation corresponds to multiplication by
$\rho(v,t)$ and integration over $v\in\mathbb{R}^{d_{v}},$ integrating by
parts shows that $\rho$ satisfies the desired equation (101), in a weak sense.

###  9.2 Diamond Integration

In this section we use notation analogous to that established in the proof of
Lemma 9.1. Consider the equation (99), and in particular the contribution
$K(d{z}^{\dagger}-d\widehat{z})$. We first choose an increasing sequence
$(t_{j})_{j=0,\dots,N-1}$ with $\Delta t:=t_{j+1}-t_{j}$ such that $N\Delta
t=T$. The Itô integral interpretation of this contribution, on time interval
$(0,T)$, is as the $L^{2}_{\mathbb{P}}$ limit of

$I_{\Delta
t}:=\sum_{j=0}^{N-1}K(v_{t_{j}};\rho_{t_{j}})(\Delta{z}^{\dagger}_{t_{j}}-\Delta\widehat{z}_{t_{j}})$

We define the Stratonovich integral interpretation as the $L^{2}_{\mathbb{P}}$
limit of

$S_{\Delta
t}:=\sum_{j=0}^{N-1}K\left(\frac{v_{t_{j+1}}+v_{t_{j}}}{2};\frac{\rho_{t_{j+1}}+\rho_{t_{j}}}{2}\right)(\Delta{z}^{\dagger}_{t_{j}}-\Delta\widehat{z}_{t_{j}});$

as is standard, we use $\circ$ to denote Stratonovich stochastic integration.
Finally we define the diamond integral interpretation as the
$L^{2}_{\mathbb{P}}$ limit of

$R_{\Delta
t}:=\sum_{j=0}^{N-1}K\left(\frac{v_{t_{j+1}}+v_{t_{j}}}{2};\rho_{t_{j}}\right)(\Delta{z}^{\dagger}_{t_{j}}-\Delta\widehat{z}_{t_{j}}).$

Note that this is akin to a Stratonovich integral, but only with respect
variation of $K$ with respect to $v$, not $\rho.$ We use $\diamond$ to denote
this unusual form of stochastic integration.

Throughout this subsection we inter-convert between these three forms of
stochastic integration. To shorten the presentation we will sometimes use the
$=$ symbol when in fact we mean equality up to an additive constant which
disappears in the $\Delta t\to 0$ limit. Note that the evolution equation (99)
for $v$ is driven by ${z}^{\dagger}$ and $\widehat{z}$, whereas the Kushner-
Stratonovich equation (89) for $\rho$ driven only by ${z}^{\dagger}.$ This
difference will have implications for the calculations that follow, in which
we compute inter-conversions between the different stochastic integrals. The
first result highlights an interesting interpretation of the contribution

$a=\nabla\cdot(K\mathsf{\Gamma}K^{\top})-K\mathsf{\Gamma}\nabla\cdot
K^{\top};$

to the drift in (99), namely that it is simply the Itô-to-diamond correction,
giving a compact re-interpretation of the mean-field model:

######  Lemma 9.5.

Let $\rho$ solve the Kushner-Stratonovich equation (89) and let
${z}^{\dagger}$ be given by (79). The Itô interpretation of equation (99) and
its interpretation with respect to the $\diamond$ form of stochastic
integration are related through the following equivalence. The system

$\displaystyle dv$
$\displaystyle=f(v)dt+\sqrt{\mathsf{\Sigma}}dW+\nabla\cdot(K\mathsf{\Gamma}K^{\top})dt-K\mathsf{\Gamma}\nabla\cdot
K^{\top}dt+K(d{z}^{\dagger}-d\widehat{z}),$ $\displaystyle d\widehat{z}$
$\displaystyle=\mathsf{h}(v)dt+\sqrt{\mathsf{\Gamma}}dB,$

where $K=K(v,\rho)$, is equivalent to

$\displaystyle dv$
$\displaystyle=f(v)dt+\sqrt{\mathsf{\Sigma}}dW+K\diamond(d{z}^{\dagger}-d\widehat{z})$
$\displaystyle d\widehat{z}$
$\displaystyle=\mathsf{h}(v)dt+\sqrt{\mathsf{\Gamma}}dB.$

######  Proof 9.6.

We first recall the evolution equations for ${z}^{\dagger}$ and $\widehat{z}$:

$\displaystyle d{z}^{\dagger}$
$\displaystyle=\mathsf{h}({v}^{\dagger})dt+\sqrt{\mathsf{\Gamma}}d{B}^{\dagger},$
$\displaystyle d\widehat{z}$
$\displaystyle=\mathsf{h}(v)dt+\sqrt{\mathsf{\Gamma}}dB.$

Noting that $v$ is driven by $({z}^{\dagger}-\widehat{z})$ and that we are
interconverting between Itô and diamond integration, so that the noisy driving
of the $\rho$ equation does not play a role, we see that it suffices to
consider the $L^{2}_{\mathbb{P}}$ limits of the two quantities

$\displaystyle I_{\Delta t}$
$\displaystyle=\sum_{j=0}^{N-1}K(v_{t_{j}};\rho_{t_{j}})\sqrt{2\mathsf{\Gamma}\Delta
t}\xi_{t_{j}},$ $\displaystyle R_{\Delta t}$
$\displaystyle=\sum_{j=0}^{N-1}K\left(\frac{v_{t_{j+1}}+v_{t_{j}}}{2};\rho_{t_{j}}\right)\sqrt{2\mathsf{\Gamma}\Delta
t}\xi_{t_{j}},$

where the $\xi_{t_{j}}$ are i.i.d. draws from a unit Gaussian. By adding and
subtracting the Itô contribution, we obtain

$\displaystyle R_{\Delta t}$
$\displaystyle=\sum_{j=0}^{N-1}K\left(\frac{v_{t_{j+1}}+v_{t_{j}}}{2};\rho_{t_{j}}\right)\sqrt{2\mathsf{\Gamma}\Delta
t}\xi_{t_{j}}$
$\displaystyle=\sum_{j=0}^{N-1}K(v_{t_{j}};\rho_{t_{j}})\sqrt{2\mathsf{\Gamma}\Delta
t}\xi_{t_{j}}+\sum_{j=0}^{N-1}\left(K\left(\frac{v_{t_{j+1}}+v_{t_{j}}}{2};\rho_{t_{j}}\right)-K(v_{t_{j}};\rho_{t_{j}})\right)\sqrt{2\mathsf{\Gamma}\Delta
t}\xi_{t_{j}}$ $\displaystyle=I_{\Delta
t}+\sum_{j=0}^{N-1}\left(\frac{1}{2}D_{v}K(v_{t_{j}};\rho_{t_{j}})(v_{t_{j+1}}-v_{t_{j}})+\mathcal{O}\left(|v_{t_{j+1}}-v_{t_{j}}|^{2}\right)\right)\sqrt{2\mathsf{\Gamma}\Delta
t}\xi_{t_{j}},$

where the last line follows from a first order Taylor expansion. Using a
discretization of the evolution for $v$, and neglecting the terms that do not
contribute to the quadratic variation when computing the $L^{2}_{\mathbb{P}}$
limit, we substitute for $v_{t_{j+1}}-v_{t_{j}}$ to obtain

(208) $I_{\Delta
t}+\sum_{j=0}^{N-1}\frac{1}{2}D_{v}K(v_{t_{j}};\rho_{t_{j}})\left(K(v_{t_{j}};\rho_{t_{j}})\sqrt{2\mathsf{\Gamma}\Delta
t}\xi_{t_{j}}\right)\sqrt{2\mathsf{\Gamma}\Delta t}\xi_{t_{j}}.$

We now consider the correction term resulting from (208), which is determined
by expectation of the summand with respect to the random increments $\xi$,
scaled by $\Delta t^{-1}.$ Dropping the $v$ and $\rho$ dependence of $K$ and
the $t_{j}$ dependence for notational convenience, its $k$’th component is
given by222222Using Einstein summation convention, as described for example in
Gonzalez & Stuart (2008), so that index repeated twice is summed over.

$\displaystyle\left[\mathbb{E}\left(D_{v}K\right)\left(K\mathsf{\Gamma}^{1/2}\xi\right)\mathsf{\Gamma}^{1/2}\xi\right]_{k}$
$\displaystyle=\left[\mathbb{E}\left(\partial_{v_{i}}K\right)\left[K\mathsf{\Gamma}^{1/2}\xi\right]_{i}\mathsf{\Gamma}^{1/2}\xi\right]_{k}$
$\displaystyle=\left[\mathbb{E}\left(\partial_{v_{i}}K\right)\left(K_{il}(\mathsf{\Gamma}^{1/2})_{lj}\xi_{j}\right)\mathsf{\Gamma}^{1/2}\xi\right]_{k}$
$\displaystyle=\mathbb{E}\left(K_{il}(\mathsf{\Gamma}^{1/2})_{lj}\xi_{j}\right)\left(\partial_{v_{i}}K\right)_{kn}(\mathsf{\Gamma}^{1/2})_{nm}\xi_{m}\;\;\;\;$
$\displaystyle=\left(\partial_{v_{i}}K_{kn}\right)\mathsf{\Gamma}_{ln}K_{il}.$

But,

$\displaystyle\left(\partial_{v_{i}}K_{kn}\right)\mathsf{\Gamma}_{ln}K_{il}$
$\displaystyle=\partial_{v_{i}}\left(K_{kn}\mathsf{\Gamma}_{ln}K_{il}\right)-K_{kn}\mathsf{\Gamma}_{ln}\left(\partial_{v_{i}}K_{il}\right)$
$\displaystyle=\left[\nabla\cdot\left(K\mathsf{\Gamma}K^{\top}\right)\right]_{k}-\left[K\mathsf{\Gamma}\nabla\cdot\left(K^{\top}\right)\right]_{k},$

Recalling that the only contributions in $R_{\Delta t}$ which do not vanish
under the $\Delta t\rightarrow 0$ limit are the ones given by (208), taking
the $\Delta t\rightarrow 0$ limit of $R_{\Delta t}$ yields

$R=I+\int_{0}^{T}\left(\nabla\cdot\left(K\mathsf{\Gamma}K^{\top}\right)-K\mathsf{\Gamma}\nabla\cdot\left(K^{\top}\right)\right)dt,$

where the convergence is in the
$L^{2}_{\mathbb{P}}\left(\Omega;\mathcal{C}([0,T];\mathbb{R}^{d_{v}})\right)$
sense; this concludes the proof.

######  Lemma 9.7.

Let $\rho$ solve the Kushner-Stratonovich equation (89) and let
${z}^{\dagger}$ be given by (79). The interpretation with respect to the
$\diamond$ form of stochastic integration of equation (99) and its
Stratonovich interpretation are related through the following equivalence. The
system

$\displaystyle dv$
$\displaystyle=f(v)dt+\sqrt{\mathsf{\Sigma}}dW+K\diamond(d{z}^{\dagger}-d\widehat{z})$
$\displaystyle d\widehat{z}$
$\displaystyle=\mathsf{h}(v)dt+\sqrt{\mathsf{\Gamma}}dB,$

where $K=K(v,\rho)$, is equivalent to

$\displaystyle dv$
$\displaystyle=f(v)dt+bdt+\sqrt{\mathsf{\Sigma}}dW+K\circ(d{z}^{\dagger}-d\widehat{z})$
$\displaystyle d\widehat{z}$
$\displaystyle=\mathsf{h}(v)dt+\sqrt{\mathsf{\Gamma}}dB,$

where the term $b=b(v,\rho)$ satisfies

$b(v,\rho)=\frac{\rho}{2}D_{\rho}K(\mathsf{h}-\mathbb{E}\mathsf{h}).$

######  Proof 9.8.

To make the interconversion between diamond and Stratonovich integration we
need only consider the quadratic variation contribution induced by
${z}^{\dagger}$ since the equation for $\rho$ is driven only by
${z}^{\dagger}$ and not by $\widehat{z}.$ Thus it suffices to compare the
$L^{2}_{\mathbb{P}}$ limits of

$\displaystyle R_{\Delta t}$
$\displaystyle=\sum_{j=0}^{N-1}K\left(\frac{v_{t_{j+1}}+v_{t_{j}}}{2};\rho_{t_{j}}\right)\sqrt{\mathsf{\Gamma}\Delta
t}\xi^{\dagger}_{t_{j}},$ $\displaystyle S_{\Delta t}$
$\displaystyle=\sum_{j=0}^{N-1}K\left(\frac{v_{t_{j+1}}+v_{t_{j}}}{2};\frac{\rho_{t_{j+1}}+\rho_{t_{j}}}{2}\right)\sqrt{\mathsf{\Gamma}\Delta
t}\xi^{\dagger}_{t_{j}};$

Here $\xi^{\dagger}_{t_{j}}$ are i.i.d. draws from a unit Gaussian. By adding
and subtracting the diamond contribution, we obtain

$\displaystyle S_{\Delta t}$
$\displaystyle=\sum_{j=0}^{N-1}K\left(\frac{v_{t_{j+1}}+v_{t_{j}}}{2};\frac{\rho_{t_{j+1}}+\rho_{t_{j}}}{2}\right)\sqrt{\mathsf{\Gamma}\Delta
t}\xi^{\dagger}_{t_{j}}$ $\displaystyle=R_{\Delta
t}+\sum_{j=0}^{N-1}\left(K\left(\frac{v_{t_{j+1}}+v_{t_{j}}}{2};\frac{\rho_{t_{j+1}}+\rho_{t_{j}}}{2}\right)-K\left(\frac{v_{t_{j+1}}+v_{t_{j}}}{2};\rho_{t_{j}}\right)\right)\sqrt{\mathsf{\Gamma}\Delta
t}\xi^{\dagger}_{t_{j}}$ $\displaystyle=R_{\Delta
t}+\sum_{j=0}^{N-1}\left(\frac{1}{2}D_{\rho}K\left(\frac{v_{t_{j+1}}+v_{t_{j}}}{2};\rho_{t_{j}}\right)(\rho_{t_{j+1}}-\rho_{t_{j}})+\mathcal{O}\left(|\rho_{t_{j+1}}-\rho_{t_{j}}|^{2}\right)\right)\sqrt{\mathsf{\Gamma}\Delta
t}\xi^{\dagger}_{t_{j}},$

where the last line follows from a first order Taylor expansion in $\rho$.
Using a discretization of the Kushner-Stratonovich equation given in Theorem
3.5, and neglecting the terms that do not contribute to the quadratic
variation in the $L^{2}_{\mathbb{P}}$ limit, we substitute for
$\rho_{t_{j+1}}-\rho_{t_{j}}$ to obtain

$S_{\Delta t}=R_{\Delta
t}+\sum_{j=0}^{N-1}\left(\frac{1}{2}D_{\rho}K\left(\frac{v_{t_{j+1}}+v_{t_{j}}}{2};\rho_{t_{j}}\right)\langle\rho_{t_{j}}\mathsf{\Gamma}^{-1}(\mathsf{h}_{t_{j}}-\mathbb{E}\mathsf{h}_{t_{j}}),\sqrt{\mathsf{\Gamma}\Delta
t}\xi^{\dagger}_{t_{j}}\rangle+\mathcal{O}(\Delta