us a sample path $\\{{z}^{\dagger}(t)\\}_{t\in\mathbb{R}^{+}}$ of the
observation coordinates of a realization of the SDE (75), from which we wish
to recover the true realization of the state
$\\{{v}^{\dagger}(t)\\}_{t\in\mathbb{R}^{+}}.$ These observation and state
sample paths are generated by
${v}^{\dagger}_{0},\\{{W}^{\dagger}\\}_{t\in\mathbb{R}^{+}}$ and
$\\{{B}^{\dagger}\\}_{t\in\mathbb{R}^{+}}$, specific realizations of the
initial condition and the Brownian motions driving the state and observation
components of (75). We also introduce
$Z^{\dagger}(t)=\\{{z}^{\dagger}(s)\\}_{0\leq s\leq t}.$ It is natural to
think about our objective in two different ways:

  * •

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

  * •

Objective 2: design an algorithm which estimates the distribution of random
variable $v(t)|Z^{\dagger}(t)$.

In both cases we are interested in Markovian formulations which update the
estimate $v(t)$, or the distribution $v(t)|Z^{\dagger}(t)$, sequentially as
the data is acquired. All of the algorithms we describe depend only on the
increments of the process ${z}^{\dagger}(t)$, hence the fixing of
${z}^{\dagger}(0)=0$ is immaterial. In the next two subsections we describe
control theoretic and probabilistic approaches to this problem which,
respectively, provide the basis for algorithms addressing Objectives 1 and 2.
Following the road-map from the previous section in the discrete time setting,
we then proceed to study exact transport leading to mean-field equations
related to these Objective 2; we then study second order approximations of
exact transport, and finally reach ensemble Kalman methods through particle
approximations.

###  3.2 Control Theory Perspective

As for the time-discrete problem, we again start with the control theoretical
approach based on the small uncertainty assumption. Specifically we assume
that the three covariances appearing in (76) are small so that the states and
observations can be well approximated as deterministic. Thus we initially set
$\mathsf{\Sigma}$ and $\mathsf{\Gamma}$ to zero. Applying the rescalings (72)
to (3) we obtain, in the deterministic setting,

(77a) $\displaystyle\widehat{v}_{n+1}$ $\displaystyle=v_{n}+\Delta tf(v_{n}),$
(77b) $\displaystyle\widehat{z}_{n+1}$ $\displaystyle=\widehat{z}_{n}+\Delta
t\mathsf{h}(\widehat{v}_{n+1}),$ (77c) $\displaystyle v_{n+1}$
$\displaystyle=\widehat{v}_{n+1}+K(\Delta{z}^{\dagger}_{n+1}-\Delta\widehat{z}_{n+1}).$

Equations (77) define a deterministic map $v_{n}\mapsto v_{n+1}$ using the
observed increments from (73b) so that

$\displaystyle\Delta{z}^{\dagger}_{n+1}$
$\displaystyle:={z}^{\dagger}\big{(}(n+1)\Delta
t\big{)}-{z}^{\dagger}\big{(}n\Delta t\big{)},$

derived from a specific fixed realization of (75). Taking the continuous time
limit, and eliminating $\widehat{z}$, we obtain the following estimator $v$
for ${v}^{\dagger}$ given $Z^{\dagger}:$

(78) ${dv}=f(v)dt+K\bigl{(}{d{z}^{\dagger}}-\mathsf{h}(v)dt\bigr{)},$

where ${z}^{\dagger}$ (the data) is obtained from a specific fixed realization
of (75):

(79a) $\displaystyle d{v}^{\dagger}$
$\displaystyle=f({v}^{\dagger})dt+\sqrt{\mathsf{\Sigma}}dW^{\dagger},\quad
v_{0}\sim\mathsf{N}(m_{0},C_{0})$ (79b) $\displaystyle d{z}^{\dagger}$
$\displaystyle=\mathsf{h}({v}^{\dagger})dt+\sqrt{\mathsf{\Gamma}}dB^{\dagger},\quad
z(0)=0.$

This is a continuous time analog of the 3DVAR algorithm (4). Equation (78) is
a controlled ordinary differential equation (ODE); typically it is initialized
with $v(0)\sim\mathsf{N}(m_{0},C_{0})$.

We now include the effect of uncertainty, allowing for non-zero covariances in
(76). Accounting for noise in the expressions for $\widehat{v}_{n+1}$ and
$\widehat{z}_{n+1}$ in (77) we obtain the following rescaling of (14):

$\displaystyle\widehat{v}_{n+1}$ $\displaystyle=v_{n}+\Delta
tf(v_{n})+\xi_{n},$ $\displaystyle\widehat{z}_{n+1}$
$\displaystyle=\widehat{z}_{n}+\Delta
t\mathsf{h}(\widehat{v}_{n+1})+\eta_{n+1},$ $\displaystyle v_{n+1}$
$\displaystyle=\widehat{v}_{n+1}+K_{n}(\Delta{z}^{\dagger}_{n+1}-\Delta\widehat{z}_{n+1}),$

for all $n\in\mathbb{Z}^{+}$, where we assume
$v_{0},\\{\xi_{n}\\}_{n\in\mathbb{Z}^{+}}$ and
$\\{\eta_{n}\\}_{n\in\mathbb{N}}$ are mutually independent Gaussians defined
by (3.1). We may now formally take the $\Delta t\to 0$ limit and obtain the
following continuous time analog of (14), namely the controlled SDE
formulation:

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

where ${z}^{\dagger}$ is given by (79). The unit Brownian motions
$W,W^{\dagger},B$ and $B^{\dagger}$, in $\mathbb{R}^{d_{v}}$,
$\mathbb{R}^{d_{v}}$, $\mathbb{R}^{d_{z}}$ and $\mathbb{R}^{d_{z}}$
respectively, are all independent of one another. As in the discrete time
analog, encapsulated in (14), the choice of a time-dependent gain matrix $K$
remains to be determined. We adopt a mean-field perspective for determination
of $K$ in continuous time, as we did in discrete time. To this end we progress
to discuss the evolution of probability measures describing the conditional
distribution of $v(t)|Z^{\dagger}(t)$.

###  3.3 Probabilistic Perspective

####  3.3.1 The Filtering Distribution

The derivation of continuous time limits in the previous two subsections is
relatively straightforward. However, there is an important practical and
theoretical issue which we need to address. Instead of continuous time
observations ${z}^{\dagger}(t)$ one typically has access to the data at
discrete time instances $\tau_{k}=k\delta$, $\delta>0$, $k=1,\ldots$ only. In
order to make use of continuous time algorithms and theory it is then useful
to construct a continuous time approximation ${z}^{{\dagger},\delta}$; to be
concrete we will use piecewise linear interpolation. With this set-up we need
to deal with two small parameters; the time-step $\Delta t$ used in (72) to
obtain a continuous time limit, and the data sampling interval $\delta$. There
are results for continuous time filtering which imply that the desired
limiting equation can be found in either Itô or Stratonovich forms by
considering different orders of the limits $\Delta t\to 0$ and $\delta\to 0$ –
see the bibliographic remarks in Subsection 3.7. However most theoretical
results are derived by first taking $\delta\to 0$ and then $\Delta t\to 0$.
From a practical and theoretical perspective it is sometimes more convenient
to first consider $\Delta t\to 0$ followed by $\delta\to 0$. We will utilize
the later sequence of limits in the following subsection in order to derive a
set of evolution equations for the conditional probability measure $\mu(v,t)$
solving Objective 2. These equations will in turn guide our choice of the gain
matrix $K$ in (80a). Finally, when the data ${z}^{\dagger}$ arises itself from
numerical simulations, then it is most convenient to set $\delta=\Delta t$,
and hence $\tau_{k}=t_{k}$; this is implicitly used in the derivation of the
continuous time sample path equations in the two preceding subsections.

Our starting point here is the iteration on measures, the filtering cycle,
defined by (19), under the scalings (72). We assume that $\delta$ is an
integer multiple of $\Delta t$ so that the $\\{\tau_{k}\\}$ are a subset of
the $\\{t_{n}\\}.$ We replace the true observation path ${z}^{\dagger}(t)$ by
its piecewise linear approximation ${z}^{{\dagger},\delta}(t)$ based on linear
interpolation of values $\\{{z}^{\dagger}(\tau_{k})\\}.$ We then have that the
implied observation increments $\Delta{z}^{\dagger}_{n+1}$ are constant over
the time intervals $[t_{n},t_{n+1}]$ and are given by

$\displaystyle\Delta{z}^{\dagger}_{n+1}=\frac{d{z}^{{\dagger},\delta}}{dt}(t_{n})\Delta
t.$

In this setting, the operators $P$ and $L_{n}$, defined by (19a) and (19b),
respectively, become

(81a) $\displaystyle(P\mu)(dv)$
$\displaystyle=\Bigl{(}\int_{u\in\mathbb{R}^{d_{v}}}p(u,v)\mu(du)\Bigr{)}dv,$
(81b) $\displaystyle L_{n}(\mu)(dv)$
$\displaystyle=q(v,\Delta{z}^{\dagger}_{n+1})\mu(dv)\Big{/}\Bigl{(}\int_{\mathbb{R}^{d_{v}}}q(v,\Delta{z}^{\dagger}_{n+1})\mu(dv)\Bigr{)},$

where

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

With these scalings, we now derive the continuous time analog of (19), for
$\mu(v,t).$ For ease of exposition we assume that $\mu$ has density $\rho$ and
derive the equation satisfied by $\rho.$ To do this we employ the split-step
principle: we find the continuous time evolution equation associated with each
of $P$ and $L_{n}$ (equations (81a) and (81b) respectively) separately, and
then add the right-hand sides of the resulting evolution equations to obtain
the desired continuous time limit resulting from the composition of $L_{n}$
and $P$. We use $r$ as a dummy variable to denote the density being evolved,
for both of the split-steps, and in both discrete ($r_{n}$) and continuous
($r(t)$) time, in what follows.

First consider the continuous time limit of the evolution associated to $P$ as
described by (81a) and (82a). Since the underlying continuous time limit of
the sample paths is given by the SDE (75a), the time evolution of the
probability density $r$ is given by the Fokker–Planck equation 111111Here, and
in what follows, we use the standard notation from continuum mechanics for the
divergence of vector and second order tensor fields, and for the gradient of
scalar and vector fields; see the bibliography Subsection 3.7 for references.

(83)
$\displaystyle\partial_{t}r=-\nabla\cdot(rf)+\frac{1}{2}\nabla\cdot(\nabla\cdot(r\mathsf{\Sigma})).$

Consider now the second component of the split-step argument: we determine a
continuous time limit of the evolution associated to $L_{n}$ as described by
(81b) and (82b). The following lemma presents an evolution equation for $r$
associated to $L_{n}$, describing how observation of the piecewise continuous
interpolated data ${z}^{{\dagger},\delta}(t)$ changes the density $r(t,v)$.

######  Lemma 3.1.

Assume that $\mathsf{\Gamma}\succ 0.$ The continuous time limit of the
evolution associated to $L_{n}$, as described by (81b) and (82b), is given by

(84)
$\displaystyle\partial_{t}r=\left\langle\mathsf{h}-\mathbb{E}\mathsf{h},\frac{d{z}^{{\dagger},\delta}}{dt}\right\rangle_{\mathsf{\Gamma}}r-\frac{1}{2}\left\\{\left|\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}-\mathbb{E}\left|\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}\right\\}r.$

######  Proof 3.2.

Consider the discrete-time evolution $r_{n+1}=L_{n}r_{n}$, where $L_{n}$ is
defined by (81b) and (82b). By Taylor expansion we have

$\exp\left(-\frac{1}{2\Delta t}\left|\Delta{z}^{\dagger}_{n+1}-\Delta
t\mathsf{h}(v)\right|_{\mathsf{\Gamma}}^{2}\right)=1-\frac{\Delta
t}{2}\left|\frac{\Delta{z}^{\dagger}_{n+1}}{\Delta
t}-\mathsf{h}(v)\right|^{2}_{\mathsf{\Gamma}}+\mathcal{O}(\Delta t^{2});$

then, using expressions (81b) and (82b), we obtain

(85) $L_{n}r_{n}=\frac{1}{C}\left(1-\frac{\Delta
t}{2}\left|\frac{\Delta{z}^{\dagger}_{n+1}}{\Delta
t}-\mathsf{h}(v)\right|^{2}_{\mathsf{\Gamma}}+\mathcal{O}(\Delta
t^{2})\right)r_{n},$

where

$C=\int_{\mathbb{R}^{d_{v}}}\left(1-\frac{\Delta
t}{2}\left|\frac{\Delta{z}^{\dagger}_{n+1}}{\Delta
t}-\mathsf{h}(v)\right|^{2}_{\mathsf{\Gamma}}+\mathcal{O}(\Delta
t^{2})\right)r_{n}dv.$

By integrating and noting that $r_{n}$ is a density, it follows that

(86) $C=1-\frac{\Delta
t}{2}\mathbb{E}\left|\frac{\Delta{z}^{\dagger}_{n+1}}{\Delta
t}-\mathsf{h}(v)\right|^{2}_{\mathsf{\Gamma}}+\mathcal{O}(\Delta t^{2}),$

where expectation is with respect to $v$ distributed as random variable with
density $r_{n}$. Hence, multiplying numerator and denominator in (85) by
$1+\frac{\Delta t}{2}\mathbb{E}\left|\frac{\Delta{z}^{\dagger}_{n+1}}{\Delta
t}-\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}$ we obtain

$r_{n+1}=r_{n}\left\\{1+\Delta
t\left\langle\mathsf{h}-\mathbb{E}\mathsf{h},\frac{d{z}^{{\dagger},\delta}}{dt}\right\rangle_{\mathsf{\Gamma}}-\frac{\Delta
t}{2}\left\\{\left|\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}-\mathbb{E}\left|\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}\right\\}+\mathcal{O}(\Delta
t^{2})\right\\},$

using that the data increments $\Delta{z}^{\dagger}_{n+1}$ are given by
(3.3.1) for fixed $\delta$. Taking the time continuous limit $\Delta t\to 0$
with fixed observation interval $\delta>0$ leads to the evolution equation
(84).

Now taking the $\delta\to 0$ limit in (84), we obtain the following nonlocal
nonlinear stochastic evolution equation for density $r(v,t):$

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

Here $\circ$ denotes Stratonovitch integration; this form of integration
arises in the limit $\delta\to 0$ because the equation is derived by making a
smooth approximation ${z}^{{\dagger},\delta}$ of ${z}^{\dagger}$ and passing
to the limit. Recall that ${z}^{\dagger}$ is given by (79). The equation is
nonlocal and nonlinear because $\mathbb{E}$ denotes expectation at time $t$
with respect to density $r(\cdot,t).$

We now wish to invoke the split-step principle and combine the evolutions (83)
and (87). However before doing this we proceed to convert the equation (87)
into its more common Itô representation. For this, the following lemma is
crucial. Here we provide a streamlined proof using known abstract results
about covariation of correlated stochastic processes; in Appendix 9 we provide
explicit calculations for the reader who is interested in understanding the
details of the conversion by means of the definitions of Itô and Stratonovitch
integrals as limits. In particular, the concepts underlying the covariation
calculations in the following are derived from first principles in Lemma 9.1.

######  Lemma 3.3.

Assume that $\mathsf{\Gamma}\succ 0$ and that ${z}^{\dagger}$ is given by
(79). The Itô and Stratonovich interpretations of the stochastic forcing term
in (87) are related through

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

######  Proof 3.4.

Using the formula for the Itô–Stratonovich conversion between two
semimartingales, we have

(88)
$\displaystyle\int_{0}^{t}r\left\langle\mathsf{h}-\mathbb{E}\mathsf{h},\circ
d{z}^{\dagger}\right\rangle_{\mathsf{\Gamma}}=\int_{0}^{t}r\left\langle\mathsf{h}-\mathbb{E}\mathsf{h},d{z}^{\dagger}\right\rangle_{\mathsf{\Gamma}}+\frac{1}{2}\int_{0}^{t}d\left[\left\langle
r(\mathsf{h}-\mathbb{E}\mathsf{h}),{z}^{\dagger}\right\rangle_{\mathsf{\Gamma}}\right],$

where $[\langle\cdot,\cdot\rangle]$ denotes covariation. Note that the
covariation between $r$ and ${z}^{\dagger}$ satisfies
$d[r,{z}^{\dagger}]=r(\mathsf{h}-\mathbb{E}\mathsf{h})dt$. Hence the
covariation in (88) satisfies

$\displaystyle d[\langle
r(\mathsf{h}-\mathbb{E}\mathsf{h}),{z}^{\dagger}\rangle_{\mathsf{\Gamma}}]$
$\displaystyle=\langle\mathsf{h}-\mathbb{E}\mathsf{h},d[r,{z}^{\dagger}]\rangle_{\mathsf{\Gamma}}-rd[\langle\mathbb{E}\mathsf{h},{z}^{\dagger}\rangle_{\mathsf{\Gamma}}]$
$\displaystyle=r\left|\mathsf{h}-\mathbb{E}\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}dt-r\int\langle\mathsf{h},d[r,{z}^{\dagger}]\rangle_{\mathsf{\Gamma}}dv$
$\displaystyle=r\left|\mathsf{h}-\mathbb{E}\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}dt-r\mathbb{E}\langle\mathsf{h},(\mathsf{h}-\mathbb{E}\mathsf{h})\rangle_{\mathsf{\Gamma}}dt$
$\displaystyle=r\left\\{\left|\mathsf{h}-\mathbb{E}\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}-\mathbb{E}\left|\mathsf{h}-\mathbb{E}\mathsf{h}\right|^{2}_{\mathsf{\Gamma}}\right\\}dt.$

Upon rearranging, we find that

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

which in turn leads to the following Itô representation of (87):

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

Using the Itô form of the equation from the preceding lemma, the Fokker-Planck
equation (83) and the split-step principle delivers:

######  Theorem 3.5.

Assume that $\mathsf{\Gamma}\succ 0$ and that ${z}^{\dagger}$ is given by
(79). The time evolution of the density $\rho(\cdot,t)$ for the random
variable $v(t)|Z^{\dagger}(t)$ is characterized by the SPDE

(89) $d\rho=-\nabla\cdot(\rho
f)dt+\frac{1}{2}\nabla\cdot(\nabla\cdot(\rho\mathsf{\Sigma}))dt+\left\langle\mathsf{h}-\mathbb{E}\mathsf{h},d{z}^{\dagger}-\mathbb{E}\mathsf{h}dt\right\rangle_{\mathsf{\Gamma}}\rho.$

This is known as the Kushner–Stratonovich equation. The equation (89) should
now be interpreted in the Itô sense with respect to the driving noise
${z}^{\dagger}(t)$. The corresponding Stratonovich formulation follows
immediately from (87) in combination with (83).

####  3.3.2 Important Repeatedly Used Notation

In discrete time expectations may be taken under the law of the (possibly
approximate) discretely evolving mean-field model, or under the predictive
distribution found from pushing this law forward under the model. This
distinction disappears in continuous time and we simply need to compute
expectations under the (possibly approximate) continuously evolving mean-field
models that we will introduce later. To this end we define, with
$m:=\mathbb{E}v$,

(90a) $\displaystyle C$
$\displaystyle=\mathbb{E}\bigl{(}(v-m)\otimes(v-m)\bigr{)},$ (90b)
$\displaystyle C^{vf}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}v-m\bigr{)}\otimes\bigl{(}f(v)-\mathbb{E}f(v)\bigr{)}\Bigr{)},$
(90c) $\displaystyle C^{v\mathsf{h}}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}v-m\bigr{)}\otimes\bigl{(}\mathsf{h}(v)-\mathbb{E}\mathsf{h}(v)\bigr{)}\Bigr{)}.$

All expectations are under the mean-field model for $v$. The covariances
should be viewed as functions of time $t$. In deriving continuous time models
from discrete time models, we will also use discrete time analogs, computed
under the law of random variable $v_{n}$ and denoted $C_{n},C_{n}^{vf}$ and
$C_{n}^{v\mathsf{h}}.$

###  3.4 Gaussian Projected Filtering Distribution

As in discrete time, the Gaussian projected filtering distribution plays an
important conceptual role in understanding later filtering algorithms. We
derive this Gaussian filter in Subsection 3.4.1, and specialize to the setting
of linear-Gaussian dynamics-observation in Subsection 3.4.2 where we obtain
the Kalman-Bucy filter.

####  3.4.1 Gaussian Projected Filtering

The evolution equations for the mean $m$ and the covariance matrix $C$ follow
naturally from a continuous time limit of the associated discrete time
filtering formulations. We summarize the resulting equations in the following:

######  Theorem 3.6.

Assume that $\mathsf{\Gamma}\succ 0.$ and that ${z}^{\dagger}$ is given by
(79). Consider the map from $(m_{n},C_{n})$ to $(m_{n+1},C_{n+1})$ defined by
choosing $v_{n}\sim\mathsf{N}(m_{n},C_{n})$ and then using (21), (23), (26)
and (31). Under the rescalings (72), and in the limit $\Delta t\to 0$, we
obtain the following continuous-time analogue of this map:

(91a) $\displaystyle dm$
$\displaystyle=\mathbb{E}f(v)dt+C^{v\mathsf{h}}\mathsf{\Gamma}^{-1}\bigl{(}d{z}^{\dagger}-\mathbb{E}\mathsf{h}(v)dt\bigr{)}$
(91b) $\displaystyle dC$
$\displaystyle=C^{vf}dt+(C^{vf})^{\top}dt+\mathsf{\Sigma}dt-C^{v\mathsf{h}}\mathsf{\Gamma}^{-1}(C^{v\mathsf{h}})^{\top}dt,$

where the expectation is computed under $\mathsf{N}\bigl{(}m(t),C(t)\bigr{)}$
and where ${z}^{\dagger}$ is given by (79).

######  Remark 3.7.

Note that the preceding equation implicitly defines gain matrix

(92) $K=C^{v\mathsf{h}}\mathsf{\Gamma}^{-1}.$

This specific choice of gain will play a central role in what follows.
$\blacksquare$

######  Proof 3.8.

(Lemma 3.6) The Gaussian projected filter is defined by evolution of mean and
covariance given by equations (23), (26) and (31), repeated and reordered here
for convenience:

$\displaystyle\widehat{m}_{n+1}$ $\displaystyle=\mathbb{E}\Psi(v_{n}),$
$\displaystyle\widehat{C}_{n+1}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}\Psi(v_{n})-\widehat{m}_{n+1})\otimes\bigl{(}\Psi(v_{n})-\widehat{m}_{n+1}\bigr{)}\Bigr{)}+\Sigma,$
$\displaystyle m_{n+1}$
$\displaystyle=\widehat{m}_{n+1}+\widehat{C}_{n+1}^{vh}(\widehat{C}_{n+1}^{hh}+\Gamma)^{-1}\bigl{(}{y}^{\dagger}_{n+1}-\mathbb{E}\widehat{h}_{n+1}\bigr{)},$
$\displaystyle C_{n+1}$
$\displaystyle=\widehat{C}_{n+1}-\widehat{C}_{n+1}^{vh}(\widehat{C}_{n+1}^{hh}+\Gamma)^{-1}\left({\widehat{C}_{n+1}^{vh}}\right)^{\top},$
$\displaystyle\widehat{C}_{n+1}^{vh}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}\widehat{v}_{n+1}-\mathbb{E}\widehat{v}_{n+1}\bigr{)}\otimes\bigl{(}\widehat{h}_{n+1}-\mathbb{E}\widehat{h}_{n+1}\bigr{)}\Bigr{)},$
$\displaystyle\widehat{C}_{n+1}^{hh}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}\widehat{h}_{n+1}-\mathbb{E}\widehat{h}_{n+1}\bigr{)}\otimes\bigl{(}\widehat{h}_{n+1}-\mathbb{E}\widehat{h}_{n+1}\bigr{)}\Bigr{)}.$

Recall that $\widehat{h}_{n+1}:=h(\widehat{v}_{n+1}).$ Expectations in the
prediction step are with respect to the law of
$v_{n}\sim\mathsf{N}(m_{n},C_{n});$ and expectations in the analysis step are
with respect to the law of $\widehat{v}_{n+1}$ given by (21a), assuming that
$v_{n}\sim\mathsf{N}(m_{n},C_{n});$ thus
$\textit{Law}\,(\widehat{v}_{n+1})=P\textit{Law}\,(v_{n})$.

We now impose the rescalings (72) on these equations. Note that, under the
rescalings and using that $\xi_{n+1}=\mathcal{O}(\Delta t^{\frac{1}{2}})$ has
mean zero,

$\displaystyle\widehat{C}_{n+1}^{vh}$ $\displaystyle=\Delta
t\mathbb{E}\left(\bigl{(}v_{n}-\mathbb{E}v_{n}+\xi_{n}+\mathcal{O}(\Delta
t)\bigr{)}\otimes\bigl{(}\mathsf{h}(v_{n})-\mathbb{E}\mathsf{h}(v_{n})+D\mathsf{h}(v_{n})\xi_{n}+\mathcal{O}(\Delta
t)\bigr{)}\right)$ $\displaystyle=\Delta
t\mathbb{E}\left(\bigl{(}v_{n}-\mathbb{E}v_{n}\bigr{)}\otimes\bigl{(}\mathsf{h}(v_{n})-\mathbb{E}\mathsf{h}(v_{n})\bigr{)}\right)+\mathcal{O}(\Delta
t^{2})$ $\displaystyle=\Delta tC^{v\mathsf{h}}_{n}+\mathcal{O}(\Delta t^{2}).$

Similarly

$\displaystyle\widehat{C}_{n+1}^{hh}+\Gamma$ $\displaystyle=\Delta
t^{2}\mathbb{E}\left(\bigl{(}\mathsf{h}(v_{n})-\mathbb{E}\mathsf{h}(v_{n})+\mathcal{O}(\Delta
t^{\frac{1}{2}})\bigr{)}\otimes\bigl{(}\mathsf{h}(v_{n})-\mathbb{E}\mathsf{h}(v_{n})+\mathcal{O}(\Delta
t^{\frac{1}{2}})\bigr{)}\right)+\Delta t\mathsf{\Gamma}$ $\displaystyle=\Delta
t\mathsf{\Gamma}+\mathcal{O}(\Delta t^{2}).$

Thus

(93a)
$\displaystyle\widehat{C}_{n+1}^{vh}(\widehat{C}_{n+1}^{hh}+\Gamma)^{-1}$
$\displaystyle=C^{v\mathsf{h}}_{n}\mathsf{\Gamma}^{-1}+\mathcal{O}(\Delta t),$
(93b)
$\displaystyle\widehat{C}_{n+1}^{vh}(\widehat{C}_{n+1}^{hh}+\Gamma)^{-1}\widehat{C}_{n+1}^{vh}$
$\displaystyle=\Delta
tC^{v\mathsf{h}}_{n}\mathsf{\Gamma}^{-1}(C^{v\mathsf{h}}_{n})^{\top}+\mathcal{O}(\Delta
t^{2}).$

Furthermore

$\displaystyle\mathbb{E}\Bigl{(}\bigl{(}\Psi(v_{n})$
$\displaystyle-\widehat{m}_{n+1}\bigr{)}\otimes\bigl{(}\Psi(v_{n})-\widehat{m}_{n+1}\bigr{)}\Bigr{)}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}v_{n}+\Delta
tf(v_{n})-\mathbb{E}v_{n}-\Delta
t\mathbb{E}f(v_{n})\bigr{)}\otimes\Bigl{(}\bigl{(}v_{n}+\Delta
tf(v_{n})-\mathbb{E}v_{n}-\Delta t\mathbb{E}f(v_{n})\bigr{)}\Bigr{)}$
$\displaystyle=C_{n}+\Delta tC_{n}^{vf}+\Delta
t(C_{n}^{vf})^{\top}+\mathcal{O}(\Delta t^{2}).$

Using these approximations, and the fact that
$\Delta{z}^{\dagger}_{n+1}=\mathcal{O}(\Delta t^{\frac{1}{2}})$, we obtain

$\displaystyle\widehat{m}_{n+1}$ $\displaystyle=m_{n}+\Delta
t\mathbb{E}f(v_{n}),$ $\displaystyle\widehat{C}_{n+1}$
$\displaystyle=C_{n}+\Delta tC_{n}^{vf}+\Delta t(C_{n}^{vf})^{\top}+\Delta
t\mathsf{\Sigma}+\mathcal{O}(\Delta t^{2}),$ $\displaystyle m_{n+1}$
$\displaystyle=\widehat{m}_{n+1}+C^{v\mathsf{h}}_{n}\mathsf{\Gamma}^{-1}\bigl{(}\Delta{z}^{\dagger}_{n+1}-\Delta
t\mathbb{E}\mathsf{h}(\widehat{v}_{n+1})\bigr{)}+\mathcal{O}(\Delta
t^{\frac{3}{2}}),$ $\displaystyle C_{n+1}$
$\displaystyle=\widehat{C}_{n+1}-\Delta
tC^{v\mathsf{h}}_{n}\mathsf{\Gamma}^{-1}\left(C^{v\mathsf{h}}_{n}\right)^{\top}+\mathcal{O}(\Delta
t^{2}),$ $\displaystyle C^{vf}_{n}$
$\displaystyle=\mathbb{E}\bigl{(}(v_{n}-m_{n})\otimes(f(v_{n})-\mathbb{E}f(v_{n}))\bigr{)},$
$\displaystyle C^{v\mathsf{h}}_{n}$
$\displaystyle=\mathbb{E}\bigl{(}(v_{n}-m_{n})\otimes(\mathsf{h}(v_{n})-\mathbb{E}\mathsf{h}(v_{n}))\bigr{)}.$

In the following $C^{v\mathsf{h}},C^{vf}$ are functions of time, defined by
(90), and $C_{n}^{v\mathsf{h}},C_{n}^{vf}$ are discrete time analogs computed
under the law of $v_{n}.$

Combining the prediction and analysis steps we find that

$\displaystyle m_{n+1}$ $\displaystyle=m_{n}+\Delta
t\mathbb{E}f(v_{n})+C_{n}^{v\mathsf{h}}\mathsf{\Gamma}^{-1}\bigl{(}\Delta{z}^{\dagger}_{n+1}-\Delta
t\mathbb{E}\mathsf{h}(\widehat{v}_{n+1})\bigr{)}+\mathcal{O}(\Delta
t^{\frac{3}{2}})$ $\displaystyle C_{n+1}$ $\displaystyle=C_{n}+\Delta
tC^{vf}_{n}+\Delta t\left(C^{vf}_{n}\right)^{\top}+\Delta
t\mathsf{\Sigma}-\Delta
tC^{v\mathsf{h}}_{n}\mathsf{\Gamma}^{-1}\left(C^{v\mathsf{h}}_{n}\right)^{\top}+\mathcal{O}(\Delta
t^{2}).$

Taking the $\Delta t\rightarrow 0$ limit, we deduce the continuous-time
analogue of the equations (23), (26) as stated.

####  3.4.2 The Kalman-Bucy Filter

We consider the setting where $f(\cdot)=F\cdot$ and
$\mathsf{h}(\cdot)=\mathsf{H}\cdot$ in equations (75) so that

(94a) $\displaystyle dv$ $\displaystyle=Fvdt+\sqrt{\mathsf{\Sigma}}dW,\quad
v_{0}\sim\mathsf{N}(m_{0},C_{0})$ (94b) $\displaystyle dz$
$\displaystyle=Hvdt+\sqrt{\mathsf{\Gamma}}dB,\quad z(0)=0.$

Now consider data ${z}^{\dagger}(t)$ generated by

(95a) $\displaystyle d{v}^{\dagger}$
$\displaystyle=F{v}^{\dagger}dt+\sqrt{\mathsf{\Sigma}}dW^{\dagger},\quad
v_{0}\sim\mathsf{N}(m_{0},C_{0})$ (95b) $\displaystyle d{z}^{\dagger}$
$\displaystyle=H{v}^{\dagger}dt+\sqrt{\mathsf{\Gamma}}dB^{\dagger},\quad
z(0)=0.$

The filtering distribution for $v(t)|{Z}^{\dagger}(t)$ is Gaussian and is
given by the Kalman-Bucy filter, which is the continuous time analog of what
appears in Subsection 2.4.2. It may be obtained from the continuous time
Gaussian projected filter, which is exact in this setting, and also yields a
solution of the Kushner-Stratonovich equation The following summarizes the
result.

######  Theorem 3.9.

Assume that $\mathsf{\Gamma}\succ 0$ in (94). The density $\rho(\cdot,t)$ for
the random variable $v(t)|Z^{\dagger}(t)$, with ${z}^{\dagger}$ given by (95),
is Gaussian with mean $m(\cdot)$ and covariance $C(\cdot)$ solving

(96a) $\displaystyle dm$
$\displaystyle=Fmdt+C\mathsf{H}^{\top}\mathsf{\Gamma}^{-1}\left(dz^{{\dagger}}-\mathsf{H}mdt\right),$
(96b) $\displaystyle dC$
$\displaystyle=FCdt+CF^{\top}dt+\mathsf{\Sigma}dt-C\mathsf{H}^{\top}\mathsf{\Gamma}^{-1}\mathsf{H}Cdt.$

Furthermore, if $\rho|_{t=0}$ is Gaussian then the Kushner-Stratonovich
equation (89) has solution given by the Gaussian
$\mathsf{N}\bigl{(}m(t),C(t)\bigr{)}.$

######  Proof 3.10.

Under (72), (3) the limiting Gaussian projected filter from Theorem 3.6 gives
(96). This equation thus reproduces the true filtering distribution in the
linear-Gaussian setting.

###  3.5 Mean Field Evolution Equations

In Subsection 3.5.1 we describe work concerning the derivation of explicit
mean-field SDEs equal in law to the filtering distribution; this is a
departure from our discussion of this topic in discrete time where no explicit
maps were identified in the general setting. Subsections 3.5.2 and 3.5.3
describe a variety of explicit approximate mean-field SDEs, based on matching
first and second order moment information, and arising from rescaling of the
discrete time setting. We conclude in Subsection 3.5.4 by identifying a mean-
field SDE which is equal in law to the Kalman-Bucy filter.

####  3.5.1 Perfect Transport

Here we seek a mean field dynamical system with law given by that of the
Kushner-Stratonovich equation. The analog of the discrete-time transport map
$T^{S}_{n}$, is to find transport evolution equations for $v$ and
$\widehat{z}$ in the form a mean-field SDE. We will seek to achieve this in
the specific form

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

where $\rho$ is the time-dependent density of $v$. The goal is to chose drift
$a$ and gain $K$ such that the induced time evolution of the density $\rho$ of
$v$ agrees with the density $\rho$ given by the Kushner-Stratonovich equation
(89). We thus find a controlled SDE, similar in form to that proposed in (80)
but with an additional drift term and with mean field dependence. We emphasize
that, as in the discrete-time transport formulation, there is a considerable
degree of non-uniqueness in the choice of transport in general, and here in
the choice of $a$ and $K$ specifically. We make specific, simple, choices in
the theorem that follows. Working in continuous time enables very explicit
identification of exact mean-field models; this is not possible in discrete
time.

######  Theorem 3.11.

Assume that $\mathsf{\Gamma}\succ 0$ and that there exists $K=K(v;\rho)$
satisfying the identity

(98) $\displaystyle-\nabla\cdot(\rho
K^{\top})=\rho\mathsf{\Gamma}^{-1}(\mathsf{h}-\mathbb{E}\mathsf{h}).$

Consider the stochastic mean-field dynamics given by

(99a) $\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}),$ (99b) $\displaystyle d\widehat{z}$
$\displaystyle=\mathsf{h}(v)dt+\sqrt{\mathsf{\Gamma}}dB,$

where ${z}^{\dagger}$ given by (79). Assume that solution $v$ has law with
smooth density and that the Kushner-Stratonovich equation (89) has smooth
density as solution. Then the law of $v$ has density given by the Kushner-
Stratonovich equation (89).

######  Proof 3.12.

To simplify calculations we will first look at choosing $a$ and $K$ to get
agreement, at the level of densities of $v$, with (3.4); a straightforward
modification, using the split-step principle again, then provides the
generalization to (89). This enables us to consider the case where
$f(\cdot)\equiv 0$ and $\mathsf{\Sigma}=0.$

Note that, in (97), ${z}^{\dagger}$ is a fixed, given, trajectory and we are
interested in the evolution of the probability density induced for $v$ by the
randomness over the distribution on trajectories $\widehat{z}$. Using (97b) in
(97a) and setting $f$ and $\mathsf{\Sigma}$ to zero we obtain

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

Applying Fokker-Planck analysis, modified to the mean-field setting, shows
that the time evolution of the density $\rho$ of $v$ under (97) is provided by
the nonlinear (because of dependence of $a,K$ on $\rho$) SPDE121212Here we use
the standard convention from continuum mechanics that the divergence of a
matrix is to be interpreted via computation of derivatives with respect to the
second index; see Subsection 3.7 for references to relevant continuum
mechanics textbooks. Strictly speaking equation (101) is an SPDE only if we
now view ${z}^{\dagger}$ as a random variable rather than a fixed realization.
to be interpreted in the Itô sense:

(101) $\displaystyle d\rho$
$\displaystyle=-\nabla\cdot(\rho(a-K\mathsf{h}))dt-\langle\nabla\cdot(\rho
K^{\top}),d{z}^{\dagger}\rangle+\nabla\cdot(\nabla\cdot(\rho
K\mathsf{\Gamma}K^{\top}))dt.$

Note that, although ${z}^{\dagger}$ is a fixed trajectory, it contributes to
the diffusion term in this Fokker-Planck equation, explaining the factor $1$
rather than $1/2.$ This arises as a contribution from the quadratic variation
of the path of ${z}^{\dagger}$ to the evolution of $\rho.$ For further details
on the derivation of (101) see Lemma 9.3 in Appendix 9. Comparing with (3.4)
we find that to get agreement $a$ and $K$ have to be chosen such that

$\displaystyle\rho\langle(\mathsf{h}-\mathbb{E}\mathsf{h}),\mathsf{\Gamma}^{-1}(d{z}^{\dagger}-\mathbb{E}\mathsf{h}dt)\rangle=$
$\displaystyle\qquad\qquad-\nabla\cdot(\rho(a-K\mathsf{h}))dt-\langle\nabla\cdot(\rho
K^{\top}),d{z}^{\dagger}\rangle+\nabla\cdot(\nabla\cdot(\rho
K\mathsf{\Gamma}K^{\top}))dt.$

Equating the two terms involving the data ${z}^{\dagger}$ shows immediately
that $K(v;\rho)$ has to satisfy the (vector-valued) PDE (98). From this,
equating the terms that do not involve the data ${z}^{\dagger}$, it follows
that $a(v;\rho)$ has to satisfy

$\displaystyle\nabla\cdot(\rho a)$ $\displaystyle=\nabla\cdot(\nabla\cdot(\rho
K\mathsf{\Gamma}K^{\top}))+\nabla\cdot(\rho
K(\mathsf{h}-\mathbb{E}\mathsf{h}))$
$\displaystyle=\nabla\cdot(\nabla\cdot(\rho
K\mathsf{\Gamma}K^{\top}))-\nabla\cdot(K\mathsf{\Gamma}\nabla\cdot(\rho
K^{\top}))$
$\displaystyle=\nabla\cdot\left(\rho\left(\nabla\cdot(K\mathsf{\Gamma}K^{\top})-K\mathsf{\Gamma}\nabla\cdot
K^{\top}\right)\right).$

A natural choice for $a$ is provided by asking that the term on which the
divergence acts is zero. This yields

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

a solution for $a$ with no explicit dependence on $\rho;$ note, however, that
$a$ does depend on $\rho$ implicitly through the dependence of $K$ on $\rho.$

With these choices of $a,K$, and applying the split-step principle so that the
mean-field dynamics are consistent with (89) rather than (3.4), we obtain a
version of the feedback particle filter; in particular we find the equation in
its stochastic mean-field formulation given by (99) where $K=K(v;\rho)$ solves
(98) and $\rho$ evolves according to the Kushner-Stratonovich equation (89),
an equation which also defines the law of $v.$

######  Remark 3.13.

There is an interesting interpretation of the contribution $a$ to the drift
term in (99): it is simply the Itô-to-Stratonovich-like correction with regard
to the $v$-dependence in $K(v;\rho)$; however there is a subtlety that the
equation for $\rho$ itself depends on the data ${z}^{\dagger}$ and the
correction does not account for the $\rho-$dependence of $K(v;\rho)$ – the
Stratonovich correction is only with respect to $v-$dependence of the drift;
an Itô interpretation is retained with respect to the $\rho$ dependence. We
refer to Appendix 9.2 for an in-depth discussion concerning this unusual form
of stochastic integration. In that section we also derive the full
Stratonovich correction (with respect to both $v$ and $\rho$ dependence) of
the exact mean-field model. $\blacksquare$

######  Remark 3.14.

The mean field equations (99) require a gain $K(v;\rho)$ which satisfies (98).
Let $\mathbb{E}$ denote expectation with respect to random variable $v$
distributed according to probability measure with density $\rho$. Using
appropriately regular test functions
$\psi:\mathbb{R}^{d_{v}}\to\mathbb{R}^{d_{v}}$, chosen to have mean-zero under
$\mathbb{E}$, equation (98) can be rephrased in the weak form 131313 Recall
that the conventions from continuum mechanics that we use to define the
divergence and gradient of vector fields are discussed in the bibliography
Subsection 3.7.

(103)
$\mathbb{E}\left(K^{\top}\nabla\psi\right)=\mathsf{\Gamma}^{-1}C^{\mathsf{h}\psi};$

here $C^{\mathsf{h}\psi}$ denotes the covariance between $\mathsf{h}$ under
$\mathbb{E}$ and, in what follows, we also let $C^{v\mathsf{h}}$ denote
covariance between $v$ and $\mathsf{h}$.

The particular choice $\psi(v)=v-\mathbb{E}v$ leads to

$\mathbb{E}K=C^{v\mathsf{h}}\mathsf{\Gamma}^{-1}.$

This relation will be central to the subsequent discussion concerning second
order transport methods in continuous time. Note, in particular, that it can
be satisfied by making the _constant gain_ ansatz that $K$ is independent of
$v$ and is then given by $K=C^{v\mathsf{h}}\mathsf{\Gamma}^{-1},$ as derived
for the Gaussian projected filter in Remark 3.7. More generally speaking, we
note that (103) is amenable to numerical approximations. This will be
discussed further in Subsection 3.6.1 below. $\blacksquare$

We close this subsection with the mean field formulation using deterministic
innovations:

(104)
$dv=f(v)dt+\sqrt{\mathsf{\Sigma}}dW+\nabla\cdot(K\mathsf{\Gamma}K^{\top})dt-K\mathsf{\Gamma}\nabla\cdot
K^{\top}dt+K\left(d{z}^{\dagger}-\frac{1}{2}(\mathsf{h}+\mathbb{E}\mathsf{h})dt\right),$

where ${z}^{\dagger}$ is given by (79). The relationship between this equation
and (99), is analogous to the relationship between innovation terms (51) and
(15) already encountered in the discrete time setting.

####  3.5.2 Second Order Transport – Stochastic Case

Recall that, in the discrete-time setting, there are an uncountable set of
maps which effect approximate transport, in the sense of matching the first
and second order statistics of the analysis map, using either stochastic or
deterministic models. These are elucidated in Appendix 8. We, however,
concentrated in the main text on a handful of examples. In this and the next
subsection we study continuous time analogs of some of these examples,
starting with the Kalman transport map (46) recalled here, and reformulated,
for convenience:141414The notation for the first two equations, which we use
variants of in what follows, is shorthand for the first two equations
appearing in (14), together with the assumptions detailed following those
equations. Thus we use $\mathsf{N}(0,\Sigma)$ to denote an i.i.d. realization
from the stated Gaussian distribution, and similarly for other variables.

(105a) $\displaystyle\widehat{v}_{n+1}$
$\displaystyle=\Psi(v_{n})+\mathsf{N}(0,\Sigma),$ (105b)
$\displaystyle\widehat{y}_{n+1}$
$\displaystyle=h(\widehat{v}_{n+1})+\mathsf{N}(0,\Gamma),$ (105c)
$\displaystyle v_{n+1}$
$\displaystyle=\widehat{v}_{n+1}+\widehat{C}^{vh}_{n+1}(\widehat{C}^{hh}_{n+1}+\Gamma)^{-1}({y}^{\dagger}_{n+1}-\widehat{y}_{n+1}).$

We now apply the rescaling (72) to obtain, using (93a),

$\displaystyle\widehat{v}_{n+1}$ $\displaystyle=v_{n}+\Delta
tf(v_{n})+\sqrt{\Delta t}\mathsf{N}(0,\mathsf{\Sigma}),$
$\displaystyle\widehat{z}_{n+1}$ $\displaystyle=\widehat{z}_{n}+\Delta
t\mathsf{h}(\widehat{v}_{n+1})+\sqrt{\Delta t}\mathsf{N}(0,\mathsf{\Gamma}),$
$\displaystyle v_{n+1}$
$\displaystyle=\widehat{v}_{n+1}+C^{v\mathsf{h}}_{n}\mathsf{\Gamma}^{-1}(\Delta{z}^{\dagger}_{n+1}-\Delta\widehat{z}_{n+1})+\mathcal{O}(\Delta
t^{3/2}).$

The preceding calculation, leading to $\mathcal{O}(\Delta t^{3/2})$ error,
uses the fact that the noises entering the equations for $\widehat{v}_{n+1},$
$\widehat{z}_{n+1}$ and $z^{\dagger}_{n+1}$ are independent. Taking the
continuous time limit we obtain

(106a) $\displaystyle dv$
$\displaystyle=f(v)dt+\sqrt{\mathsf{\Sigma}}dW+C^{v\mathsf{h}}\mathsf{\Gamma}^{-1}(d{z}^{\dagger}-d\widehat{z}),$
(106b) $\displaystyle d\widehat{z}$
$\displaystyle=\mathsf{h}(v)dt+\sqrt{\mathsf{\Gamma}}dB.$

Again $W,B$ are independent unit Brownian motions of appropriate dimensions
and ${z}^{\dagger}$ is given by (79).

######  Remark 3.15.

Note the recurrence of the continuous time gain $K$ first identified in Remark
3.7. We have derived the mean field equations (106) from the associated
stochastic discrete time formulation (105). However there is another way to
derive them. If we make the assumption that the gain $K=K(\rho)$ is
independent of $v$ then the equations (106) follow directly from the general
framework (99) if (3.14) is invoked. This is because both divergence terms
drop out in (99) if $K$ is independent of $v$, and (3.14) simply gives
$K=C^{v\mathsf{h}}\mathsf{\Gamma}^{-1}$. $\blacksquare$

####  3.5.3 Second Order Transport – Deterministic Case

We may apply a similar analysis to (46), rewritten and reformulated here for
convenience:

$\displaystyle\widehat{v}_{n+1}$
$\displaystyle=\Psi(v_{n})+\mathsf{N}(0,\Sigma),\quad\widehat{h}_{n+1}=h(\widehat{v}_{n+1}),$
$\displaystyle v_{n+1}$
$\displaystyle=\widehat{v}_{n+1}-\widetilde{K}_{n}(\widehat{h}_{n+1}-\mathbb{E}\widehat{h}_{n+1})+K_{n}({y}^{\dagger}_{n+1}-\mathbb{E}\widehat{h}_{n+1}),$

where

$\displaystyle K_{n}$
$\displaystyle=\widehat{C}_{n+1}^{vh}\left(\widehat{C}_{n+1}^{hh}+\Gamma\right)^{-1},\quad\widetilde{K}_{n}=\widehat{C}^{vh}_{n+1}\Bigl{(}(\widehat{C}^{hh}_{n+1}+\Gamma)+\Gamma^{1/2}(\widehat{C}^{hh}_{n+1}+\Gamma)^{1/2}\Bigr{)}^{-1}.$

Applying the rescalings (72) gives

$\displaystyle\widehat{v}_{n+1}$ $\displaystyle=v_{n}+\Delta
tf(v_{n})+\sqrt{\Delta t}\mathsf{N}(0,\mathsf{\Sigma}),$ $\displaystyle
v_{n+1}$
$\displaystyle=\widehat{v}_{n+1}+C_{n}^{v\mathsf{h}}\mathsf{\Gamma}^{-1}\Bigl{(}\Delta{z}^{\dagger}_{n+1}-\frac{1}{2}(\mathsf{h}(\widehat{v}_{n+1})+\mathbb{E}\mathsf{h}(\widehat{v}_{n+1})\Bigr{)}+{\mathcal{O}}(\Delta
t^{2}),$

where we have used

$\displaystyle K_{n}$
$\displaystyle=C_{n}^{v\mathsf{h}}\mathsf{\Gamma}^{-1}+\mathcal{O}(\Delta
t^{3/2}),\quad\widetilde{K}_{n}=\frac{1}{2}C_{n}^{v\mathsf{h}}\mathsf{\Gamma}^{-1}+\mathcal{O}(\Delta
t^{3/2}).$

Taking the continuous time limit we obtain

(107)
$dv=f(v)dt+\sqrt{\Sigma}dW+C^{v\mathsf{h}}\mathsf{\Gamma}^{-1}\Bigl{(}d{z}^{\dagger}-\frac{1}{2}\bigl{(}\mathsf{h}(v)+\mathbb{E}\mathsf{h}(v)\bigr{)}dt\Bigr{)},$

where, once again, ${z}^{\dagger}$ is given by (79).

These mean field equations can also be derived directly from the general mean
field equation (104), assuming a deterministic innovation, invoking the
constant gain approximation $K=K(\rho)$ and using (3.14). The reasoning is
similar to that in Remark 3.15.

####  3.5.4 Mean Field Formulation of the Kalman-Bucy Filter

######  Theorem 3.16.

Assume that $\mathsf{\Gamma}\succ 0$ and consider the evolution of random
variable $v$ initialized at a Gaussian and satisfying

(108a) $\displaystyle dv$
$\displaystyle=Fvdt+\sqrt{\Sigma}dW+CH^{\top}\mathsf{\Gamma}^{-1}(d{z}^{\dagger}-d\widehat{z}),$
(108b) $\displaystyle d\widehat{z}$ $\displaystyle=Hvdt+\sqrt{\Gamma}dB,$

where $C$ is the covariance of $v$ and where ${z}^{\dagger}$ is given by (95).
Then the mean and covariance of $v$ satisfy equations (96).

######  Proof 3.17.

We consider the setting where $f(\cdot)=F\cdot$ and $\mathsf{h}(\cdot)=H\cdot$
with $F$ and $H$ as defined in (72). The mean field models become exact in
this setting and have mean and covariance satisfying (96). With this in mind
we note that under the linearity assumptions (106) gives (108).

Similar considerations, starting at (107), give the following variant of the
preceding theorem:

######  Theorem 3.18.

Assume that $\mathsf{\Gamma}\succ 0$ and consider the evolution of random
variable $v$ initialized at a Gaussian and satisfying

(109)
$dv=Fvdt+\sqrt{\Sigma}dW+CH^{\top}\mathsf{\Gamma}^{-1}\Bigl{(}d{z}^{\dagger}-\frac{1}{2}H\bigl{(}v+m\bigr{)}dt\Bigr{)},$

where $(m,C)$ are the mean and covariance of $v$ and where ${z}^{\dagger}$ is
given by (95). Then the mean and covariance of $v$ satisfy equations (96).

###  3.6 Ensemble Kalman Methods

We now discuss several particle approximations to the mean field equations
derived in the preceding subsections. We start with the mean field equations
(99), based on perfect transport, before considering particle approximations
for time continuous approximate transport formulations. Throughout this
subsection ${z}^{\dagger}$ is given by (95).

####  3.6.1 Perfect Particle Filters

As in Subsection 2.6.1 we define $\mathsf{J}:=\\{1,\cdots,J\\}.$ The desired
approximation to the mean field model (99) evolves particle ensemble
$\\{{v}^{(j)}\\}_{j\in\mathsf{J}}$ and its associated empirical measure

(110) $\mu^{\mathsf{J}}(t)=\frac{1}{J}\sum_{j=1}^{J}\delta_{{v}^{(j)}}$

according to the interacting particle system

(111a) $\displaystyle d{v}^{(j)}$
$\displaystyle=f({v}^{(j)})dt+a^{(j)}dt+\sqrt{\Sigma}d{W}^{(j)}+K^{(j)}(d{z}^{\dagger}-d\widehat{z}^{(j)}),$
(111b) $\displaystyle d\widehat{z}^{(j)}$
$\displaystyle=\mathsf{h}({v}^{(j)})dt+\sqrt{\Gamma}d{B}^{(j)},$

where ${z}^{\dagger}$ is given by (79) and the
$\\{{W}^{(j)}\\}_{j\in\mathsf{J}}$ and $\\{{B}^{(j)}\\}_{j\in\mathsf{J}}$ are
mutually independent collections of i.i.d. Brownian motions in
$\mathbb{R}^{d_{v}}$ and $\mathbb{R}^{d_{y}}$ respectively. The interaction
between the particles arises from the gain matrices

(112) $K^{(j)}:=K^{\mathsf{J}}({v}^{(j)};\mu^{\mathsf{J}}),$

$j\in\mathsf{J}$, where the matrix valued function $K^{\mathsf{J}}$
approximates the solution $K$ of (98); the drift term $a^{(j)}$ is defined by
(102) with $K(v)$ replaced by $K^{\mathsf{J}}(v;\mu^{\mathsf{J}})$ and then
evaluated at $v=v^{(j)}$. To be self-consistent we require that the empirical
measure (110) becomes the desired particle approximation for the law of $v$
satisfying (99).

The key challenge for implementing a particle approximation of the mean field
model (99) lies in the numerical approximation of (98), or its weak
formulation (103). We start with the numerical approximation of the weak
formulation (103). Let $\mathcal{F}$ denote a space of vector-valued functions
mapping $\mathbb{R}^{d_{v}}$ into $\mathbb{R}^{d}$ and assumed to have mean
zero with respect to expectation $\mathbb{E}$ under density $\rho$.151515We
will consider different choices for $d$ but use the same notation for the
space. We make the ansatz that $K^{\top}=\nabla\Psi$, for some
$\Psi\in\mathcal{F}.$ Then (103) can be rewritten as

$\mathbb{E}(\nabla\Psi\nabla\psi)=\mathsf{\Gamma}^{-1}C^{\mathsf{h}\psi},$

assumed to hold for all $\psi\in\mathcal{F}.$ In order to obtain the desired
numerical approximation $K^{\mathsf{J}}$ let $\mathbb{E}^{\mathsf{J}}$ denote
expectation with respect to the empirical measure (110) and consider
$(K^{\mathsf{J}})^{\top}=\nabla\Psi^{\mathsf{J}}$ and

(113)
$\mathbb{E}^{\mathsf{J}}(\nabla\Psi^{\mathsf{J}}\nabla\psi)=\mathsf{\Gamma}^{-1}C^{\mathsf{h}\psi}$

with the correlation $C^{\mathsf{h}\psi}$ approximated by

$C^{\mathsf{h}\psi}=\mathbb{E}^{\mathsf{J}}\Bigl{(}\bigl{(}\mathsf{h}(v)-\mathbb{E}^{\mathsf{J}}\mathsf{h}(v)\bigr{)}\otimes\bigl{(}\psi(v)-\mathbb{E}^{\mathsf{J}}\psi(v)\bigr{)}\Bigr{)}.$

The final step in the numerical approximation is to choose an appropriate
finite dimensional subspace $\mathcal{F}^{L}\subset\mathcal{F}$ of dimension
$L\ll J$ and to determine $\Psi^{\mathsf{J}}\in\mathcal{F}^{L}$ such that
(113) holds for all $\psi\in\mathcal{F}^{L}$.

We return to the strong formulation (98) again making the ansatz
$K^{\top}=\nabla\Psi$ and giving rise to the equivalent formulation

$\mathcal{L}_{\rho}\Psi=-\mathsf{\Gamma}^{-1}(\mathsf{h}-\mathbb{E}\mathsf{h})$

with differential operator $\mathcal{L}_{\rho}$ defined by

$\mathcal{L}_{\rho}\Psi=\rho^{-1}\nabla\cdot(\rho\nabla\Psi)=\nabla^{2}\Psi+\nabla\Psi\nabla\log\rho.$

Note that in the case of scalar observations, $\mathcal{L}_{\rho}$ is the
infinitesimal generator of a diffusion process with invariant density $\rho$;
in the vector case the same statement holds component-by-components. This fact
may be used to approximate the action of $\mathcal{L}_{\rho}$ via its heat
semigroup, and this may be used as the basis for numerical approximations. See
the bibliography Subsection 3.7 for references to the relevant literature.

While we have discussed two general approaches to approximate the gain matrix
in (112) numerically so far, we now return to the the constant gain
approximation; that is, we make $K^{\mathsf{J}}(\cdot;\mu^{\mathsf{J}})$
independent of its first argument so that $K^{(j)}=K^{\mathsf{J}}$ for some
$K^{\mathsf{J}}$ depending only on $\mu^{\mathsf{J}}.$ In doing so we recover
an empirical version of (92), the continuous time gain first identified in
discussion of the Gaussian projected filter in Remark 3.7, and arising again
in second order transport as discussed in Remark 3.15. In the current context
the identification of $K^{\mathsf{J}}$ arises from (113), by making the choice
$\Psi(v)=(K^{\mathsf{J}})^{\top}(v-\mathbb{E}^{\mathsf{J}}v)$ and
$\psi(v)=v-\mathbb{E}^{\mathsf{J}}v$. These choices result in

(114) $K^{\mathsf{J}}=C^{v\mathsf{h}}\mathsf{\Gamma}^{-1},$

where

$C^{v\mathsf{h}}=\mathbb{E}^{\mathsf{J}}\Bigl{(}\bigl{(}v-\mathbb{E}^{\mathsf{J}}v\bigr{)}\otimes\bigl{(}\mathsf{h}(v)-\mathbb{E}^{\mathsf{J}}\mathsf{h}(v)\bigr{)}\Bigr{)}.$

Note that the constant gain approximation (114) also implies that the drift
term $a^{(j)}$ in (111) vanishes. We summarize numerical implementation
details in the following two subsections.

####  3.6.2 Stochastic Ensemble Kalman Filters

We now consider the mean-field model (106) and its numerical approximation.
Since the drift correction $a$ does not appear here, the only significant
difference to the interacting particle approximation (111) arises from the
choice of the interaction term, that is, $K^{(j)}$, which is now independent
of $j$ and determined by (114). In summary, we obtain the following

(115a) $\displaystyle d{v}^{(j)}$
$\displaystyle=f({v}^{(j)})dt+\sqrt{\Sigma}d{W}^{(j)}+C^{v\mathsf{h}}\mathsf{\Gamma}^{-1}(d{z}^{\dagger}-d\widehat{z}^{(j)}),$
(115b) $\displaystyle d\widehat{z}^{(j)}$
$\displaystyle=\mathsf{h}({v}^{(j)})dt+\sqrt{\Gamma}d{B}^{(j)},$

####  3.6.3 Deterministic Ensemble Kalman Filters

We make an empirical approximation of the mean-field model (107) as follows:
for $j\in\mathsf{J}:=\\{1,\dots,J\\}$ we consider

(116)
$d{v}^{(j)}=f({v}^{(j)})dt+\sqrt{\Sigma}d{W}^{(j)}+C^{v\mathsf{h}}\mathsf{\Gamma}^{-1}\Bigl{(}d{z}^{\dagger}-\frac{1}{2}\bigl{(}\mathsf{h}({v}^{(j)})+\mathbb{E}^{\mathsf{J}}\mathsf{h}(v)\bigr{)}dt\Bigr{)}.$

The notation is as in the preceding subsection and, in particular, the formula
(110) again gives the particle approximation of the approximate filter; now
the relevant approximate filter is defined by the distribution of (116).

###  3.7 Bibliographical Notes

This chapter is framed in the language of SDEs; see Evans (2012), Oksendal
(2013) for background in this area. In passing from discrete to continuous
time we often invoke ideas from the numerical solution of SDEs; see Kloeden &
Platen (1991) and Higham (2001) for introductions to this area. The
conventions from continuum mechanics, that we use to define the divergence and
gradient of vector fields, are the same as those adopted, and described in
detail, in Gonzalez & Stuart (2008), Gurtin (1982).

The Kalman-Bucy filter (Kalman & Bucy, 1961) contains what is perhaps the
first systematic derivation and analysis of an algorithm for the incorporation
of continuous time data into estimation of a sample path of an SDE. Its
extension to nonlinear and non-Gaussian distributions is provided by the
Kushner-Stratonovich equation (89). An heuristic derivation of both the
Kalman-Bucy filter as well as the Kushner-Stratonovich equation can be found
in Jazwinski (2007) while Bain & Crisan (2008) covers the field of continuous
time filtering in full detail.

The idea of Strang-splitting, which we use to derive the Kushner-Stratonovich
equation, originates in Strang (1968); for an overview of splitting methods
see McLachlan & Quispel (2002). Furthermore, we rely on robustness results for
continuous time filtering (Clark & Crisan, 2005), which imply that smooth
approximations $z^{\dagger,\delta}$ to stochastic observations ${z}^{\dagger}$
are justified and that the order of taking limits $\Delta t\to 0$ and
$\delta\to 0$ can be accounted for by appropriate Stratonovich to Itô
correction terms. An introduction to the required covariation formulae used in
proof of Lemma 3.3 can be found in Eberle (2013). We note that robustness
results do not carry over to associated mean field equations and filtering
problems with correlated noise (Coghi et al., 2021).

A mean field approach to the Kushner-Stratonovich equation (89) appeared first
in the work of Crisan & Xiong (2010), which utilizes robustness results and
smoothed data $z^{\dagger,\delta}$ in the $\delta\to 0$ limit. The mean field
equations (107) were proposed in Yang et al. (2013) while the stochastic
counterpart appeared first in Reich (2019). The mathematical relationship
between the various mean field formulations has been analyzed in Pathiraja et
al. (2021). The continuous time ensemble Kalman filter formulations (115) and
(116) appeared first in Bergemann & Reich (2012). These formulations are based
on earlier work on homotopy formulations of the Bayesian inference step by
Bergemann & Reich (2010b) and Reich (2011). Rigorous derivation of continuous
time ensemble Kalman filter formulations from their discrete time counterparts
can be found in Lange & Stannat (2019), Lange & Stannat (2021), Blömker et al.
(2018) and Blömker et al. (2021). Derivation and analysis of the properties of
continuous time limits in the context of solving inverse problems may be found
in Schillings & Stuart (2017). Numerical time-stepping methods for continuous
time ensemble Kalman filter formulations are analyzed in Amezcua et al.
(2014).

The numerical approximation of the Kushner-Stratonovich equation (89) has a
long history. The paper Crisan & Lyons (1999) demonstrated how (89) can be
approximated by a particle method; this is a generalization of the bootstrap
particle filter to continuous time Crisan et al. (1998). See also Bain &
Crisan (2008) for a detailed discussion of alternative approximation
techniques. The paper Hu et al. (2002) discussed solution of the unnormalized
and linear version of the Kushner-Stratonovich equation known as the Zakai
equation. Numerical implementations of the mean field equations (99) are
discussed in Taghvaei et al. (2017) while Taghvaei et al. (2020) provides a
detailed analysis of the diffusion map approximation discussed after (3.6.1).
The constant gain approximation $K=C^{v\mathsf{h}}\mathsf{\Gamma}$, which
corresponds to the ensemble Kalman filter, arises as a particular scaling
limit from the diffusion map approach as discussed, for example, in Taghvaei
et al. (2017). 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 in continuous time, mostly focused on the solution of inverse
problems; however the methods developed are more widely applicable. Similar to
the stochastic ensemble Kalman filter, particle implementation (115) leads to
undesirable correlations via the synthetic data $\widehat{z}^{(j)}$, which
appear both in the gain $K^{\mathsf{J}}=C^{v\mathsf{h}}\mathsf{\Gamma}^{-1}$
and the innovation term $dI=d{z}^{\dagger}-d\widehat{z}^{(j)}$ effectively
giving rise to colored noise. These numerically induced correlations vanish in
the limit $J\to\infty$.

Theoretical analysis of the continuous time 3DVAR algorithm (78) may be found
in Law et al. (2012), Blömker et al. (2013), Azouani et al. (2014), Gesho et
al. (2016), Olson & Titi (2003), Mondaini & Titi (2018) and Larios & Pei
(2017). Well-posedness, stability, and accuracy results for the ensemble
Kalman filter have been derived, for example, in Kelly et al. (2014), de
Wiljes et al. (2018), Del Moral & Tugaut (2018), and de Wiljes & Tong (2020).

##  4 Inverse Problems

Inverse problems may be solved by using ideas from filtering in stochastic
dynamical systems, by the application of ideas arising in ensemble Kalman
methods, and by the conjunction of the two. The power of the resulting methods
is that they are derivative-free, requiring the forward model only as a black
box. After commencing, in Subsection 4.1, with definitions of the classical
(optimization) and Bayesian (probabilistic) perspectives on the inverse
problems, we proceed to study three approaches to inversion based on ideas
from filtering and ensemble Kalman methodology. These are:

  * •

In Subsection 4.2 we show how the Bayesian formulation of the inverse problem
may be tackled by introducing a transport from prior to posterior, underpinned
by a stochastic dynamical system with filtering distribution given by the
posterior, at step $N$, when initiated at the prior, at step $0$. Algorithms
to approximate the filtering distribution may then be employed, leading to a
variety of approaches to Bayesian inversion, and in particular to ensemble
Kalman based filtering approaches to inversion.

  * •

Subsection 4.3 establishes a connection between the Kalman type mean field
equations from Subsection 4.2 and gradient descent methods used in
optimization. In particular, continuous time Kalman mean field equations
naturally lead to preconditioned mean field gradient descent formulations
which enjoy the property of affine invariance. Furthermore, it is demonstrated
that these mean field equations possess a gradient flow structure in the space
of probability measures. This leads to a generalization of the well-known
Wasserstein gradient flow structure, which we refer to as a _Kalman-
Wasserstein_ gradient flow. While these structures rely on an exact
computation of gradients, an attractive feature of Kalman type gradient
descent methods lies in their ability to approximate gradients using only
evaluations of the forward problem, and no adjoint calculations. Such
approximations can be justified via statistical linearization and are now used
in a wide range of applications.

  * •

Subsection 4.4 is based on finding partially and noisily observed stochastic
dynamical systems with filtering distribution that is ergodic and delivers a
solution of the inverse problem in statistical steady state; ensemble Kalman
methods may then be iterated to find this statistical steady state. The
starting point is a stochastic dynamical system for which the steady state of
the filtering distribution is equal, for linear Gaussian inverse problems, to
a Gaussian with mean given by a regularized least squares solution to the
inverse problem, or with Gaussian equal to the posterior distribution. This
property does not hold exactly in the nonlinear case, but the assumptions
under which it does hold, namely linear and Gaussian structure, are precisely
those under which the ensemble Kalman second-order transport filters are
exact. It is thus self-consistent to apply ensemble Kalman methods to filter
the stochastic dynamical system to solve problems beyond the linear Gaussian
setting, by studying large time asymptotic behavior of the filtering
distribution.

Throughout these three subsections we present ideas based on the road-map of
Figure 1. Furthermore in each subsection we study both discrete and continuous
time; the latter is useful for insight whilst the former provides the basis
for practical algorithms. Since Sections 2 and 3 show in detail how to derive
a variety of ensemble Kalman methods from mean-field models based on second
order transport, the focus of the work in this section is on: a) the
identification of appropriate dynamical systems with filtering distribution
linked to the inverse problem; b) the properties of mean-field models derived
from the filtering problem arising from the dynamical systems; c) the
derivation of other novel mean-field models arising from the concept of
statistical linearization. Once mean-field models have been identified in b)
and c), application of ideas from Sections 2 and 3 is straightforward.
However, in Subsection 4.5 we present applications of ensemble Kalman
methodology to some specific mean-field models identified for inversion; and
we include a simple illustrative example, building on the Lorenz ’96 model
considered earlier, to show the methods is use for inverse problems arising in
parameter estimation for dynamical systems. We conclude in Subsection 4.6,
containing bibliographic notes.

###  4.1 Set-up

Consider the inverse problem of finding unknown parameter
$u\in\mathbb{R}^{d_{u}}$ from data $w\in\mathbb{R}^{d_{w}}$, when $w$ is
related to $u$ via the equation

(117) $w=G(u)+\gamma.$

Here $G:\mathbb{R}^{d_{u}}\to\mathbb{R}^{d_{w}}$ is the _forward model_ and
$\gamma$ represents noise polluting the data. We assume that $G$ is bounded on
compact subsets of $\mathbb{R}^{d_{u}}.$

We consider the _Bayesian_ approach to this inverse problem. To be concrete we
assume as _prior_ $u\sim\mathsf{N}(m_{0},C_{0})$, that
$\gamma\sim\mathsf{N}(0,\mathsf{\Gamma})$ and that $u$ and $\gamma$ are
independent; thus the likelihood
$w|u\sim\mathsf{N}\bigl{(}G(u),\mathsf{\Gamma}\bigr{)}$. Throughout this
section we assume that $C_{0}\succ 0$ and $\mathsf{\Gamma}\succ 0.$ We may
then define

$\displaystyle\Phi(u)$
$\displaystyle=\frac{1}{2}|w-G(u)|_{\mathsf{\Gamma}}^{2},$
$\displaystyle\Phi_{R}(u)$
$\displaystyle=\Phi(u)+\frac{1}{2}|u-m_{0}|_{C_{0}}^{2}.$

Under the stated assumptions, application of Bayes Theorem shows that the
_posterior distribution_ on $u|w$ is measure $\mu^{w}$ given by

(118a) $\displaystyle\mu^{w}(du)$ $\displaystyle=\frac{1}{\cal
Z}\exp\bigl{(}-\Phi_{R}(u)\bigr{)}du,$ (118b) $\displaystyle\cal Z$
$\displaystyle=\int_{\mathbb{R}^{d_{u}}}\exp\bigl{(}-\Phi_{R}(u)\bigr{)}du.$

We will also discuss least-squares based optimization approaches to the
inverse problem. These revolve around minimization of $\Phi$ or $\Phi_{R}$.
Minimizing $\Phi$ is a _maximum likelihood approach_ ; minimizing $\Phi_{R}$
is a _maximum a posteriori approach_. In the latter context it is useful to
observe that we may write

$\Phi_{R}(u)=\frac{1}{2}|w_{R}-G_{R}(u)|^{2}_{\mathsf{\Gamma}_{R}},$

where

(119) $w_{R}:=\begin{pmatrix}w\\\ m_{0}\end{pmatrix},\quad\quad\quad
G_{R}(u):=\begin{pmatrix}G(u)\\\
u\end{pmatrix},\quad\quad\quad\mathsf{\Gamma}_{R}:=\begin{pmatrix}\mathsf{\Gamma}&0\\\
0&C_{0}\end{pmatrix}.$

Thus the regularized least squares problem associated with the inverse problem
(117) may be viewed as an unregularized least squares problem arising from the
modified inverse problem

$w_{R}=G_{R}(u)+\gamma_{R},$

with $\gamma_{R}\sim\mathsf{N}(0,\mathsf{\Gamma}_{R}).$

In the case of linear $G(\cdot)=L\cdot$ we define

(120) $L_{R}:=\begin{pmatrix}L\\\ I\end{pmatrix}$

and note that, with this definition, $G_{R}(\cdot)=L_{R}\cdot.$ Note also
that, since $C_{0}\succ 0$, $\Phi_{R}$ has a unique minimizer if $G$ is
linear. Since, in this linear setting,

$\Phi_{R}(u)=\frac{1}{2}|w_{R}-L_{R}u|^{2}_{\mathsf{\Gamma}_{R}},$

we deduce that the posterior covariance $C_{R}$ of the Gaussian posterior
$\mu^{w}$ is given via its precision, the Hessian of $\Phi_{R}:$

(121) $C_{R}^{-1}=L_{R}^{\top}\mathsf{\Gamma}_{R}^{-1}L_{R}.$

The equilibrium mean $m_{R}$ solves the _normal equations_

(122) $C_{R}^{-1}m_{R}=L_{R}^{\top}\mathsf{\Gamma}_{R}^{-1}w_{R}.$

The following subsections are devoted to the derivation and study of mean
field methods, which are therefore the basis for ensemble Kalman methods, to
sample from $\mu^{w}$, or to study the related least squares optimization
problems of minimizing $\Phi(\cdot)$ or $\Phi_{R}(\cdot)$.

###  4.2 Transporting Prior To Posterior

The basic idea used in this subsection, to address the solution of inverse
problems, is rooted in a sequential formulation of Bayesian inference. To
understand this, define for integer $N>1$,

(123)
$\Phi_{R,n}(u)=\frac{1}{N}\sum_{j=1}^{n}\Phi(u)+\frac{1}{2}|u-m_{0}|_{C_{0}}^{2},$

noting that $\Phi_{R,N}(u)=\Phi_{R}(u).$ Now consider the sequence $\mu_{n}$
of probability measures with negative log density given (up to an additive
constant with respect to variation of $u$) by $\Phi_{R,n}(u).$ Then the
Bayesian inference problem (118) can be reformulated as a sequence of $N$
Bayesian inference steps where the prior $\mu_{n}$ is morphed into posterior
$\mu_{n+1}$ using the data likelihood $\exp(-\Phi(u)/N)$ in each step:

$\mu_{n+1}(du)\propto\exp\Bigl{(}-\frac{1}{N}\Phi(u)\Bigr{)}\mu_{n}(du).$

The initial prior is set to $\mu_{0}=\mathsf{N}(m_{0},C_{0})$ and the $N^{th}$
posterior $\mu_{N}$ delivers the desired Bayesian solution to the inverse
problem, given in (118). A variant on this idea is to morph from prior to
posterior by (possibly artificially) considering the data $w$ as sequentially
acquired and incrementally including components of the data at each step $n$.

These sequential formulations of Bayesian inversion are well known and have,
for example, been exploited in the use of sequential Monte Carlo methods for
Bayesian inference; see Subsection 4.6 for details. However these sequential
Monte Carlo methods do not always scale well to high dimensions, as for state
estimation. Consequently ensemble Kalman variants of them have an important
place in the field. In the following subsections we use (123) as the basis for
ensemble Kalman filtering frameworks in discrete time, and subsequently (in
the limit $N\to\infty$) in continuous time.

####  4.2.1 Bayesian Inversion Using Discrete Time Filtering

Consider the combined state-observation system in the form, for
$n\in\\{0,\cdots,N-1\\}$,

(124a) $\displaystyle u_{n+1}$ $\displaystyle=u_{n},$ (124b) $\displaystyle
w_{n+1}$ $\displaystyle=G(u_{n+1})+\frac{1}{\sqrt{\Delta t}}\gamma_{n+1},$

where $\\{\gamma_{n+1}\\}_{n=0}^{N-1}$ is an i.i.d. sequence with variance
$\mathsf{N}(0,\mathsf{\Gamma}).$ It is intuitive that, if $N\Delta t=1$,
$u_{0}\sim\mathsf{N}(m_{0},C_{0})$ and $w_{n+1}=w$ for all
$n\in\\{0,\cdots,N-1\\}$, then $u_{N}$ conditioned on
$\\{w_{n+1}\\}_{n=0}^{N-1}$ is distributed as $\mu^{w}.$ This follows since
use of (124b) for $n\in\\{0,\cdots,N-1\\}$ corresponds to making $N$
independent noisy observations of $G(u_{0})$, all with noise variance $\Delta
t^{-1}\mathsf{\Gamma};$ this is statistically equivalent to a single noisy
observation of $G(u_{0})$ with noise variance $\mathsf{\Gamma}.$ Since $u_{0}$
is initialized as $\mathsf{N}(m_{0},C_{0})$ (the prior) the problem reduces to
the Bayesian inverse problem for $u|w.$

This intuition may be substantiated using the discussion above, following
equation (123). We prove the desired result but, before doing so, and in order
to avail ourselves of the results from Sections 2 and 3, we rescale the
observation equation in (124) to obtain

(125a) $\displaystyle u_{n+1}$ $\displaystyle=u_{n},$ (125b) $\displaystyle
y_{n+1}$ $\displaystyle=\Delta tG(u_{n+1})+\eta_{n+1},$

with $\eta_{n+1}\sim\mathsf{N}(0,\Delta t\mathsf{\Gamma})$ and $y_{n+1}=\Delta
tw_{n+1}.$ Let ${Y}^{\dagger}_{n}=\\{y_{\ell}\\}_{\ell=1}^{n}.$ We may then
show the following:

######  Theorem 4.1.

Assume that $C_{0}\succ 0$, $\mathsf{\Gamma}\succ 0$ and $N\Delta t=1$. Assume
also that $u_{0}\sim\mathsf{N}(m_{0},C_{0})$ and
$\eta_{n+1}\sim\mathsf{N}(0,\Delta t\mathsf{\Gamma});$ furthermore, assume
that $\\{\eta_{n}\\}_{n=1}^{N}$ forms an i.i.d. sequence, independent of
$u_{0}.$ Then $\mu_{N}$, the law of $u_{N}|{Y}^{\dagger}_{N}$, is equal to the
posterior distribution $\mu^{w}$ if the data is chosen as $y_{n}=\Delta t\,w,$
for $n\in\\{1,\cdots,N\\}.$

######  Proof 4.2.

Let $\mu_{n}$ be the law of $u_{n}|{Y}^{\dagger}_{n}$. Since (in the notation
of Section 2) $\widehat{\mu}_{n+1}=\mu_{n}$ we see that the mapping $\mu_{n}$
to $\mu_{n+1}$ is simply given by Bayes theorem: $\mu_{n+1}=L_{n}(\mu_{n}).$
This observation yields the identity, expressed in terms of $\rho_{n}$ the
Lebesgue density of measure $\mu_{n}:$

$\displaystyle\log\rho_{n+1}-\log\rho_{n}$ $\displaystyle=-\frac{1}{2\Delta
t}|y_{n}-\Delta tG(u)|_{\mathsf{\Gamma}}^{2}+{\rm const},$
$\displaystyle=-\frac{\Delta t}{2}|w-G(u)|_{\mathsf{\Gamma}}^{2}+{\rm const}.$

Summing over $n\in\\{0,\cdots,N-1\\}$, using the fact that

$\displaystyle\log\rho_{0}$
$\displaystyle=-\frac{1}{2}|u-m_{0}|_{C_{0}}^{2}+{\rm const},$

and exponentiating we obtain the desired result. We also deduce that

$\log\rho_{n}=-\Phi_{R,n}+{\rm const},$

with $\Phi_{R,n}$ given by (123).

Thus we may approach the problem of sampling from $\mu^{w}$ by (approximately)
solving the filtering problem for $\mu_{n}$, until discrete time $n=N$, using
ensemble Kalman methods. To illustrate this idea we describe the Gaussian
projected filter in this setting, and then the mean field model of Kalman
transport form, both applied to the specific state-observation model (125);
particle approximations of the mean field models lead to Kalman based
algorithms for the inverse problem. We discuss both the discrete and
continuous time settings and, to this end, we first study the continuous limit
of (125) found by sending $\Delta t\to 0$, with $N\Delta t=1$ fixed.

####  4.2.2 Bayesian Inversion Using Continuous Time Filtering

The reparametrization (73a) yields the continuum limit of (125) in the form of
the SDE

(126a) $\displaystyle du$ $\displaystyle=0,$ (126b) $\displaystyle dz$
$\displaystyle=G(u)dt+\sqrt{\mathsf{\Gamma}}dB.$

with $B$ a unit Brownian motion in $\mathbb{R}^{d_{z}}.$ We define
$Z^{\dagger}(t)=\\{{z}^{\dagger}(s)\\}_{0\leq s\leq t},$ and consider the
filtering distribution for the random variable $u(t)|Z^{\dagger}(t).$ However
there is a twist on the standard filtering setting: we are interested in the
case where the data has constant derivative: $d{z}^{\dagger}(t)/dt=w.$ Since
the path ${z}^{\dagger}$ has zero quadratic variation the probability
distribution is found by setting ${z}^{\dagger}(t)=tw$ within the Stratonovich
formulation of the non-local evolution equation for the density. Referring to
(87) we see that this yields the following evolution for density $\rho(u,t)$
of $u(t)|Z^{\dagger}(t):$

(127) $\displaystyle\partial_{t}\rho=\left\langle
G-\mathbb{E}G,w\right\rangle_{\mathsf{\Gamma}}\rho-\frac{1}{2}\left\\{\left|G\right|^{2}_{\mathsf{\Gamma}}-\mathbb{E}\left|G\right|^{2}_{\mathsf{\Gamma}}\right\\}\rho.$

Here $\mathbb{E}$ denotes integration with respect to density $\rho(\cdot,t)$
so that the equation is non-local with respect to variable $u$. This is the
analogue of the Kushner-Stratonovich equation for the filtering problem
defined by (126), since the unconditioned variable $u$ has trivial dynamics
and since we are studying the case where the data has zero quadratic variation
and is in fact differentiable.

We may then show the following:

######  Theorem 4.3.

Assume that $C_{0}\succ 0$, $\mathsf{\Gamma}\succ 0$ and
$u_{0}\sim\mathsf{N}(m_{0},C_{0})$. Let $\mu(t)$ be the law of
$u(t)|Z^{\dagger}(t)$, with data chosen as $z^{\dagger}(t)=tw,$ for
$t\in(0,1).$ Then $\mu(t)$ has density $\rho(\cdot,t)$ satisfying (127), with
solution given by the formulae

(128a) $\displaystyle\rho(u,t)$ $\displaystyle=\frac{1}{{\cal
Z}(t)}\exp\bigl{(}-t\Phi(u)\bigr{)}\rho_{0}(u),$ (128b) $\displaystyle{\cal
Z}(t)$
$\displaystyle=\int_{\mathbb{R}^{d_{u}}}\exp\bigl{(}-t\Phi(u)\bigr{)}\rho_{0}(u)du;$

in particular $\mu(1)$ is equal to the posterior distribution $\mu^{w}$.

######  Proof 4.4.

Let $\rho_{0}$ denote the probability density function of the prior
$\mathsf{N}(m_{0},C_{0}).$ Recall that $\mu(t)$ has density given by equation
(127) with initial condition $\rho|_{t=0}=\rho_{0}.$ First note that (128b)
gives

$\frac{d{\cal Z}}{dt}=-\mathbb{E}\Phi\,{\cal Z}$

and hence that (128a) gives

$\displaystyle\partial_{t}\rho$ $\displaystyle=-(\Phi-\mathbb{E}\Phi)\rho,$
$\displaystyle=\bigl{(}\langle
w,G\rangle_{\mathsf{\Gamma}}-\frac{1}{2}|G|_{\mathsf{\Gamma}}^{2}\bigr{)}\rho-\mathbb{E}\bigl{(}\langle
w,G\rangle_{\mathsf{\Gamma}}-\frac{1}{2}|G|_{\mathsf{\Gamma}}^{2}\bigr{)}\rho,$
$\displaystyle=\left\langle
G-\mathbb{E}G,w\right\rangle_{\mathsf{\Gamma}}\rho-\frac{1}{2}\left\\{\left|G\right|^{2}_{\mathsf{\Gamma}}-\mathbb{E}\left|G\right|^{2}_{\mathsf{\Gamma}}\right\\}\rho.$

Thus $\rho$ given by (128) solves (127), and since (127) characterizes the law
of $u(t)|Z^{\dagger}(t)$ when $\frac{dz^{\dagger}}{dt}=w$, the result is
proven.

####  4.2.3 Gaussian Projected Filter: Discrete Time

We now apply the ideas from Subsection 2.4, on the Gaussian projected filter
in the general setting, to the specific setting of (4.2). Let $\mathbb{E}$
denote expectation under $u\sim\mathsf{N}(m_{n},C_{n})$ and define

(129a) $\displaystyle C_{n}^{uG}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}u-\mathbb{E}u\bigr{)}\otimes\bigl{(}G(u)-\mathbb{E}G(u)\bigr{)}\Bigr{)},$
(129b) $\displaystyle C_{n}^{GG}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}G(u)-\mathbb{E}G(u)\bigr{)}\otimes\bigl{(}G(u)-\mathbb{E}G(u)\bigr{)}\Bigr{)}.$

Noting that prediction under (4.2a) is trivial it follows that the predicted
mean and covariance satisfy $\widehat{m}_{n+1}=m_{n}$ and
$\widehat{C}_{n+1}=C_{n}.$ Hence, using (31) in the specific setting of (4.2),
yields

(130a) $\displaystyle m_{n+1}$ $\displaystyle=m_{n}+\Delta
tC_{n}^{uG}(\mathsf{\Gamma}+\Delta
tC_{n}^{GG})^{-1}\bigl{(}w-\mathbb{E}G(u)\bigr{)},$ (130b) $\displaystyle
C_{n+1}$ $\displaystyle=C_{n}-\Delta tC_{n}^{uG}(\mathsf{\Gamma}+\Delta
tC_{n}^{GG})^{-1}\bigl{(}C_{n}^{uG}\bigr{)}^{\top}.$

Note that the difference between the data $w$ and the mean of $G(u)$ under the
Gaussian at step $n$ acts as a forcing term in the evolution of the mean from
$n$ to $n+1$, promoting a Gaussian which agrees with the data. This forcing
term is weighted by covariance information. The covariance of the Gaussian
projected filter decreases from step to step since $\langle
u,C_{n+1}u\rangle\leq\langle u,C_{n}u\rangle$ for all
$u\in\mathbb{R}^{d_{u}};$ this reflects the fact that more information is
received at each step $n\mapsto n+1$ as the unknown $u$ is repeatedly
observed.

In the setting of linear inverse problems, where $G(u)=Lu$, the equations
reduce to

(131a) $\displaystyle m_{n+1}$ $\displaystyle=m_{n}+\Delta
tC_{n}L^{\top}(\mathsf{\Gamma}+\Delta
tLC_{n}L^{\top})^{-1}\bigl{(}w-Lm_{n}\bigr{)},$ (131b) $\displaystyle C_{n+1}$
$\displaystyle=C_{n}-\Delta tC_{n}L^{\top}(\mathsf{\Gamma}+\Delta
tLC_{n}L^{\top})^{-1}LC_{n}.$

We know from Section 2 that this is simply the Kalman filter applied to our
specific linear and Gaussian dynamical system, and that it exactly solves the
filtering problem:

######  Lemma 4.5.

Assume that $C_{0}\succ 0$, $\mathsf{\Gamma}\succ 0$ and that $G(u)=Lu$ for
$L\in\mathbb{R}^{d_{w}\times d_{u}}$. Then the posterior measure $\mu^{w}$ is
Gaussian with mean $m_{N}$ and covariance matrix $C_{N}$ given by (131).

######  Proof 4.6.

Theorem 4.1 shows that, if $w_{n}\equiv w$ for $n\in\\{1,\cdots,N\\}$, then
the filtering distribution for (124) at time $N$ equals the posterior
distribution $\mu^{w}$. Since $G$ is linear this posterior distribution is
Gaussian and the filtering distribution for (124) is Gaussian. In Section 2 we
show that the Gaussian projected filter gives the exact filtering distribution
in the linear Gaussian setting, and hence the desired result follows.

####  4.2.4 Gaussian Projected Filter: Continuous Time

To further elucidate the structure of Gaussian projected filtering for the
inverse problem, we study its continuous time formulation from Subsection 3.4
when applied to the specific state-observation model (125), and its continuous
time limit (126). Taking the limit $\Delta t\to 0$ in (130) we obtain the
following evolution equations for mean and covariance:

(132a) $\displaystyle\frac{dm}{dt}$
$\displaystyle=C^{uG}\mathsf{\Gamma}^{-1}\bigl{(}w-\mathbb{E}G(u)\bigr{)}$
(132b) $\displaystyle\frac{dC}{dt}$
$\displaystyle=-C^{uG}\mathsf{\Gamma}^{-1}(C^{uG})^{\top},$ (132c)
$\displaystyle C^{uG}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}u-\mathbb{E}u\bigr{)}\otimes\bigl{(}G(u)-\mathbb{E}G(u)\bigr{)}\Bigr{)}$

and in which all where expectations are computed under
$\mathsf{N}\bigl{(}m(t),C(t)\bigr{)}.$

We complete this subsection by solving equations (132) explicitly in the
setting of linear $G$. Consider $\rho$ solving the Kushner-Stratonovich
equation (127), with solution given by (128). In the general linear setting we
know that the Gaussian projected filter is exact. In the current inverse
problem setting this leads to Lemma 4.7 below. For this we note that, in the
linear setting $G(u)=Lu$, the Gaussian projected filter (132) becomes

(133a) $\displaystyle\frac{dm}{dt}$
$\displaystyle=CL^{\top}\mathsf{\Gamma}^{-1}(w-Lm)$ (133b)
$\displaystyle\frac{dC}{dt}$ $\displaystyle=-CL^{\top}\mathsf{\Gamma}^{-1}LC,$
(133c) $\displaystyle C$
$\displaystyle=\mathbb{E}\bigl{(}(u-\mathbb{E}u)\otimes(u-\mathbb{E}u)\bigr{)},$

and again all expectations are computed under
$\mathsf{N}\bigl{(}m(t),C(t)\bigr{)}.$

######  Lemma 4.7.

Assume that $C_{0}\succ 0$, $\mathsf{\Gamma}\succ 0$ and that $G(u)=Lu$ for
$L\in\mathbb{R}^{d_{w}\times d_{u}}$. Then $\rho$ solving equation (127) is
given by the Gaussian $\mathsf{N}\bigl{(}m(t),C(t)\bigr{)}$ where $m(t),C(t)$
solve the Gaussian projected filter equations (133). In particular the
posterior measure $\mu^{w}$ is Gaussian and given by
$\mathsf{N}\bigl{(}m(1),C(1)\bigr{)}.$

######  Proof 4.8.

From Theorem 4.3, and in particular equation (128), we see that, in the linear
case $G(\cdot)=L\cdot$, the solution of (127) is given by

(134)
$\rho(u,t)\propto\exp\Bigl{(}-\frac{t}{2}|w-Lu|_{\mathsf{\Gamma}}^{2}-\frac{1}{2}|u-m_{0}|_{C_{0}}^{2}\Bigr{)}.$

Completing the square shows that this density corresponds to Gaussian
$\mathsf{N}\bigl{(}m(t),C(t)\bigr{)}$ with mean and covariance satisfying

(135a) $\displaystyle C(t)^{-1}m(t)$
$\displaystyle=tL^{\top}\mathsf{\Gamma}^{-1}w+C_{0}^{-1}m_{0},$ (135b)
$\displaystyle C(t)^{-1}$
$\displaystyle=C_{0}^{-1}+tL^{\top}\mathsf{\Gamma}^{-1}L.$

Note that since $C_{0}\succ 0$ it follows that $C_{0}^{-1}\succ 0$ and hence
(135b) shows that, for all $t\geq 0$, $C(t)^{-1}\succ 0$ and hence that
$C(t)\succ 0$ for all $t\geq 0$; hence $C(t)$ is well-defined by (135b) and
$m(t)$ is well-defined by (135a). It simply remains to show that $m(t)$ and
$C(t)$ given by these formulae solve (133) when $m(0)=m_{0}$ and $C(0)=C_{0}.$

We thus turn our attention to the equations (133). Note that $C(t)$ solving
(133b) satisfies $C(0)=C_{0}\succ 0.$ Hence $C(0)^{-1}\succ 0.$ Thus, by
continuity, $C(t)$ remains invertible for some positive interval of time
$t\in[0,\tau)$ and, on this interval, direct computation with (133b) shows
that

(136) $\frac{dC^{-1}}{dt}=L^{\top}\Gamma^{-1}L.$

From this it follows by integration that $C^{-1}(t)\succeq C_{0}^{-1}\succ 0$
for all $t$ and hence that we may take $\tau=\infty.$ Furthermore, the
integration also shows that the solution of (136), solving (133b), delivers
(135b) as desired.

We then notice that, from (133a),

$\displaystyle C^{-1}\frac{dm}{dt}$
$\displaystyle=L^{\top}\mathsf{\Gamma}^{-1}w-L^{\top}\mathsf{\Gamma}^{-1}Lm,$
$\displaystyle=L^{T}\mathsf{\Gamma}^{-1}w-\frac{dC^{-1}}{dt}m.$

It follows that

$\frac{d}{dt}\bigl{(}C^{-1}m)=L^{\top}\mathsf{\Gamma}^{-1}w$

and integration, together with use of the initial conditions, shows that
(133a) delivers the desired identity (135a).

######  Remark 4.9.

It is also possible to establish the preceding result by studying the discrete
time problem and taking a limit. Recall that the solution $\rho$ of (127) is
given by (134) in the linear setting. By Theorem 4.3 we know that
$\rho(\cdot,1)$ is the density of the posterior distribution. Furthermore, for
any $t>0$, $\rho(u,t)$ is the posterior distribution for the Bayesian
inference problem with scaled likelihood function
$\exp(-\frac{t}{2}|w-Lu|_{\mathsf{\Gamma}}^{2})$. The associated posterior
mean $m(t)$ and covariance matrix $C(t)$ satisfy (135); application of the
Woodbury matrix identity shows that

$\displaystyle m(t)$
$\displaystyle=m_{0}+C_{0}L^{\top}(LC_{0}L^{\top}+t^{-1}\mathsf{\Gamma})^{-1}(w-Lm_{0}),$
$\displaystyle C(t)$
$\displaystyle=C_{0}-C_{0}L^{\top}(LC_{0}L^{\top}+t^{-1}\Gamma)^{-1}LC_{0};$

these expressions also follow from the Kalman filter update formulae (34) with
$\Gamma$ being replaced by $\mathsf{\Gamma}/t$, and identifications of the
predictive distribution in (34) with the prior here, and observation operator
$H$ in (34) with $L$ here. Furthermore, we may take
$\mathsf{N}\bigl{(}m(t),C(t)\bigr{)}$ as prior and apply the scaled likelihood
$\exp(-\frac{\Delta t}{2}|w-Ly|_{\mathsf{\Gamma}}^{2})$ to obtain posterior
$\mathsf{N}\bigl{(}m(t+\Delta t),C(t+\Delta t)\bigr{)}$ with associated Kalman
update formula

$\displaystyle m(t+\Delta t)$
$\displaystyle=m(t)+C(t)L^{\top}(LC(t)L^{\top}+\Delta
t^{-1}\mathsf{\Gamma})^{-1}(w-Lm(t)),$ $\displaystyle C(t+\Delta t)$
$\displaystyle=C(t)-C(t)L^{\top}(LC(t)L^{\top}+\Delta
t^{-1}\Gamma)^{-1}LC(t),$

again using (34). It simply remains to show that taking the limit $\Delta t\to
0$ yields evolution equations (133), which indeed follow trivially from
observing

$C(t)L^{\top}(LC(t)L^{\top}+\Delta t^{-1}\mathsf{\Gamma})^{-1}=\Delta
tC(t)L^{\top}\mathsf{\Gamma}^{-1}+\mathcal{O}(\Delta t^{2}).\quad\blacksquare$

####  4.2.5 Kalman Transport: Discrete Time

Recall that we introduce Gaussian projected filtering in Section 2 as a
stepping stone to mean field dynamical systems. These mean field models in
turn lead to ensemble Kalman methods through particle approximation. In this
Section 4, devoted to inverse problems, we simply highlight use of one of the
mean field models, namely the stochastic _Kalman transport_ model. We leave
details of particle approximations of it to the reader, to the brief
discussions of the topic in Subsection 4.5 and to Appendix 6 for related
pseudo code. The adaptation of deterministic Kalman transport filter models to
inverse problems is briefly discussed in Subsection 4.2.7 below.

Employing the state-observation model (125) within (46) we obtain the mean-
field dynamical system, for i.i.d. unit Gaussian sequence $\\{\xi_{n}\\}$ in
$\mathbb{R}^{d_{w}}$,

(137a) $\displaystyle u_{n+1}$ $\displaystyle=u_{n}+\Delta tC^{uG}_{n}(\Delta
tC^{GG}_{n}+\mathsf{\Gamma})^{-1}\Bigl{(}w-G(u_{n})-\sqrt{\frac{\mathsf{\Gamma}}{\Delta
t}}\xi_{n}\Bigr{)},$ (137b) $\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{)},$
(137c) $\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{)}.$

Note that, here, expectation $\mathbb{E}$ is computed under the law of $u_{n}$
itself. The ideas of Section 2 show that random variable $u_{N}$ provides an
approximation to the posterior distribution $\mu^{w}$, provided that
$u_{0}\sim\mathsf{N}(m_{0},C_{0})$, the prior distribution. More can be said
in the linear case:

######  Lemma 4.10.

Assume that $C_{0}\succ 0$ and $\mathsf{\Gamma}\succ 0$ and that $G(u)=Lu$ for
$L\in\mathbb{R}^{d_{w}\times d_{u}}$. Then the solution of the mean field
model (137) satisfies $u_{N}\sim\mu^{w}.$

######  Proof 4.11.

In the linear setting equation (137) defines a closed evolution in the set of
Gaussian probability measures. Furthermore, some calculation confirms that the
mean and covariance of $u_{N}$ then coincide with the Gaussian projected
filter, in this linear case, given by equations (131). Indeed, by taking the
expectation under the law of $u_{n}$ of (137a), it is readily checked that in
the linear setting

(138) $m_{n+1}=m_{n}+\Delta tC_{n}L^{\top}(\mathsf{\Gamma}+\Delta
tLC_{n}L^{\top})^{-1}\bigl{(}w-Lm_{n}\bigr{)}.$

To obtain the evolution equation of the covariance, we first recall that

(139)
$C_{n+1}=\mathbb{E}\Bigl{(}(u_{n+1}-m_{n+1})\otimes(u_{n+1}-m_{n+1})\Bigr{)}.$

Substituting into (139) the expression for $u_{n+1}$, given by (137a) in the
linear setting, and the expression for $m_{n+1}$, given by (138), and then
computing the expectation yields

(140) $C_{n+1}=C_{n}-\Delta tC_{n}L^{\top}(\mathsf{\Gamma}+\Delta
tLC_{n}L^{\top})^{-1}LC_{n}.$

Thus equations (138), (140) are identical to (131) and it follows from Lemma
4.5 that $u_{N}\sim\mu^{w}.$

####  4.2.6 Kalman Transport: Continuous Time

We now study the inverse problem using Kalman transport in continuous time.
Upon taking the continuous time limit of (137) we formally obtain

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

Here $B\in\mathbb{R}^{d_{w}}$ is a unit Brownian motion and expectation is
under the law of $u$ itself. The SDE (141) may also be seen as a consequence
of (106) applied to the specific case of state-observation model (125).
However, special care is required since the observations $z^{\dagger}$ in
Section 3 were assumed to have non vanishing quadratic variation while
$dz^{\dagger}/dt=w$ in this section.

To obtain further insight into this mean-field dynamical system we once again
consider the linear setting:

######  Lemma 4.12.

Assume that $C_{0}\succ 0$ and $\mathsf{\Gamma}\succ 0$. Consider the SDE
(141) in the setting where $G(u)=Lu$ for matrix $L\in\mathbb{R}^{d_{w}\times
d_{u}}$. Then $u(1)\sim\mu^{w}.$

######  Proof 4.13.

The Gaussian projected filter is exact in the linear setting, by Lemma 4.7,
and hence delivers the desired posterior at time $t=1$, by Theorem 4.3. Thus
it suffices to show that the mean and covariance of $u$ from (141) satisfy the
Gaussian projected filter in the linear setting; this is given by (133).

We first note that

$du=CL^{\top}\mathsf{\Gamma}^{-1}\bigl{(}w\,dt-
Lu\,dt-\sqrt{\mathsf{\Gamma}}dB\bigr{)}$

where $C$ is the covariance of $u$ under its law. By the Itô formula,
$m=\mathbb{E}u$ satisfies (133a). It follows that $e=u-m$ satisfies

$de=-CL^{\top}\mathsf{\Gamma}^{-1}Le\,dt-CL^{\top}\mathsf{\Gamma}^{-1/2}dB.$

A second use of the Itô formula shows that $C=\mathbb{E}(e\otimes e)$
satisfies (133b). The proof is complete.

####  4.2.7 Non-Stochastic Transports

We briefly summarize corresponding deterministic transport formulations, both
in discrete and continuous time. Following the approximation (2.10), which
holds in our case provided $\Delta t$ is small enough, and choosing $K_{n}$ as
implicitly defined in (137), we obtain

$u_{n+1}=u_{n}+\Delta tC^{uG}_{n}(\Delta
tC^{GG}_{n}+\mathsf{\Gamma})^{-1}\Bigl{(}w-\frac{1}{2}\bigl{(}G(u_{n})+\mathbb{E}G(u_{n})\bigr{)}\Bigr{)}.$

Furthermore, taking the $\Delta t\to 0$ limit results in the mean field ODE
formulation

(142)
$\frac{du}{dt}=C^{uG}\mathsf{\Gamma}^{-1}\Bigl{(}w-\frac{1}{2}\bigl{(}G(u)+\mathbb{E}G(u)\bigr{)}\Bigr{)}.$

Analogously to Lemma 4.12 we have

######  Lemma 4.14.

Assume that $C_{0}\succ 0$ and $\mathsf{\Gamma}\succ 0$. Consider the mean
field model (142) in the setting where $G(u)=Lu$ for matrix
$L\in\mathbb{R}^{d_{w}\times d_{u}}$. Then $u(1)\sim\mu^{w}.$

######  Proof 4.15.

The Gaussian projected filter is exact in the linear setting, by Lemma 4.7,
and hence delivers the desired posterior at time $t=1$, by Theorem 4.3. Thus
it suffices to show that the mean and covariance of $u$ from (142) satisfy the
Gaussian projected filter in the linear setting; this is given by (133).

We first note that the mean under (142) satisfies

$\frac{dm}{dt}=CL^{T}\mathsf{\Gamma}^{-1}(w-Lm),$

which is (133a). Using this it also follows that

$\frac{d}{dt}(u-m)=-\frac{1}{2}CL^{T}\mathsf{\Gamma}^{-1}L(u-m)$

from which it follows that the variance satisfies (133b).

###  4.3 Geometric Structure in Kalman Transport, Optimization and Sampling

So far we have taken a Bayesian perspective on solving inverse problems,
considering both discrete and continuous time transport approaches in
Subsection 4.2. The methods were characterized by the requirement that either
$n=1,\ldots,N$ subject to $N\Delta t=1$ or $t\in[0,1]$, and must be
initialized at the prior. In this subsection, we lift these restrictions and
consider the behavior of continuous time Kalman transport methods in the
asymptotic limit $t\to\infty$, leading to algorithms which can be initialized
freely. Analogous investigations could be conducted for the discrete time
formulations with $\Delta t$ fixed and iteration index $N\to\infty$ and are
left to the reader; furthermore the following Subsection 4.4 also covers
discrete time methods based on iterating to infinity.

We show that Kalman transport methods, when viewed from such an asymptotic
perspective, are closely related to gradient descent methods for optimization.
However, the resulting methodologies lead to new perspectives: firstly, Kalman
transport methods possess the property of _affine invariance_ leading, for
linear inverse problems, to methods with convergence rate that is independent
of the conditioning of the problem; second, they constitute mean field
equations and give rise to a novel gradient flow structure in the space of
probability measures; third, the mean field equations are gradient free while
still being closely related to gradient based optimization methods; and fourth
the ideas can be generalized from optimization to sampling.

In Subsection 4.3.1 we study gradient descent, making connections to affine
invariance via a mean field preconditioning. In Subsection 4.3.2 we show how
_statistical linearization_ can be used to define ensemble Kalman variants on
numerous iterative methods for least squares problems. And in Subsection 4.3.3
we combine the two preceding subsections by introducing a derivative-free
approximation of mean field preconditioned gradient flow, obtained by
statistical linearization. In Subsection 4.3.4 we generalize the foregoing
considerations from gradient descent to the Langevin equation, via development
of Kalman-Wasserstein gradient flow in the space of probability measures.

####  4.3.1 Gradient Flow Structures In Parameter Space

We start by recalling some properties of standard gradient descent when
applied to an objective function $\Psi:\mathbb{R}^{d}\to\mathbb{R}^{+}$ and
given symmetric positive-definite preconditioner $B\in\mathbb{R}^{d_{u}\times
d_{u}}:$

(143) $\frac{du}{dt}=-B\nabla\Psi(u).$

Note that, along solutions of (143),

$\displaystyle\frac{d}{dt}\Psi(u)$
$\displaystyle=\Big{\langle}\nabla\Psi(u),\frac{du}{dt}\Big{\rangle}$
$\displaystyle=-\Big{\langle}B^{-1}\frac{du}{dt},\frac{du}{dt}\Big{\rangle}$
$\displaystyle=-|B^{\frac{1}{2}}\nabla\Psi(u)|^{2}\leq 0.$

In particular we highlight that

$\frac{d}{dt}\Psi(u)=-\Bigl{|}\frac{du}{dt}\Bigr{|}_{B}^{2}.$

Equation (143) is said to possess a gradient flow structure in parameter space
$\mathbb{R}^{d_{u}}$ because the vector field driving the evolution of $u$ is
tangential to the gradient of the objective function in a Euclidean metric
weighted by $B$; this is behind the non-increasing property of $\Psi(u)$ along
trajectories.

Consider now the affine transformation $u\mapsto\tilde{u}$ given by

(144) $\tilde{u}=Au+b,$

where $A$ is an invertible matrix and $b$ a vector. An important issue in all
vector space optimization problems is the relative scaling of the components
of the vector. A highly desirable feature of an algorithm is that it be
insensitive to such scaling issues. This can be addressed by looking at
differences between: (i) the algorithm for $u$ rewritten in terms of
$\tilde{u}$ given by (144); (ii) the same algorithm applied directly to
variable $\tilde{u}$ optimizing $\tilde{\Psi}(\tilde{u})$, with the latter
defined by

(145) $\tilde{\Psi}(\tilde{u})=\Psi\Bigl{(}A^{-1}(\tilde{u}-b)\Bigr{)}.$

When these two ways of using reparametrization (144) lead to the same
algorithm, for all choices of $A,b$, we say that the algorithm is _affine
invariant_. Such algorithms are highly desirable as they are not sensitive to
the scaling of the variables which, at the optimum, may be unknown _a priori_.

Applying the transformation (144) to (143) leads to

(146)
$\frac{d\tilde{u}}{dt}=-AB\nabla\Psi\Bigl{(}A^{-1}(\tilde{u}-b)\Bigr{)},$

the descent approach underlying algorithm viewed as in (i). In contrast,
applying the same gradient descent to $\tilde{\Psi}$ given by (145), leads to
the the descent approach underlying algorithms viewed as in (ii):

(147)
$\frac{d\tilde{u}}{dt}=-BA^{-\top}\nabla\Psi\Bigl{(}A^{-1}(\tilde{u}-b)\Bigr{)}.$

The two equations (146), (147) only agree if

$AB=BA^{-\top}$

and such an identity cannot hold for all $A$, for a fixed $B$. Thus the basic
gradient descent (143) is not affine invariant. However it is a remarkable
fact that, by generalizing (143) to allow for mean-field dependence, we can
achieve affine invariance.

To this end consider the mean-field generalization of (143)

(148) $\frac{du}{dt}=-B(\rho)\nabla\Psi(u),$

where $\rho(\cdot,t)$ is the probability density function associated with the
law of $u$, assuming that $u$ is initialized from probability density function
$\rho_{0}(\cdot).$ If we assume that $B(\cdot)$ is positive-definite symmetric
for all possible input densities, then similar arguments to before show that
$\Psi(u)$ is non-increasing along trajectories of (148). Furthermore similar
arguments to before show that the algorithm is affine invariant if, for all
invertible $A\in\mathbb{R}^{d_{u}\times d_{u}}$,

$AB(\rho)=B(\tilde{\rho})A^{-\top},$

where $\tilde{\rho}$ is the density of $\tilde{u}$ related to $u$, with
density $\rho$, by (144). This identity holds if $B(\cdot)$ is chosen to be
the covariance associated with its argument. We have thus discovered the
affine invariant mean field gradient descent

(149) $\frac{du}{dt}=-C\nabla\Psi(u),$

where $C$ is the covariance operator under the law of $u(t)$ and $u(0)$ is
chosen at random from probability measure with density $\rho_{0}.$

Algorithms based on solving (149) are not themselves ensemble Kalman methods,
although we have drawn inspiration from the power of mean field methods to
motivate the approach. In the next subsection we introduce the idea of
statistical linearization, leading to a variety of ensemble Kalman methods for
optimization; and in the subsection following it we use this idea to
approximate (149) by an ensemble Kalman version of gradient descent which
obviates the need for computing adjoints of the forward model.

####  4.3.2 Statistical Linearization

A basic building block in the Gaussian projected filter and mean-field models
for inverse problems that we have presented in Subsection 4.2 is the object

(150)
$C^{uG}=\mathbb{E}\Bigl{(}\bigl{(}u-\mathbb{E}u\bigr{)}\otimes\bigl{(}G(u)-\mathbb{E}G(u)\bigr{)}\Bigr{)},$

evolving in both discrete and continuous time, and its regularized analog
arising when $G_{R}$, as defined in (119), is used in place of $G$. In a
methodology based on exact properties only for first and second order
statistics, it is natural that $C^{uG}$ should appear when solving the inverse
problem: the correlation between the parameter $u$, which we wish to estimate,
and $G(u)$ to which we have access when evaluated at a noise-polluted version
of the true value, is a fundamental object. One way of understanding its role
in algorithms for inversion is through the idea of statistical linearization,
providing a link between ensemble methods and derivatives of the objective
function. The underlying principle is that the differences used in ensemble
methods, and covariances in particular, act as a surrogate for derivatives.

The expectation defining (150) is computed under the distribution of a
Gaussian (for the Gaussian projected filter) or a more general distribution
(mean-field model). To get some insight into the connection between ensemble
methods and derivatives, we first consider the setting where
$u\sim\mathsf{N}\bigl{(}m,C\bigr{)}.$ Note that such $u$ can be written as
$u=m+\sqrt{C}\xi$ where $\xi\sim\mathsf{N}(0,I).$ Thus (150) may be
reformulated as

(151)
$C^{uG}=\mathbb{E}\Bigl{(}\bigl{(}\sqrt{C}\xi\bigr{)}\otimes\bigl{(}G(m+\sqrt{C}\xi)-\mathbb{E}G(m+\sqrt{C}\xi)\bigr{)}\Bigr{)},\quad\xi\sim\mathsf{N}(0,I).$

Using this we obtain the following connection between $C^{uG}$ and derivatives
of $G:$

######  Lemma 4.16.

Assume that the second derivative of $G$ is small: there is $\epsilon\ll 1$
such that

$\sup_{u\in\mathbb{R}^{d_{u}}}|D^{2}G(u)[\zeta,\zeta]|\leq\epsilon|\zeta|^{2}.$

Then $C^{uG}$ given by (151) satisfies

$C^{uG}=CDG(m)^{\top}+{\mathcal{O}}(\epsilon).$

Thus, when $C^{-1}\succ\lambda I$, for some $\lambda>0$ independent of
$\epsilon$,

(152) $DG(m)=(C^{uG})^{\top}C^{-1}+{\mathcal{O}}(\epsilon).$

######  Proof 4.17.

Note that

$\displaystyle G(m+\sqrt{C}\xi)$
$\displaystyle=G(m)+DG(m)\sqrt{C}\xi+{\mathcal{O}}(\epsilon),$
$\displaystyle\mathbb{E}G(m+\sqrt{C}\xi)$
$\displaystyle=G(m)+{\mathcal{O}}(\epsilon).$

From (151)

$C^{uG}=\mathbb{E}\Bigl{(}\bigl{(}\sqrt{C}\xi\bigr{)}\otimes\bigl{(}DG(m)\sqrt{C}\xi+{\mathcal{O}}(\epsilon)\bigr{)}\Bigr{)},\quad\xi\sim\mathsf{N}(0,I),$

and the desired result follows.

######  Remark 4.18.

Another perspective on the preceding lemma is via Stein’s identity. This
states that

$C^{uG}=C(\mathbb{E}DG)^{\top}$

for (151) when expectation is computed under a Gaussian measure
$\mathsf{N}(m,C)$. The identity can be verified via integration by parts.
Given the assumption stated in Lemma 4.16,
$\mathbb{E}DG(u)=DG(m)+\mathcal{O}(\epsilon)$ and the approximation result
(152) also follows. $\blacksquare$

When $D^{2}G$ is indeed small it is not unreasonable to use (152) as the basis
for an approximation to $DG(\cdot)$ at points in an ${\mathcal{O}}(1)$ ball
around the mean. We now take this idea further and consider random variable
$u\in\mathbb{R}^{d_{u}}$ (not necessarily Gaussian) and compute $C$ and
$C^{uG}$ as covariance of $u$ and cross-covariance of $u$ with $G(u)$
respectively. We refer to use of the approximation

(153) $DG(u)\approx(C^{uG})^{\top}C^{-1},\quad\text{for
all}\,u\in\mathbb{R}^{d_{u}},$

as _statistical linearization._ The approximation can be invoked to replace
the derivative within any standard optimization or sampling algorithm to solve
the inverse problem (117). Doing so results in a mean-field algorithm; that
algorithm in turn can be approximated by particle methods. This philosophy
allows the conversion of standard single particle optimization and sampling
algorithms for inverse problems, with dependence on the derivative of the
forward model, into derivative-free interacting particle system optimizers and
samplers. An important application of this methodological approach is in the
development of ensemble Kalman approximations of Gauss-Newton and Levenberg-
Marquadt algorithms; pointers to the literature will be given in the
bibliography Subsection 4.6. We now turn to the use of statistical
linearization in gradient descent, perhaps the most basic setting in which it
can be used for optimization.

####  4.3.3 Gradient Descent and Statistical Linearization

We provide further insight into the statistical linearization approach from
the previous subsection by applying it in the context of gradient descent. We
consider the unregularized least squares function so that $\Psi(u)=\Phi(u)$ in
(143); but similar ideas may be developed in the regularized setting where
$\Psi(u)=\Phi_{R}(u).$ In the unregularized setting

(154) $\nabla\Phi(u)=DG(u)^{\top}\mathsf{\Gamma}^{-1}\bigl{(}G(u)-w\big{)},$

and we obtain, from (143),

$\frac{du}{dt}=DG(u)^{\top}\mathsf{\Gamma}^{-1}\bigl{(}w-G(u)\big{)}.$

Similarly, the covariance preconditioned gradient flow (149) becomes

(155) $\frac{du}{dt}=CDG(u)^{\top}\mathsf{\Gamma}^{-1}\bigl{(}w-G(u)\bigr{)},$

We now approximate this equation using statistical linearization.

From (153) we deduce the equivalent (assuming $C$ is invertible) approximation

$CDG(u)^{\top}\approx C^{uG},\quad\text{for all}\,u\in\mathbb{R}^{d_{u}}.$

Combining this with (154) we obtain

(156) $C\nabla\Phi(u)\approx
C^{uG}\mathsf{\Gamma}^{-1}\bigl{(}G(u)-w\bigr{)}.$

Making this approximation in (155) gives the following ensemble Kalman
approximation of mean field gradient descent:

(157) $\frac{du}{dt}=C^{uG}\mathsf{\Gamma}^{-1}\bigl{(}w-G(u)\bigr{)}.$

######  Remark 4.19.

We note that equation (157) is identical to (141) with the Brownian motion
simply set to zero. This provides a new perspective on Kalman transport for
inverse problems, namely that it may be thought of as a noisy version of an
approximation of mean-field gradient descent based on statistical
linearization; building on this idea we note that equation (142) may also be
viewed as a different approximation of mean-field gradient descent, also using
statistical linearization. Furthermore the Gaussian projected filter, in
continuous time, has a related structure which is apparent from comparing
(132) and (157).

It may be verified that the affine invariance of (149) is preserved under the
statistical linearization ansatz and, hence, (141), (142) and (157) are also
affine invariant. In fact, this property holds for all Kalman transport
formulations from Subsection 4.2. $\blacksquare$

We now observe that statistical linearization, and (156) in particular, become
exact for linear $G(u)=Lu$. In this setting we may hence study the asymptotic
behavior of the mean field gradient descent model (149) by studying (142) or
(157). To carry this out we study the deterministic continuous time Kalman
transport formulation (142), noting that it delivers perfect transport in the
linear setting. Furthermore, when initialized with
$u(0)\sim\mathsf{N}\bigl{(}m_{0},C_{0}\bigr{)},$ we know that $u(t)$ is
Gaussian distributed with mean $m(t)$ and covariance matrix $C(t)$ satisfying
the modified Kalman formulae (4.9), by Lemma 4.14. If we assume that $L$ is
invertible and hence that the desired minimizer $u^{\dagger}$ of $\Phi$ is
given by $u^{\dagger}=L^{-1}w$ then (4.9) reduce to

(158a) $\displaystyle m(t)$
$\displaystyle=m_{0}-C_{0}L^{\top}(LC_{0}L^{\top}+t^{-1}\mathsf{\Gamma})^{-1}L(m_{0}-u^{\dagger}),$
(158b) $\displaystyle C(t)$
$\displaystyle=C_{0}-C_{0}L^{\top}(LC_{0}L^{\top}+t^{-1}\mathsf{\Gamma})^{-1}LC_{0}.$

The following lemma characterizes the asymptotic behavior of the associated
mean field equations as $t\to\infty$.

######  Lemma 4.20.

Assume that $\mathsf{\Gamma}\succ 0$. Consider the mean-field ODE (142) and
assume further that $G(u)=Lu$ for invertible matrix
$L\in\mathbb{R}^{d_{u}\times d_{y}}$, so that $\Phi$ is given by

$\Phi_{L}(u)=\frac{1}{2}|Lu^{\dagger}-Lu|_{\mathsf{\Gamma}}^{2}$

with $u^{\dagger}=L^{-1}w$. Then the equation (142) for $u(t)$ may be written
in the form

(159)
$\frac{du}{dt}=-CL^{\top}\mathsf{\Gamma}^{-1}L\Bigl{(}\frac{1}{2}(u+m)-u^{\dagger}\Bigr{)},$

where $C$ is the covariance of $u(t)$ and $m$ its mean. Furthermore, if $u(0)$
is Gaussian then $u(t)$ is Gaussian with mean $m(t)$ and covariance $C(t)$
satisfying (158). Since $B:=L^{\top}\mathsf{\Gamma}^{-1}L\succ 0$, $u(t)$
converges to the mean $m(t)$ at an algebraic rate $t^{-1/2}$:
$d(t)=u(t)-m(t)$, satisfies

(160)
$\bigl{|}B^{1/2}d(t)\bigr{|}^{2}\leq\frac{1}{t}\bigl{|}B^{1/2}d(0)\bigr{|}^{2}.$

The mean, $m(t)$, converges to the minimizer $u^{\dagger}$ at an algebraic
rate $t^{-1}$:

(161)
$m(t)-u^{\dagger}=t^{-1}L^{-1}\mathsf{\Gamma}(LC_{0}L^{\top}+t^{-1}\mathsf{\Gamma})^{-1}L(m(0)-u^{\dagger}).$

######  Proof 4.21.

The mean $m(t)$ satisfies (158a) and straightforward algebraic manipulations
lead to (161). Furthermore, the mean field equation (159) and the
corresponding evolution equation for the mean

$\frac{dm}{dt}=-CL^{\top}\mathsf{\Gamma}^{-1}L(m-u^{\dagger})$

imply that

$\frac{d}{dt}(u-m)=-\frac{1}{2}CL^{\top}\mathsf{\Gamma}^{-1}L(u-m).$

Hence

(162)
$\frac{d}{dt}\bigl{|}B^{1/2}d\bigr{|}^{2}=-d^{\top}BCBd\leq-\frac{1}{t}\bigl{|}B^{1/2}d\bigr{|}^{2},$

where $d=u-m$ and we have used

$t^{-1}I\succ B^{1/2}CB^{1/2};$

this last matrix inequality follows from (158b) and its equivalent form

$C(t)^{-1}=C_{0}^{-1}+tB,$

which may be established using the Woodbury matrix identity. The estimate
(160) follows from the differential inequality (162) by integration.

######  Remark 4.22.

The same algebraic rates also apply to the mean field gradient descent
formulation (157). This algebraic rate of convergence is an undesirable
feature of algorithms based on the mean field model (157); algorithms based on
discretization of such continuous-time algorithms are not recommended in
practice unless used in combination with stochastic approximations to the cost
function $\Psi(u)$ (_stochastic gradient descent_), or if used with adaptive
time-stepping and a stopping criterion; see discussions in Subsection 4.6.
However by considering a stochastic analogue of the methodology we arrive, in
the next Subsection 4.3.4, at the mean field model (167); this has exponential
convergence rates to the posterior. Furthermore in Subsection 4.4.1 we will
also exhibit discrete-time algorithms for optimization which exhibit
exponential convergence. $\blacksquare$

####  4.3.4 Gradient Flow Structures in Spaces of Probability Measures

In this subsection, we investigate the geometric properties of the mean field
gradient descent equations (149) in some more detail. We assume that $u(t)$
has smooth probability density $\rho(u,t)$ for all $t\geq 0$ and denote the
manifold of all smooth probability density functions on $\mathbb{R}^{d_{u}}$
by ${\mathcal{P}}_{+}$. Then $\rho$ satisfies the Liouville equation

(163) $\partial_{t}\rho=\nabla\cdot(\rho C\nabla\Psi).$

Note that the appearance of the covariance matrix $C=\mathcal{C}(\rho)$
renders (163) a nonlinear and nonlocal partial differential equation on
${\mathcal{P}}_{+}$. We will show that the evolution of $\rho$ on
${\mathcal{P}}_{+}$ has gradient flow structure. In order to do so, we first
need to identify the functional $\mathcal{E}(\rho)$ on ${\mathcal{P}}_{+}$
which is being minimized. In case of (149), this functional is simply

$\mathcal{E}(\rho)=\int_{\mathbb{R}^{d_{u}}}\Psi(u)\rho(u)du,$

namely the expected value of $\Psi(u).$ Note that the Fréchet
derivative161616Identified by writing
$\mathcal{E}(\rho+\sigma)-\mathcal{E}(\rho)$ as a linear operator acting on
$\sigma\in T_{\rho}{\mathcal{P}}_{+}$ plus higher order terms. of
$\mathcal{E}$ is given by

$\frac{\delta\mathcal{E}}{\delta\rho}=\Psi.$

Hence we can rewrite (163) as

(164)
$\partial_{t}\rho=\nabla\cdot\Bigl{(}\rho\mathcal{C}(\rho)\nabla\frac{\delta\mathcal{E}}{\delta\rho}\Bigr{)}.$

The desired gradient flow structure on ${\mathcal{P}}_{+}$ requires an
appropriate metric tensor
$g_{\rho,\mathcal{C}}\,:\,T_{\rho}{\mathcal{P}}_{+}\times
T_{\rho}{\mathcal{P}}_{+}\to\mathbb{R}$ on the the tangent space
$T_{\rho}{\mathcal{P}}_{+}$ at $\rho\in{\mathcal{P}}_{+}$:

$\displaystyle T_{\rho}{\mathcal{P}}_{+}$ $\displaystyle=\left\\{\sigma\in
C^{\infty}(\mathbb{R}^{d_{u}})\,:\,\int_{\mathbb{R}^{d_{u}}}\sigma
du=0\right\\}\,.$

The metric tensor itself is defined by

$g_{\rho,\mathcal{C}}(\sigma_{1},\sigma_{2}):=\int_{\mathbb{R}^{d_{u}}}\left\langle\nabla\psi_{1}\,,\,\mathcal{C}(\rho)\nabla\psi_{2}\right\rangle\,\rho\,du,$

where
$\sigma_{i}=-\nabla\cdot\left(\rho\mathcal{C}(\rho)\nabla\psi_{i}\right)\in
T_{\rho}\mathcal{P}_{+}$ for $i=1,2$.171717A precise mathematical treatment
requires further assumptions on the considered set ${\mathcal{P}}_{+}$ of
probability functions ensuring that the required potentials $\psi$ indeed
exists for all $\sigma\in T_{\rho}{\mathcal{P}}_{+}$. See the bibliography for
references and discussion of open questions in this area. The equation (164)
takes gradient flow structure in $(\mathcal{P}_{+},g_{\rho,\mathcal{C}})$:

(165)
$\frac{d}{dt}\mathcal{E}(\rho)=-g_{\rho,\mathcal{C}}(\partial_{t}\rho,\partial_{t}\rho).$

It is interesting to compare the gradient structure (165) on
${\mathcal{P}}_{+}$ to the gradient flow structure (4.3.1) on
$\mathbb{R}^{d_{u}}$. To establish (165) note that

$\displaystyle\frac{d}{dt}\mathcal{E}(\rho)$
$\displaystyle=\int_{\mathbb{R}^{d_{u}}}\frac{\delta\mathcal{E}}{\delta\rho}\partial_{t}\rho\,du\,$
$\displaystyle=-\int_{\mathbb{R}^{d_{u}}}\rho\,\Bigl{|}\mathcal{C}(\rho)^{\frac{1}{2}}\nabla\frac{\delta\mathcal{E}}{\delta\rho}\Bigr{|}^{2}\,du\,$
$\displaystyle=-\int_{\mathbb{R}^{d_{u}}}\left\langle\nabla\psi\,,\,\mathcal{C}(\rho)\nabla\psi\right\rangle\,\rho\,du\,$
$\displaystyle=-g_{\rho,\mathcal{C}}(\sigma,\sigma)\,,$

where

$\sigma=-\nabla\cdot\left(\rho\mathcal{C}(\rho)\nabla\psi\right),\quad\psi=\frac{\delta\mathcal{E}}{\delta\rho}.$

Thus $\sigma=-\partial_{t}\rho$ and the gradient flow identity is established.

The established machinery is very powerful and demonstrates that any evolution
equation of type (164) with appropriate potential $\mathcal{E}$ induces a
gradient flow on ${\mathcal{P}}_{+}.$ A particularly significant example
follows from the choice

$\mathcal{E}(\rho)=\int\left(\Phi_{R}+\ln\rho\right)\rho\,du$

with Fréchet derivative given by 181818Formal calculations give rise to an
additional constant term which is irrelevant and in fact zero if the Fréchet
derivative is computed using $\sigma\in T_{\rho}{\mathcal{P}}_{+}$.

$\frac{\delta\mathcal{E}}{\delta\rho}=\Phi_{R}+\ln\rho.$

In this case the associated evolution equation (164) becomes

(166)
$\partial_{t}\rho=\nabla\cdot(\rho\mathcal{C}(\rho)\nabla\Phi_{R})+\nabla\cdot(\mathcal{C}(\rho)\nabla\rho).$

The mean-field SDE model for which this is the nonlinear and nonlocal Fokker-
Planck equation is

(167) $du=-C\nabla\Phi_{R}(u)dt+\sqrt{2C}dW,$

where we reverted to the notation $C$, for the covariance, as used previously
in all mean field models. Theory concerning this equation is contained in the
following bibliography section. We note here that it may be verified that
(167) is invariant under affine transformations of type (144) and that the
posterior distribution $\rho^{w}$ is stationary, that is, $\rho=\rho^{w}$
leads to $\partial_{t}\rho=0$ in (166). Thus (167) provides an attractive
generalization of standard Langevin dynamics for sampling from the posterior
distribution; furthermore it may be approximated by statistical linearization
to obtain a derivative-free methodology. For (166) the desired posterior
distribution is approached only in the limit $t\to\infty$. The relative
efficiency of methods based on this approach, in comparison with the transport
methods considered in Subsection 4.2, depends on the nature of this
convergence; the fact that it is exponential in many settings make the
approach potentially competitive. However discrete time versions of the
approach of iterating mean field models to infinity can exhibit greater
competitiveness and are explored in the following subsection.

We conclude this subsection by noting that the Kalman-Wasserstein gradient
flow structure for the Liouville and Fokker-Planck equations considered here
is not maintained under statistical linearization: they are not of gradient
descent type in ${\mathcal{P}}_{+}$ unless $G(u)$ is linear.

###  4.4 Steady State Filtering Distribution and Inversion

The transport of prior to posterior in $N$ steps, or one time unit, the
subject of Subsection 4.2, provides an attractive framework within which to
derive algorithms for approximate Bayesian inversion from state estimation and
filtering of a dynamical system. However it is a somewhat rigid approach,
requiring specific initialization of the variable $u_{0}$ at the prior; and it
does not possess any strong stability property that might be used to control
the errors arising from using second order, rather than exact, transport, or
inexact initialization. On the other hand, Subsection 4.3 shows that Kalman-
inspired methods may be used to solve optimization and sampling problems,
which may be initialized arbitrarily but require iteration to infinity to
address their respective goals. In this subsection we further explore the idea
of iterating to infinity with Kalman methods. As in Subsection 4.2, the
methods also employ state estimation and filtering of a dynamical system.
However, instead of transporting an initial prior to the posterior at a fixed
time (discrete or continuous), we instead seek an underlying partially and
noisily observed dynamical system with filtering distribution that has as
steady state a solution to the inverse problem; in brief we work in settings
where the filtering distribution is ergodic and use the statistical steady
state to extract information about the inverse problem.

Subsections 4.4.1 and 4.4.2 approach this problem from the optimization
perspective on inversion, and provide an intuitive introduction to the use of
ergodic filtering distributions. However the approach we adopt does not lead
to affine invariant methods and so, in Subsections 4.4.3 and 4.4.4 we
introduce a different filtering distribution. Identifying filtering
distributions with statistical steady state which solves the inverse problem
is difficult in general, and so we build the methodology by adopting the more
limited goal of using a dynamical system which achieves this goal exactly in
the linear Gaussian setting; we then deploy it beyond the linear Gaussian
setting. The approach is consistent with the use of second order transport,
rather than exact transport, underlying the application of ensemble Kalman
methods: second order transport is also exact only in the linear Gaussian
setting, but is used beyond this setting. In the first two subsections
concerned with optimization the underlying partially and noisily observed
dynamical system used is, like (1) or (75), not of mean-field form. However,
in the subsequent two subsections concerned with the Bayesian approach, the
underlying partially and noisily observed dynamical system used is itself of
mean-field form; indeed it is an unusual mean-field model that depends on its
own filtering distribution and it is this fact that facilitates the derivation
of an affine invariant methods. Note that in Sections 2 and 3 respectively,
the underlying dynamical system ((1) or (75) respectively), to which ensemble
methods are applied, was not of mean-field form: mean-field was only
introduced to construct ensemble Kalman methods.

####  4.4.1 Steady State Filtering Distribution and Classical Inversion:
Discrete Time

With the goal of solving least-squares based optimization problems using ideas
from filtering, consider now the stochastic dynamical system

(168a) $\displaystyle u_{n+1}$ $\displaystyle=\alpha
u_{n}+(1-\alpha)r_{0}+\xi_{n},$ (168b) $\displaystyle w_{n+1}$
$\displaystyle=G(u_{n+1})+\eta_{n+1}.$

Here $\alpha\in(0,1]$,
$\xi_{n}\sim\mathsf{N}(0,\sigma^{\prime}\mathsf{\Sigma})$,
$\sigma^{\prime}\geq 0$ and
$\eta_{n+1}\sim\mathsf{N}(0,\gamma^{\prime}\mathsf{\Gamma})$,
$\gamma^{\prime}>0$; we assume that $\mathsf{\Sigma},\mathsf{\Gamma}\succ 0$.
Unlike the transport setting of Subsection 4.2 we do not require that $u_{0}$
is initialized at the prior, although this is a reasonable choice of
initialization. Let $W_{n}=\\{w_{\ell}\\}_{\ell=1}^{n}$ and consider the
filtering distribution defined by random variable $u_{n}|W_{n}.$ The following
theorem provides motivation for the study of this distribution, and its
limiting behavior as $n\to\infty$, to obtain information about the least
squares optimization approach to the inverse problem.

######  Lemma 4.23.

Consider the filtering distribution $u_{n}|W_{n}$ of (168), and assume that
$u_{0}$ is initialized at a Gaussian $\mathsf{N}(m_{0},C_{0}).$191919For
simplicity of exposition we denote this as $\mathsf{N}(m_{0},C_{0})$ but
emphasize that here it need not be the prior distribution. Assume that
$C_{0},\mathsf{\Gamma},\mathsf{\Sigma}\succ 0$, $\sigma^{\prime}\geq 0$,
$\gamma^{\prime}>0$, and $\alpha\in(0,1]$ and and consider the setting where
$G(u)=Lu$ for matrix $L\in\mathbb{R}^{d_{w}\times d_{u}}$. Then the filtering
distribution is Gaussian $\mathsf{N}(m_{n},C_{n}).$ Furthermore, if
$w_{n}\equiv w$ then any positive-definite steady solution $C_{\infty}$ of the
covariance update equation satisfies

(169a) $\displaystyle C_{\infty}^{-1}$
$\displaystyle=\frac{1}{\gamma^{\prime}}L^{\top}\mathsf{\Gamma}^{-1}L+\widehat{C}_{\infty}^{-1},$
(169b) $\displaystyle\widehat{C}_{\infty}$
$\displaystyle=\alpha^{2}C_{\infty}+\sigma^{\prime}\mathsf{\Sigma},$

where, note, $\widehat{C}_{\infty}$ is also positive-definite. Furthermore,
the associated steady state mean $m_{\infty}$ is the minimizer of the
Tikhonov-Phillips regularized least squares functional $\Phi_{TP}$ given by

$\Phi_{TP}(m):=\frac{1}{2\gamma^{\prime}}|w-Lm|_{\mathsf{\Gamma}}^{2}+\frac{1-\alpha}{2}|m-r_{0}|_{\widehat{C}_{\infty}}^{2}.$

######  Proof 4.24.

We first identify the update equations for $(m_{n},C_{n})$. Note that the
predictive mean $\widehat{m}_{n+1}$ and covariance $\widehat{C}_{n+1}$ defined
by (168a) satisfy

$\displaystyle\widehat{m}_{n+1}$ $\displaystyle=\alpha m_{n}+(1-\alpha)r_{0},$
$\displaystyle\widehat{C}_{n+1}$
$\displaystyle=\alpha^{2}C_{n}+\sigma^{\prime}\mathsf{\Sigma}.$

To find $(m_{n+1},C_{n+1})$ it is possible to simply use (34); however it is
convenient to derive the formulae in a different form, using precisions rather
than covariance. To this end we view the Gaussian
$\mathsf{N}(\widehat{m}_{n+1},\widehat{C}_{n+1})$ as prior distribution for
the linear inverse problem defined by (168b) with $w_{n+1}=w.$ Note that the
likelihood, since linear and Gaussian, is conjugate to the prior so that the
posterior on $u_{n+1}|W_{n+1}$ is Gaussian with mean and covariance
$(m_{n+1},C_{n+1})$ which can be found by completing the square:

$\displaystyle C_{n+1}^{-1}$
$\displaystyle=\frac{1}{\gamma^{\prime}}L^{\top}\mathsf{\Gamma}^{-1}L+\widehat{C}_{n+1}^{-1},$
$\displaystyle C_{n+1}^{-1}m_{n+1}$
$\displaystyle=\frac{1}{\gamma^{\prime}}L^{\top}\mathsf{\Gamma}^{-1}w+\widehat{C}_{n+1}^{-1}\widehat{m}_{n+1}.$

The equation satisfied by the steady state covariance $C_{\infty}$ follows,
along with the definition of $\widehat{C}_{\infty}$. The steady state mean
satisfies

$C_{\infty}^{-1}m_{\infty}=\frac{1}{\gamma^{\prime}}L^{\top}\mathsf{\Gamma}^{-1}w+\widehat{C}_{\infty}^{-1}\bigl{(}\alpha
m_{\infty}+(1-\alpha)r_{0}\bigr{)}.$

Using the equation (169a) to replace $C_{\infty}$ by $\widehat{C}_{\infty}$ we
may rewrite this as

$\widehat{C}_{\infty}^{-1}m_{\infty}=\frac{1}{\gamma^{\prime}}L^{\top}\mathsf{\Gamma}^{-1}(w-Lm_{\infty})+\widehat{C}_{\infty}^{-1}\bigl{(}\alpha
m_{\infty}+(1-\alpha)r_{0}\bigr{)}.$

This in turn may be reorganized as

$(1-\alpha)\widehat{C}_{\infty}^{-1}(m_{\infty}-r_{0})=\frac{1}{\gamma^{\prime}}L^{\top}\mathsf{\Gamma}^{-1}(w-Lm_{\infty})$

which is patently the equation found by setting the gradient of $\Phi_{TP}$ to
$0.$ Finally note that the Hessian of $\Phi_{TP}$ is
$(1-\alpha)\widehat{C}_{\infty}^{-1}+\frac{1}{\gamma^{\prime}}L^{\top}\mathsf{\Gamma}^{-1}L\succ
0.$ Thus the critical point is indeed a minimizer as asserted.

######  Remark 4.25.

Assume either that $\alpha\in(0,1)$ or that, if $\alpha=1$ then assume
additionally that the range of $G^{\top}$ is $\mathbb{R}^{d_{u}}$, the entire
parameter space, and that the noise in the prediction is positive
$\sigma^{\prime}>0$. Then, although we do not prove it here, the equilibrium
solutions of the filtering distribution are, in the linear setting, unique
(within the class of solutions with positive-definite covariance) and
exponentially attracting. This gives methodology based on application of
filtering to (168) a robustness that the transport methods of Subsection 4.2
do not possess. We give details in the bibliography Subsection 4.6. These
results suggest that we should use $\alpha\in(0,1)$ except in situations where
the number of observations is greater than or equal to the dimension of the
parameter space $\mathbb{R}^{d_{w}}\geq\mathbb{R}^{d_{u}}$ when it is
plausible to use $\alpha=1;$ indeed empirical evidence (see papers cited in
Subsection 4.6) suggests that such a choice often works well for over-
determined inverse problems. $\blacksquare$

######  Remark 4.26.

The Tikhonov-Phillips regularized least squares function $\Phi_{TP}$
corresponds to regularization of the inverse problem (117), with
$G(\cdot)=L\cdot$, stemming from Gaussian prior
$\mathsf{N}\bigl{(}r_{0},(1-\alpha)^{-1}\widehat{C}_{\infty}\bigr{)}.$ This is
non-standard from the point of view of inverse problems because the prior
covariance is proportional to $\widehat{C}_{\infty}$ which depends on the
forward model $L.$ $\blacksquare$

Returning to the general nonlinear inverse problem (168), and applying the
Kalman transport mean field model (46) we obtain

(170a) $\displaystyle\widehat{u}_{n+1}$ $\displaystyle=\alpha
u_{n}+(1-\alpha)r_{0}+\xi_{n},$ (170b) $\displaystyle\widehat{w}_{n+1}$
$\displaystyle=G(u_{n+1})+\eta_{n+1},$ (170c) $\displaystyle u_{n+1}$
$\displaystyle=\widehat{u}_{n+1}+\widehat{C}^{uG}_{n+1}(\widehat{C}^{GG}_{n+1}+\gamma^{\prime}\mathsf{\Gamma})^{-1}(w-\widehat{w}_{n+1}).$

Here

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

and, recall, $\alpha\in(0,1)$,
$\xi_{n}\sim\mathsf{N}(0,\sigma^{\prime}\mathsf{\Sigma})$ and
$\eta_{n+1}\sim\mathsf{N}(0,\gamma^{\prime}\mathsf{\Gamma}).$

####  4.4.2 Steady State Filtering Distribution and Classical Inversion:
Continuous Time

We now study the inverse problem using steady state optimization in a
continuous time limit. We set $\alpha=1-a\Delta t$, $\sigma^{\prime}=\Delta t$
and $\gamma^{\prime}=\Delta t^{-1}$ in (168) to obtain the continuous time
limit

(171a) $\displaystyle du$
$\displaystyle=-a(u-r_{0})dt+\sqrt{\mathsf{\Sigma}}dW,$ (171b) $\displaystyle
dz$ $\displaystyle=G(u)dt+\sqrt{\mathsf{\Gamma}}dB$

where $W,B$ are independent unit Brownian motions of appropriate dimensions.
We are interested in the filtering distribution $u(t)|{Z}^{\dagger}(t)$ in the
case where $\frac{d{z}^{\dagger}}{dt}\equiv w.$

Taking the continuous time limit in equation (170), we obtain the following
mean field Kalman transport model for the inverse problem:

(172a) $\displaystyle du$
$\displaystyle=-a(u-r_{0})dt+C^{uG}\mathsf{\Gamma}^{-1}(w\,dt-d\widehat{z})+\sqrt{\mathsf{\Sigma}}dW,$
(172b) $\displaystyle d\widehat{z}$
$\displaystyle=G(u)dt+\sqrt{\mathsf{\Gamma}}dB,$ (172c) $\displaystyle C^{uG}$
$\displaystyle=\mathbb{E}\Bigl{(}\bigl{(}u-\mathbb{E}u\bigr{)}\otimes\bigl{(}G(u)-\mathbb{E}G(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. The resulting
algorithm nudges $G(u)$ towards the data, does so in direction informed by
$C^{uG}$, and has additional noise and reversion to a prescribed point
$r_{0}.$ It is of interest to consider this continuous time limit in the
context of Remark 4.19. It too has at its core an approximate mean field
gradient descent, defined via statistical linearization. However the
additional noise term in state-space, and the reversion to $r_{0}$, lead to a
loss of affine invariance.

We now generalize the ideas in the previous subsection from the optimization
approach to the Bayesian approach to inversion. Following the same philosophy
as the last subsection, we first identify a dynamical system with steady state
filtering distribution given by solution of the Bayesian inverse problem. As
in the optimization setting of the preceding subsection this is difficult in
general, but may be achieved for linear Gaussian inverse problems. And, as in
the optimization setting, we can then use the approach in the nonlinear
setting noting that this invokes an approximation that is consistent with the
use of second order, rather than exact, transport to derive the mean-field
models. Since these mean-field models are the basis for our ensemble Kalman
algorithms this approach to inversion is based on a self-consistent set of
approximations. The resulting methodology leads to affine invariant dynamics.

####  4.4.3 Steady State Filtering Distribution and Bayesian Inversion:
Discrete Time

The previous two subsections study the idea of iterating to infinity to solve
optimization approaches to the inverse problem. The methods do not solves the
Bayesian inverse problem defined by (119), (4.1) and, furthermore, are not
affine invariant. In this and the next subsection we address these two issues,
generalizing the ideas of the preceding two subsections. To this end, consider
the dynamical system

(173a) $\displaystyle u_{n+1}$ $\displaystyle=u_{n}+\xi_{n},$ (173b)
$\displaystyle w_{n+1}$ $\displaystyle=G_{R}(u_{n+1})+\eta_{n+1}.$

Let $W_{n}=\\{w_{\ell}\\}_{\ell=1}^{n}$ and consider the filtering
distribution defined by random variable $u_{n}|W_{n},$ and let $C_{n}$ denote
the covariance of this random variable. Choose $\alpha>0$ and asssume that
$\xi_{n}\sim\mathsf{N}(0,\alpha C_{n})$ and that
$\eta_{n+1}\sim\mathsf{N}\bigl{(}0,(1+\alpha^{-1})\mathsf{\Gamma}_{R}).$ Note
a _fundamental difference with all the other state-observation models_
considered in this paper: this is a mean-field model. Everywhere else in the
paper we only introduce mean-field models when we construct transport maps for
approximate filtering. Furthermore, (173) with $\xi_{n}\sim\mathsf{N}(0,\alpha
C_{n})$ and $C_{n}$ as defined is a non-standard mean-field stochastic
dynamical system which depends not on its own law, but rather on the law of
its own filtering distribution $\mu_{n}$. Although the filtering distribution
depends on the history $W_{n}$, equation (173) can be rendered Markovian by
coupling it to the Markovian evolution of the filtering distribution
$\mu_{n}\mapsto\mu_{n+1}$, and noting that $C_{n}$ is computed under
$\mu_{n}.$

Recall $w_{R}$ defined by (119). The non-standard mean-field dynamical system
we have just introduced has a highly desirable property in the linear setting:

######  Lemma 4.27.

Assume that $C_{0}\succ 0$ and $\mathsf{\Gamma}\succ 0$. Consider the
filtering distribution $u_{n}|W^{\dagger}_{n}$ for the mean-field dynamical
system (173) in the setting where $w^{\dagger}_{n}\equiv w_{R}$ and in the
setting where $u_{0}$ is initialized at a Gaussian.202020 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_{n}|W^{\dagger}_{n}$ is Gaussian
$\mathsf{N}(m_{n},C_{n}).$ Then the steady state mean $m_{R}$ and precision
matrix $C^{-1}_{R}$ of the filter are those of the posterior distribution
(118) arising from the Bayesian inverse problem defined by (117), in this
linear setting.

######  Proof 4.28.

Under the dynamics (173a) the predictive mean and covariance are given by

$\widehat{m}_{n+1}=m_{n},\quad\widehat{C}_{n+1}=(1+\alpha)C_{n}.$

Taking $(\widehat{m}_{n+1},\widehat{C}_{n+1})$ as prior for the inverse
problem defined by (173b) we deduce that, by completing the square analogously
to the same step in the proof of Lemma 4.23,

$\displaystyle C_{n+1}^{-1}$
$\displaystyle=\widehat{C}_{n+1}^{-1}+\frac{\alpha}{1+\alpha}L_{R}^{\top}\mathsf{\Gamma}_{R}^{-1}L_{R},$
$\displaystyle C_{n+1}^{-1}m_{n+1}$
$\displaystyle=\widehat{C}_{n+1}^{-1}\widehat{m}_{n+1}+\frac{\alpha}{1+\alpha}L_{R}^{\top}\mathsf{\Gamma}_{R}^{-1}w_{R};$

here $L_{R}$ is defined in (120). The equilibrium covariance is defined by