$\displaystyle=\mathopen{}\mathclose{{}\left(\bar{w}^{(n)}_{j}}\right)^{2}\mathopen{}\mathclose{{}\left[\int_{\mathbb{R}^{d}}\frac{\mathopen{}\mathclose{{}\left(\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1})\operatorname{\mathbb{P}}(v_{j+1}|v^{(n)}_{j})}\right)^{2}}{\pi(v_{j+1}|v_{j}^{(n)},Y_{j+1})}dv_{j+1}-\operatorname{\mathbb{P}}(y_{j+1}|v^{(n)}_{j})^{2}}\right].$

Choosing
$\pi(v_{j+1}|v_{j}^{(n)},Y_{j+1})=\operatorname{\mathbb{P}}(v_{j+1}|v_{j}^{(n)},y_{j+1})$,
as in the OPF, we obtain

$\displaystyle{\rm
Var}_{\operatorname{\mathbb{P}}(v_{j+1}|v_{j}^{(n)},Y_{j+1})}[\bar{w}^{(n)}_{j+1}]$
$\displaystyle=\mathopen{}\mathclose{{}\left(\bar{w}^{(n)}_{j}}\right)^{2}\mathopen{}\mathclose{{}\left[\int_{\mathbb{R}^{d}}\frac{\mathopen{}\mathclose{{}\left(\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1})\operatorname{\mathbb{P}}(v_{j+1}|v^{(n)}_{j})}\right)^{2}}{\operatorname{\mathbb{P}}(v_{j+1}|v_{j}^{(n)},y_{j+1})}dv_{j+1}-\operatorname{\mathbb{P}}(y_{j+1}|v^{(n)}_{j})^{2}}\right]$
$\displaystyle=\mathopen{}\mathclose{{}\left(\bar{w}^{(n)}_{j}}\right)^{2}\mathopen{}\mathclose{{}\left[\operatorname{\mathbb{P}}(y_{j+1}|v^{(n)}_{j})^{2}-\operatorname{\mathbb{P}}(y_{j+1}|v^{(n)}_{j})^{2}}\right]$
$\displaystyle=0.$

∎

######  Remark 12.2.

The optimal particle filter is optimal in the very precise sense of the
theorem. Note in particular that no optimality criterion is asserted by this
theorem with respect to iterating the particle updates, and in particular when
resampling is included. The nomenclature “optimal” should thus be treated with
caution.

####  12.5 Particle Filters for High Dimensions

Particle filters often perform poorly for high-dimensional systems due to a
collapse of the particle weights: only a few particles carry most of the
weight; the optimal particle filter can ameliorate this issue because the
proposal uses the data, meaning that particle predictions are more likely to
be weighted highly. There have been attempts to formulate update steps that
help to mitigate this weight collapse. Essentially, these methods aim to push
the particles towards the region of high likelihood, such that all the
particles will be representative of the distribution. As will show in the next
lecture, there are interesting particle based methods which, whilst not
statistically consistent, do perform well as signal estimators. This
perspective on particle filter type methods, namely to use them for smart
signal estimation rather than as estimators of the filtering distribution, may
become increasingly useful in high dimensional systems.

####  12.6 Filtering Overview

Having introduced a range of filtering methods in this and the preceding four
chapters, it is helpful to summarize what these methods achieve, in a
comparitive fashion.

#####  12.6.1 Dynamical Model

The stochastic dynamics model is defined by

$v_{j+1}=\Psi(v_{j})+\xi_{j},\quad j\geq 0,$

where $\xi_{j}\sim\mathcal{N}(0,\Sigma)$ are independent and identically
distributed random variables, also independent of the initial condition
$v_{0}\sim\mathcal{N}(m_{0},C_{0})$.

Figure 12: Prediction step.

The filtering distribution of $v_{j}$ is denoted by
$\pi_{j}(v_{j})=\operatorname{\mathbb{P}}(v_{j}|Y_{j}).$ This is propagated by
the stochastic dynamics model according to the formula
$\widehat{\pi}_{j+1}={\mathcal{P}}\pi_{j}$ where ${\mathcal{P}}$ is defined in
(11.1). This is shown schematically in Figure 12. Note that ${\mathcal{P}}$ is
independent of step $j$ because the Markov chain defined by the stochastic
dynamics model is time-homogeneous. In the absence of data the probability
distribution simply evolves through repeated application of ${\mathcal{P}}$.

#####  12.6.2 Data Model

The data model is given by

$y_{j+1}=h(v_{j+1})+\eta_{j+1},$

where $\eta_{j+1}\sim{N}(0,\Gamma)$ are independent and identically
distributed random variables, independent of both $v_{0}$ and the i.i.d.
sequence $\\{\xi_{j}\\}.$ Given data, we now view the prediction from the
dynamical model as a prior which we condition on the data. We write this prior
at time $j+1$ as $\widehat{\pi}_{j+1}$, and prediction using the stochastic
dynamics model gives $\widehat{\pi}_{j+1}={\mathcal{P}}\pi_{j}.$

Using the data to update this prior by Bayes theorem gives an improved
estimate of the distribution $\pi_{j+1}(v_{j+1})$ via the formula

$\pi_{j+1}={\mathcal{A}}_{j}\widehat{\pi}_{j+1}$

where ${\mathcal{A}}_{j}$ is given by (11.2). This is shown schematically in
Figure 13. Note that ${\mathcal{A}}_{j}$ is nonlinear because
$\widehat{\pi}_{j+1}$ appears in both the numerator and the denominator. It
depends on $j$ because the data $y_{j+1}$ appears in the equation and this
will change with each set of measurements.

Figure 13: Update step.

#####  12.6.3 Interaction of the Dynamics and Data

Recall the notation $Y_{j}=\\{y_{1},\dots,y_{j}\\}$ for the collection of all
data until time $j+1$. The dynamics and data models described in the two
previous subsections can be combined as follows. Using the dynamics, the
probability distribution is propagated forward in time by the prediction step.
The data is then used to update the probability distribution at that time in
the analysis step. The two can be combined to obtain

${\pi}_{j+1}={\mathcal{A}}_{j}{\mathcal{P}}\pi_{j}.$

The combination of the two steps is shown schematically in Figure 14.

Figure 14: Prediction and update step combined.

#####  12.6.4 Summary of Applicability of Discrete Filtering Methods

There are several filtering methods for performing the prediction and analysis
steps. Some methods, such as the bootstrap particle filter, can be applied
generally to nonlinear problems. However, others require a linear model
$(\Psi(\cdot)=M\cdot)$ and/or linear observations $(h(\cdot)=H\cdot)$. Some of
the methods we have described provably approximate the probability
distribution updates. Some just estimate the state and simply use covariance
information to weight the relative importance of the predictions from the
model and of the data.

The applicability of the methods introduced is summarized in the table below,
with respect to linearity/nonlinearity of the dynamics and the observation
model. Furthermore ${\bf P}$ is used to denote methods which provably
approximate the filtering distributions $\pi_{j}$ in the large particle limit;
${\bf S}$ denotes methods which only attempt to estimate the state, using the
data.

Kalman Filter

$\Psi(\cdot)=M\cdot,h(\cdot)=H\cdot$. ${\bf P}.$

3DVAR

$\text{general}\,\Psi,h(\cdot)=H\cdot$. ${\bf S}.$

Bootstrap Particle Filter

$\text{general}\,\Psi,h$. ${\bf P}.$

Optimal Particle Filter

$\text{general}\,\Psi,h(\cdot)=H\cdot$. ${\bf P}$.

Extended Kalman Filter

$\text{general}\,\Psi,h(\cdot)=H\cdot$. ${\bf S}.$

Ensemble Kalman Filter

$\text{general}\,\Psi,h(\cdot)=H\cdot$. ${\bf S}.$

Some of these constraints can be relaxed on the setting in which they apply
can be relaxed, but the list above describes the methods as we present them in
these notes. Furthermore the last two methods can accurately predict
probability distributions in situations where approximate Gaussianity holds;
this may be induced by small noise and by large data. However it is important
to appreciate that the raison d’etre of the EnKF is to facilitate the solution
of problems with high dimensional state space; typically small ensemble sizes
are used and the algorithm is employed far from the regime in which it is able
to provably approximate the distributions, even if they are close to Gaussian.
It is for this reason that we prefer to think of the EnKF as an ensemble state
estimator.

####  12.7 Discussion and Bibliography

Particle filters often perform poorly for high-dimensional systems due to the
fact that the particle weight typically concentrates on one, or a small
number, of particles –see the work of Bickel and Snyder in [12, 107, 106].
This is the issue that the optimal particle filter tries to ameliorate; the
paper [4] shows calculations which demonstrate the extent to which this
amelioration is manifest in theory. The optimal particle filter is discussed,
and further references given, in the very clear paper [28]; see section IID.
Throughout much of this chapter we considered the case of Gaussian additive
noise and linear observation operator in which case the prediction step is
tractable; the paper [28] discusses the more general setting. The order in
which the prediction and resampling is performed can be commuted in this case
and a discussion of this fact may be found in [95]; this leads to the
distinction between what we term the GOPF and the OPF. The convergence of the
optimal particle filter is studied in [62]. The formulation of the bootstrap
and optimal particle filters as random dyamical systems may be found in [70].
An attempt to alleviate weight collapse by introducing alternative update
steps can be found in [106].

###  13 Filtering Approach to the Inverse Problem

In this final chapter we demonstrate how the two separate themes that underpin
this course, inverse problems and data assimilation, may be linked. This opens
up the possibility of transferring ideas from filtering into the setting of
quite general inverse problems. In the first section we describe the general,
abstract, connection and introduce the revolutionary idea of sequential Monte
Carlo methods (SMC). In the second section we analyze the concrete case of
applying the EnKF to solve an inverse problem, leading to ensemble Kalman
inversion (EKI). In the final section we link the EKI methodology to SMC.

####  13.1 General Formulation

Recall the inverse problem of finding $u\in\mathbb{R}^{d}$ from
$y\in\mathbb{R}^{k}$ where

$\displaystyle y=G(u)+\eta,\;\eta\sim\mathcal{N}(0,\Gamma_{0})$ (13.1)

and the related loss function

$\mathsf{L}_{0}(u)=\frac{1}{2}|y-G(u)|^{2}_{\Gamma_{0}}.$

The reason for writing $\Gamma_{0}$, rather than $\Gamma$, will become
apparent below and will also be exploited in the third section of this chapter
where we study EKI.

If we put a prior $\rho$ on the unknown $u$ then the posterior takes the form

$\pi(u)=\frac{1}{Z}\exp\bigl{(}-\mathsf{L}_{0}(u)\bigr{)}\rho(u).$

Let $J\in\mathbb{N}$ and choose $h$ so that $Jh=1$. Then define the family of
pdfs $\\{\pi_{j}\\}_{j=0}^{J}$ by

$\pi_{j}(u)=\frac{1}{Z_{j}}\exp\bigl{(}-jh\mathsf{L}_{0}(u)\bigr{)}\rho(u).$

It follows that $\pi_{0}=\rho$ and $\pi_{J}=\pi$ and, furthermore, we may
update the sequence $\\{\pi_{j}\\}_{j=0}^{J}$ sequentially using the formula

$\pi_{j+1}(u)=\frac{Z_{j}}{Z_{j+1}}\exp\bigl{(}-h\mathsf{L}_{0}(u)\bigr{)}\pi_{j}(u).$

This simply corresponds to application of Bayes theorem to the inverse problem

$\displaystyle y=G(u)+\eta,\;\eta\sim\mathcal{N}(0,\Gamma),$ (13.2)

with $\Gamma=\frac{1}{h}\Gamma_{0}$ and with prior $\pi_{j}.$ With this choice
of $\Gamma$ the identity (13.2) gives the likelihood proportional to
$\exp(-h\mathsf{L}_{0}(u))$. The update may be written

$\pi_{j+1}(u)={\mathcal{A}}\pi_{j},$

noting that the analysis operator ${\mathcal{A}}$ is nonlinear, corresponding
to multiplication by $\exp(-h\mathsf{L}_{0}(u))$ and then normalization to a
pdf.

If we let ${\mathcal{P}}_{j}$ denote any Markov kernel for which $\pi_{j}$ is
invariant then we obtain

$\pi_{j+1}={\mathcal{A}}{\mathcal{P}}_{j}\pi_{j}.$

This update formula should be compared with the filtering update formula

$\displaystyle\pi_{j+1}={\mathcal{A}}_{j}{\mathcal{P}}\pi_{j}$

introduced and used in earlier chapters.

Whilst the $j-$dependence in Bayes rule and the Markov kernel is interchanged
between the inverse problem and the filtering problem, this makes little
material difference to implementation. Firstly, once a Markov kernel
${\mathcal{P}}$ is identified under which $\pi=\pi_{J}$ is invariant, it can
be easily adapted to find a family of kernels ${\mathcal{P}}_{j}$ under which
$\pi_{j}$ is invariant, simply by rescaling the observation covariance.
Secondly the fact that the data $y$ is fixed, rather than changing at each
step, makes little difference to the implementation. This perspective opens up
an entire field, known as sequential Monte Carlo, in which ideas from
filtering may be transferred to other quite different problems, including
Bayesian inversion, as explained here.

Note however (see the discussion at the end of section 12.6) that of the
filtering methods we have introduced so far only the bootstrap filter applies
directly to the case of nonlinear observation operators; since $G$ is in
general nonlinear this means that extra ideas are required to implement the
optimal particle filter, 3DVAR, ExKF and EnKF. One approach is to replace the
sequential optimization principles by the non-quadratic optimization problem
required for nonlinear $G$; we do not discuss this idea in any detail but it
is a viable option. Another approach is to use the linearization technique
which we now outline in the context of EKI.

####  13.2 Ensemble Kalman Inversion

We now make a detour into the subject of how to use filters to estimate
parameters $u$ from data $y$ satisfying (13.2). The approach we study, and
which we will relate to sequential Monte Carlo in the next section, is to
introduce an artificial time dynamic. This can be done quite simply as
follows: we write

$\displaystyle u_{j+1}$ $\displaystyle=u_{j},$ $\displaystyle y_{j+1}$
$\displaystyle=G(u_{j+1})+\eta_{j+1},$

and we can think of finding the filtering distribution on $u_{j}|Y_{j}.$ We
discuss how to relate the $y_{j}$ to the one data point $y$ below. For now we
take $\eta_{j}\sim\mathcal{N}(0,\Gamma)$ but we will revisit this choice in
the next section, linking the problem to solution of (13.1). Because the
observation operator $G$ is, in general, nonlinear, this does not render our
system in a form where we can readily apply the EnKF. To this end we introduce
a new variable $w_{j}$ and rewrite the filter as:

$\displaystyle u_{j+1}$ $\displaystyle=u_{j},$ $\displaystyle w_{j+1}$
$\displaystyle=G(u_{j}),$ $\displaystyle y_{j+1}$
$\displaystyle=w_{j+1}+\eta_{j+1}.$

We introduce the new variables $v=(u,w)^{T}$, nonlinear map
$\Psi(v)=\bigl{(}u,G(u)\bigr{)}^{T}$ and linear operators
$H=[0,I],H^{\perp}=[I,0].$ Then if we write $v_{j}=(u_{j},w_{j})^{T}$ we may
write the dynamical system in the form

$\displaystyle v_{j+1}$ $\displaystyle=\Psi(v_{j}),$ (13.3a) $\displaystyle
y_{j+1}$ $\displaystyle=Hv_{j+1}+\eta_{j+1}.$ (13.3b)

We note that

$Hv=w,\quad H^{\perp}v=u.$

######  Remark 13.1.

Typically either of the following are used to construct artificial data
$\\{y_{j}\\}$ for the filtering algorithm, given a single instance of data
$y$:

$y_{j+1}=\begin{cases}&y\text{ (unperturbed obervations).}\\\
&y+\bar{\eta}_{j+1},\bar{\eta}_{j+1}\sim\mathcal{N}(0,\Gamma)\text{ (perturbed
observations).}\end{cases}$

The first choice is natural if viewing the algorithm as a sequential
optimizer; the latter, in the linear case, is natural when seeking to draw
samples from the posterior.∎

We now apply the EnKF to the dynamics/data model (13.3). We obtain, for
$n=1,\dots,N$,

$\displaystyle\widehat{v}_{j+1}^{(n)}$ $\displaystyle=\Psi(v_{j}^{(n)}),$
(13.4) $\displaystyle\bar{v}_{j+1}$
$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}\widehat{v}_{j+1}^{(n)},$ (13.5)
$\displaystyle\widehat{C}_{j+1}$
$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}\bigl{(}\widehat{v}_{j+1}^{(n)}-\bar{v}_{j+1}\bigr{)}\otimes\bigl{(}\widehat{v}_{j+1}^{(n)}-\bar{v}_{j+1}\bigr{)},$
(13.6) $\displaystyle v_{j+1}^{(n)}$
$\displaystyle=(I-K_{j+1}H)\widehat{v}_{j+1}^{(n)}+K_{j+1}y_{j+1}^{(n)},$
(13.7)

with the Kalman gain

$\displaystyle K_{j+1}$ $\displaystyle=\widehat{C}_{j+1}H^{T}S_{j+1},$
$\displaystyle S_{j+1}$ $\displaystyle=(HC_{j+1}H^{T}+\Gamma)^{-1}.$

Now we may simplify these expressions by using the specific $\Psi$, $v$, $H$
arising in the inverse problem. Writing

$\displaystyle\widehat{C}_{j+1}=\begin{bmatrix}C_{j+1}^{uu}&C_{j+1}^{uw}\\\
(C_{j+1}^{uw})^{T}&C_{j+1}^{ww}\end{bmatrix},\quad\quad\;\;\bar{v}_{j+1}=\begin{pmatrix}\bar{u}_{j+1}\\\
\bar{w}_{j+1}\end{pmatrix},$

we have

$\bar{u}_{j+1}=\frac{1}{N}\sum_{n=1}^{N}u_{j}^{(n)},\;\quad\quad\bar{w}_{j+1}=\frac{1}{N}\sum_{n=1}^{N}G(u_{j}^{(n)}):=\bar{G}_{j},$

and

$\displaystyle C_{j+1}^{uw}$
$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}\bigl{(}u_{j}^{(n)}-\bar{u}_{j+1}\bigr{)}\otimes\bigl{(}G(u_{j}^{(n)})-\bar{G}_{j}\bigr{)},$
$\displaystyle C_{j+1}^{ww}$
$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}\bigl{(}G(u_{j}^{(n)})-\bar{G}_{j}\bigr{)}\otimes\bigl{(}G(u_{j}^{(n)})-\bar{G}_{j}\bigr{)}.$

There is a similar expression for $C_{j+1}^{uu}$, but as we will show in what
follows it is not needed for the unknown parameter $u$ update formula. Noting
that, because of the structure of $H$, $S_{j+1}=(C_{j+1}^{ww}+\Gamma)^{-1}$ we
obtain

$\displaystyle
K_{j+1}=\begin{pmatrix}C_{j+1}^{uw}(C_{j+1}^{ww}+\Gamma)^{-1}\\\
(C_{j+1}^{ww})(C_{j+1}^{ww}+\Gamma)^{-1}\end{pmatrix}.$ (13.8)

Combining equation (13.8) with the update equation within (13.4) it follows
that

$\displaystyle\bigl{\\{}v_{j}^{(n)}\bigr{\\}}_{n=1}^{N}\rightarrow\bigl{\\{}v_{j+1}^{(n)}\bigr{\\}}_{n=1}^{N}$

and

$\displaystyle\bigl{\\{}H^{\perp}v_{j}^{(n)}\bigr{\\}}_{n=1}^{N}\rightarrow\bigl{\\{}H^{\perp}v_{j+1}^{(n)}\bigr{\\}}_{n=1}^{N},$

and hence that

$\displaystyle
u_{j+1}^{(n)}=H^{\perp}v_{j+1}^{(n)}=u_{j}^{(n)}+C_{j+1}^{uw}(C_{j+1}^{ww}+\Gamma)^{-1}\Bigl{(}y_{j+1}^{(n)}-G(u_{j}^{(n)})\Bigr{)}.$

Thus we have derived the EKI step

$u_{j+1}^{(n)}=u_{j}^{(n)}+C_{j+1}^{uw}(C_{j+1}^{ww}+\Gamma)^{-1}\Bigl{(}y_{j+1}^{(n)}-G(u_{j}^{(n)})\Bigr{)}.$
(13.9)

The full algorithm is described below:

Algorithm 5 Algorithm for Ensemble Kalman Inversion

1: Input: Initial distribution $\pi_{0}=\rho$, observations $Y_{J}$, number of
particles $N.$

2: Initial Sampling: Draw $N$ particles $u_{0}^{(n)}\sim\rho$ so that
$\pi_{0}^{N}=S^{N}\rho.$

3: Subsequent Sampling For $j=0,1,\dots,J-1$, perform

  1. 1.

Set $\bar{u}_{j+1}=\frac{1}{N}\sum_{n=1}^{N}u_{j}^{(n)}.$

  2. 2.

Set $\bar{G}_{j}=\frac{1}{N}\sum_{n=1}^{N}G(u_{j}^{(n)}).$

  3. 3.

Set
$C_{j+1}^{uw}=\frac{1}{N}\sum_{n=1}^{N}\bigl{(}u_{j}^{(n)}-\bar{u}_{j+1}\bigr{)}\otimes\bigl{(}G(u_{j}^{(n)})-\bar{G}_{j}\bigr{)}.$

  4. 4.

Set
$C_{j+1}^{ww}=\frac{1}{N}\sum_{n=1}^{N}\bigl{(}G(u_{j}^{(n)})-\bar{G}_{j}\bigr{)}\otimes\bigl{(}G(u_{j}^{(n)})-\bar{G}_{j}\bigr{)}.$

  5. 5.

$u_{j+1}^{(n)}=u_{j}^{(n)}+C_{j+1}^{uw}\bigl{(}C_{j+1}^{ww}+\Gamma)^{-1}(y_{j+1}^{(n)}-G(u_{j}^{(n)})\bigr{)}.$

4: Output: $N$ particles $u_{J}^{(1)},u_{J}^{(2)},\dots,u_{J}^{(N)}.$

######  Remark 13.2.

This algorithm may be viewed as a derivative-free optimization method, within
the broad class that contains genetic, swarm or ant colony optimization. ∎

The algorithm has the following invariant subspace property.

######  Theorem 13.3.

(Space of Ensemble) Define
$\mathcal{A}={\text{Span}}\bigl{\\{}u^{(n)}_{0}:n=1,\ldots,N\bigr{\\}}$, then
for all $0\leq j\leq J$ and all $1\leq n\leq N$, $u^{(n)}_{j}$ defined by the
iteration (13.9) lie in $\mathcal{A}$.

###### Proof.

The proof proceeds by induction. We let
$\mathfrak{u}_{j}\in({\mathbb{R}^{d}})^{N}$ denote the collection of
$\\{u^{(n)}_{j}\\}_{n=1}^{N}.$ At first, let $0\leq j\leq J$ and let $1\leq
m\leq N$ and define the following two quantities:

$\displaystyle d^{(n)}_{j}$
$\displaystyle:=(C^{ww}_{j+1}+\Gamma)^{-1}\Bigl{(}y^{(n)}_{j+1}-G(u^{(n)}_{j})\Bigr{)},$
$\displaystyle D_{mn}(\mathfrak{u}_{j})$ $\displaystyle:=-\langle
d^{(n)}_{j},G(u^{(m)}_{j})-\bar{G}_{j}\rangle.$

Notice that the $d^{(n)}_{j}$ are elements of the data space
${\mathbb{R}}^{J}$, whereas the $D_{mn}(\mathfrak{u}_{j})$ are scalar
quantities. Furthermore notice that

$\displaystyle u_{j+1}^{(n)}$
$\displaystyle=u_{j}^{(n)}+C_{j+1}^{uw}(C_{j+1}^{ww}+\Gamma)^{-1}\Bigl{(}y_{j+1}^{(n)}-G(u_{j}^{(n)})\Bigr{)}$
$\displaystyle=u_{j}^{(n)}+C_{j+1}^{uw}d^{(n)}_{j}$
$\displaystyle=u_{j}^{(n)}+\frac{1}{N}\sum_{m=1}^{N}\bigl{(}u_{j}^{(m)}-\bar{u}_{j+1}\bigr{)}\otimes\bigl{(}G(u_{j}^{(m)})-\bar{G}_{j}\bigr{)}d^{(n)}_{j}$
$\displaystyle=u^{(n)}_{j}-\frac{1}{N}\sum\limits_{m=1}^{N}D_{mn}(\mathfrak{u}_{j})(u^{(m)}_{j}-\bar{u}_{j+1}).$

From the definition of $\bar{G}_{j}$ and the bilinearity of the inner product
it follows that $\sum\limits_{m=1}^{N}D_{mn}(\mathfrak{u}_{j})=0$. Therefore
the update expression may be rewritten, for $0\leq j\leq J$ and $1\leq n\leq
N,$ as

$u^{(n)}_{j+1}=u^{(n)}_{j}-\frac{1}{N}\sum\limits_{m=1}^{N}D_{mn}(\mathfrak{u}_{j})u^{(m)}_{j}.$

Hence if the property holds for all the particles of the time step $j$, it
will clearly be the case for all the particles at time step $j+1$. ∎

######  Remark 13.4.

The choice of the initial ensemble of particles $\\{u^{(n)}_{0}\\}_{n=1}^{N}$
is thus crucial to the performance of the EnKF, since the algorithm remains in
the initial space of the initial ensemble. In the setting of Bayesian inverse
problems the initial ensemble is frequently created by drawing from the prior.
Alternatively, if the prior is Gaussian, the first $N$ covariance eigenvectors
may be used, ordered by decreasing variance contribution. More generally any
truncated basis for the space $\mathbb{R}^{d}$ is a natural initial ensemble.
However the question of how to adaptively learn a good choice of ensemble
subspace, in response to observed data, is unexplored and potentially
fruitful. ∎

######  Remark 13.5.

We describe an alternative way to approach the derivation of the EKI update
formulae. We apply Theorem 10.1 with the specific structure arising from the
dynamical system used in EKI. To this end we define

${\mathsf{I}}_{n}(b):=\frac{1}{2}\Big{\|}y_{j+1}^{(n)}-G({u}_{j}^{(n)})-\frac{1}{N}\sum_{n=1}^{N}b_{n}\bigl{(}G({u}_{j}^{(n)})-\bar{G}_{j}\bigr{)}\Big{\|}^{2}_{\Gamma}+\frac{1}{2N}\sum_{n=1}^{N}b_{n}^{2}.$
(13.10)

Once this quadratic form has been minimized with respect to $b$ then the upate
formula (10.18) gives

$\displaystyle u_{j+1}^{(n)}$
$\displaystyle={u}_{j}^{(n)}+\frac{1}{N}\sum_{n=1}^{N}b_{n}\bigl{(}{u}^{(n)}_{j}-\bar{u}_{j}\bigr{)},$
$\displaystyle w_{j+1}^{(n)}$
$\displaystyle=G({u}_{j}^{(n)})+\frac{1}{N}\sum_{n=1}^{N}b_{n}\bigl{(}G({u}^{(n)}_{j})-\bar{G}_{j}\bigr{)}.$

(Note that the vector $b$ depends on $n$; we have supressed this dependence
for notational convenience.) Theorem 13.3 is an immediate consequence of this
structure. ∎

####  13.3 Linking Ensemble Kalman Inversion and SMC

Here we link the two preceding sections. In the first subsection we describe
how EKI may be linked with SMC through a particular scaling of $\Gamma$ with
$h$. We then take the limit $h\to 0$ and obtain a system or ordinary
differential equations. In the second subsection we study the resulting
algorithm in the case in which $G$ is linear, leading to insight into the EKI
algorithm in the iterated context.

#####  13.3.1 Continuous Time Limit

Although we have emphasized the optimization perspective on ensemble methods,
we may think of one-step of ensemble Kalman inversion as approximating the
filtering mapping
$\pi_{j}=\operatorname{\mathbb{P}}(u_{j}|Y_{j})\mapsto\pi_{j+1}=\operatorname{\mathbb{P}}(u_{j+1}|Y_{j+1});$
this will be a good approximation when the measures are close to Gaussian and
the number of ensemble members $N$ is large. In order to link this to the
mapping $\pi_{j+1}={\mathcal{A}}\pi_{j}$ in the presentation of SMC in the
first section of this chapter, we set $\Gamma=\frac{1}{h}\Gamma_{0},$ and we
assume unperturbed observations $y_{j}^{(n)}=y$ for all $j.$ Let
${\mathfrak{u}}_{j}\in(\mathbb{R}^{d})^{N}$ denote, as in the previous
section, the collection of particles $\\{u^{(n)}_{j}\\}_{n=1}^{N}$ at discrete
time $j.$ Now define

$D^{0}_{mn}(\mathfrak{u}_{j}):=\Bigl{\langle}G(u^{(n)}_{j})-y,G(u^{(m)}_{j})-\bar{G}_{j}\Bigr{\rangle}_{\Gamma_{0}}.$

Note that then, to leading order in $h\ll 1$,

$\displaystyle D_{mn}(\mathfrak{u}_{j})$
$\displaystyle=-\Bigl{\langle}(C^{ww}_{j+1}+\Gamma)^{-1}(y-G(u^{(n)}_{j})),G(u^{(m)}_{j})-\bar{G}_{j}\Bigr{\rangle}$
$\displaystyle\approx h\langle
G(u_{j}^{(n)})-y,G(u^{(m)})-\bar{G}_{j}\rangle_{\Gamma_{0}}$
$\displaystyle=hD^{0}_{mn}(\mathfrak{u}_{j}).$

It follows that, also to leading order in $h$,

$u^{(n)}_{j+1}\approx
u^{(n)}_{j}-\frac{h}{N}\sum\limits_{m=1}^{N}D^{0}_{mn}(\mathfrak{u}_{j})u^{(m)}_{j}.$

Then letting $h\to 0$ yields

$\frac{du^{(n)}}{dt}=-\frac{1}{N}\sum\limits_{m=1}^{N}D^{0}_{mn}(\mathfrak{u})u^{(m)},$
(13.11)

where, similarly as above, $\mathfrak{u}$ represents the collection of $N$
particles.

######  Remark 13.6.

Notice that equation (13.11) has families of fixed points where: (i) either
the particles fit the data exactly (which corresponds to the left-hand side in
the inner product defining $D^{0}_{mn}(u)$ being zero); or (ii) all particles
collapse on their mean value (which corresponds to the right-hand side of the
same inner product). This suggests that the system of ordinary differential
equations which describes the behaviour of ensemble Kalman inversion is driven
by two desirable attributes: matching the data and achieving consensus. ∎

#####  13.3.2 Linear Setting

Now suppose $G(\cdot)$ is a linear map denoted $A\cdot$, and define

$\bar{u}=\frac{1}{N}\sum\limits_{n=1}^{N}u^{(n)}.$

In this setting we have

$\displaystyle\frac{du^{(n)}}{dt}$
$\displaystyle=-\frac{1}{N}\sum\limits_{m=1}^{N}\bigl{\langle}Au^{(n)}-y,A(u^{(m)}-\bar{u})\bigr{\rangle}_{\Gamma_{0}}u^{(m)}$
$\displaystyle=-\frac{1}{N}\sum\limits_{m=1}^{N}\bigl{\langle}Au^{(n)}-y,A(u^{(m)}-\bar{u})\bigr{\rangle}_{\Gamma_{0}}(u^{(m)}-\bar{u})$
$\displaystyle=-\frac{1}{N}\sum\limits_{m=1}^{N}\bigl{\langle}A^{T}\Gamma_{0}^{-1}(Au^{(n)}-y),u^{(m)}-\bar{u}\bigr{\rangle}(u^{(m)}-\bar{u})$
$\displaystyle=-\frac{1}{N}\sum\limits_{m=1}^{N}\bigl{(}u^{(m)}-\bar{u}\bigr{)}\otimes\bigl{(}u^{(m)}-\bar{u}\bigr{)}\bigl{(}A^{T}\Gamma_{0}^{-1}(Au^{(n)}-y)\bigr{)}.$

If we define

$C(\mathfrak{u})=\frac{1}{N}\sum\limits_{n=1}^{N}\bigl{(}u^{(n)}-\bar{u}\bigr{)}\otimes\bigl{(}u^{(n)}-\bar{u}\bigr{)}$
(13.12)

then we find that, for $n=1,\dots,N$,

$\displaystyle\frac{du^{(n)}}{dt}$
$\displaystyle=-C(\mathfrak{u})A^{T}\Gamma_{0}^{-1}(Au^{(n)}-y)$
$\displaystyle=-C(\mathfrak{u})\nabla\mathsf{L}_{0}(u^{(n)}),$

where

$\mathsf{L}_{0}(u)=\frac{1}{2}|y-Au|_{\Gamma_{0}}^{2}.$

This corresponds to a gradient descent in the subspace defined by the initial
ensemble. In particular

$\displaystyle\frac{d}{dt}\mathsf{L}\bigl{(}u^{(n)}(t)\bigr{)}$
$\displaystyle=\Bigl{\langle}\nabla\mathsf{L}\bigl{(}u^{(n)}(t)\bigr{)},\frac{du^{(n)}}{dt}(t)\Bigr{\rangle}$
$\displaystyle=-\Bigl{|}C(u)^{1/2}\nabla\mathsf{L}\bigl{(}u^{(n)}(t)\bigr{)}\Bigr{|}^{2}\leq
0,$

demonstrating that the loss function $\mathsf{L}$ is decreasing along the
trajectory associated to each ensemble member.

We note also that, if we define $e^{(n)}=u^{(n)}-\bar{u}$, then

$\frac{de^{(n)}}{dt}=-C(\mathfrak{u})A^{T}\Gamma_{0}^{-1}Ae^{(n)}.$

Because

$C:=C(\mathfrak{u})=\frac{1}{N}\sum_{n=1}^{N}e^{(n)}\otimes e^{(n)}$

it follows that

$\frac{dC}{dt}=-2CA^{T}\Gamma_{0}^{-1}AC.$

Note further that if $C$ is invertible (requires $N\geq d$) then the inverse
$P$ satisfies

$\frac{dP}{dt}=2A^{T}\Gamma_{0}^{-1}A$

so that $\|P\|\to\infty$ as $t\to\infty.$ This fact, which can be suitably
interpreted when $C$ is only invertible on a subspace, suggests that the
covariance shrinks with time, causing ensemble collapse – consensus – while at
the same time driving minimization of the loss $\mathsf{L}$ over an
appropriate subspace.

####  13.4 Discussion and Bibliography

The idea of using particle filters to sample general distributions, including
those arising in Bayesian inversion, maybe be found in [26]; a recent
application to a Bayesian inverse problem, which demonstrated the potential of
the methodology in that context, is [67]. A simple proof of convergence of the
method may be found in [11]; it is based on the proof described in [97] for
the standard bootstrap particle filter.

The use of the ensemble Kalman filter for parameter estimation was introduced
in the papers [79, 6] in which a physical dynamical model was appended with
trivial dynamics for the parameters in order to estimate them; the idea was
extended to learn an entire field of parameter values in [88]. The paper [104]
was the first to do what we do here, namely to consider all the data at once
and a single mapping of unknown parameters to data. In that paper only one
iteration is used; the papers [19, 30] demonstrated how iteration could be
useful. See also the book [91] for the use of ensemble Kalman methods in oil
reservoir simulation.

Development of the method for general inverse problems was undertaken in [58],
and further development of iterative methods is described in [57, 60].
Ensemble inversion was analyzed in [103], including the continuous-time limit
and gradient flow structure described here. The potential for ensemble
inversion in the context of hierarchical Bayesian methods is demonstrated in
[18]. The idea that we explain here, of obtaining an evolution equation for
the covariance which is satisfied exactly by ensemble methods, appears in the
remarkable papers [98, 10], and in the other papers of Bergemann and Reich
referenced therein. Their work is at the heart of the analysis in [103] which
demonstrates the ensemble collapse and approximation properties of the
ensemble method when applied to linear inverse problems.

### References

  * Abarbanel [2013] H. Abarbanel.  _Predicting The Future: Completing Models Of Observed Complex Systems_.  Springer, 2013.
  * Agapiou et al. [2013] S. Agapiou, S. Larsson, and A. M. Stuart.  Posterior contraction rates for the Bayesian approach to linear ill-posed inverse problems.  _Stochastic Processes and their Applications_ , 123(10):3828–3860, 2013.
  * Agapiou et al. [2017a] S. Agapiou, M. Burger, M. Dashti, and T. Helin.  Sparsity-promoting and edge-preserving maximum a posteriori estimators in non-parametric Bayesian inverse problems.  _Inverse Problems_ , 34(4), 2017a.
  * Agapiou et al. [2017b] S. Agapiou, O. Papaspiliopoulos, D. Sanz-Alonso, and A. M. Stuart.  Importance sampling: Intrinsic dimension and computational cost.  _Statistical Science_ , 32(3):405–431, 2017b.
  * Anderson and Moore [1979] B. Anderson and J. Moore.  Optimal filtering.  _Englewood Cliffs_ , 21:22–95, 1979.
  * Anderson [2001] C. Anderson.  An ensemble adjustment Kalman filter for data assimilation.  _Monthly Weather Review_ , 129:2284–2903, 2001.
  * Anderson [2014] C. Anderson.  Monte Carlo methods and importance sampling.  _Lecture Notes for Statistical Genetics_ , 2014.
  * Bassiri et al. [2016] P. Bassiri, C. Holmes, and S. Walker.  A general framework for updating belief distributions.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology_ , 78(5):1103–1130, 2016.
  * Bell [2001] B. Bell.  The iterated Kalman smoother as a Gauss–Newton method.  _SIAM Journal on Optimization_ , 4(3):626–636., 2001.
  * Bergemann and Reich [2012] K. Bergemann and S. Reich.  An ensemble kalman-bucy filter for continuous data assimilation.  _Meteorologische Zeitschrift_ , 127(5):1417–1440, 2012.
  * Beskos et al. [1994] A. Beskos, A. Jasra, K. Law, R. Tempone, and Y. Zhou.  Multilevel Sequential Monte Carlo Samplers.  _Stochastic Processes and their Applications_ , 127(5):1417–1440, 1994.
  * Bickel et al. [1994] P. Bickel, B. Li, and T. Bengtsson.  Pushing the limits of contemporary statistics: Contributions in honor of jayanta k. ghosh: Sharp failure rates for the bootstrap particle filter in high dimensions.  _Institute of Mathematical Statistics_ , pages 318–329, 1994.
  * Bishop [2006] C. Bishop.  _Pattern Recognition and Machine Learning_.  Springer, 2006.
  * Brett et al. [2013] C. Brett, K. Lam, K. Law, D. McCormick, M. Scott, and A. M. Stuart.  Accuracy and stability of filters for dissipative PDE’s.  _Physica D: Nonlinear Phenomena_ , 2013.
  * Bröcker [2013] J. Bröcker.  Existence and uniqueness for four-dimensional variational data assimilation in discrete time.  _SIAM Journal on Applied Dynamical Systems_ , 16(1):361–374, 2013.
  * Brooks et al. [2011] S. Brooks, A. Gelman, G. Jones, and X. Meng.  Handbook of Markov chain Monte Carlo.  _CRC press_ , 2011.
  * Carrassi et al. [2018] A. Carrassi, M. Bocquet, L. Bertino, and G. Geir.  Data assimilation in the geosciences: An overview of methods, issues, and perspectives.  _Wiley Interdisciplinary Reviews: Climate Change_ , 9(5), 2018.
  * Chada et al. [2017] N. Chada, M. Iglesias, L. Roininen, and A. M. Stuart.  Parameterizations for ensemble Kalman inversion.  _Inverse Problems_ , 34(2018), 2017.
  * Chen and Oliver [2002] Y. Chen and D. Oliver.  Ensemble randomized maximum likelihood method as an iterative ensemble smoother.  _Mathematical Geosciences_ , 44(1):1–26, 2002\.
  * Cotter et al. [2010a] S. Cotter, M. Dashti, and A. M. Stuart.  Approximation of bayesian inverse problems for pde’s.  _SIAM Journal on Numerical Analysis_ , 48(1):322–345, 2010a.
  * Cotter et al. [2010b] S. Cotter, M. Dashti, and A. M. Stuart.  MCMC methods for functions: modifying old algorithms to make them faster.  _SIAM Journal on Numerical Analysis_ , 48(1):322–345, 2010b.
  * Crisan et al. [1998] D. Crisan, P. Moral, and T. Lyons.  Discrete filtering using branching and interacting particle systems.  _Université de Toulouse. Laboratoire de Statistique et Probabilités [LSP]_ , 1998.
  * Dashti and Stuart [2017] M. Dashti and A. M. Stuart.  Bayesian approach to inverse problems.  _Handbook of Uncertainty Quantification_ , pages 311–428, 2017.
  * Dashti et al. [2013] M. Dashti, K. Law, A. M. Stuart, and J. Voss.  Map estimators and their consistency in Bayesian nonparametric inverse problems.  _Inverse Problems_ , 29(9), 2013.
  * del Moral [2004] P. del Moral.  Feynman-Kac Formulae.  _Genealogical and Interacting Particle Systems with Applications_ , pages 47–93, 2004.
  * del Moral [2006] P. del Moral.  Sequential Monte Carlo samplers.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 68(3):411–436, 2006.
  * der Vaart [1998] A. V. der Vaart.  _Asymptotic Statistics_.  Cambridge University Press, 1998.
  * Doucet et al. [2000] A. Doucet, S. Godsill, and C. Andrieu.  On sequential Monte Carlo sampling methods for Bayesian filtering.  _Statistics and Computing_ , 10(3), 2000.
  * Doucet et al. [2001] A. Doucet, N. Freitas, and N. Gordon.  _An introduction to sequential Monte Carlo methods_.  Springer, 2001.
  * Emerick and Reynolds [2013] A. Emerick and A. Reynolds.  Investigation of the sampling performance of ensemble-based methods with a simple reservoir model.  _Computational Geosciences_ , 17(2):325–350, 2013.
  * England et al. [1996] H. England, M. Hanke, and A. Neubauer.  _Regularization of inverse problems_.  Springer Science and Business Media, 1996.
  * Ernst et al. [2015] O. Ernst, B. Sprungk, and H. Starkloff.  Analysis of the ensemble and polynomial chaos Kalman filters in Bayesian inverse problems.  _SIAM/ASA Journal on Uncertainty Quantification_ , 3(1):823–851, 2015.
  * Evans [1995] G. Evans.  Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics.  _Journal of Geophysical Research: Oceans_ , 99(c5):10143–10162, 1995.
  * Evans and Leeuwen [1995] G. Evans and P. V. Leeuwen.  Methods for ensemble prediction.  _Monthly Weather Review_ , 123(7):2181–2196, 1995.
  * Evans and Leeuwen [1996] G. Evans and P. V. Leeuwen.  Assimilation of Geosat altimeter data for the Agulhas current using the ensemble Kalman filter with a quasigeostrophic model.  _Monthly Weather Review_ , 124(1):85–96, 1996\.
  * Evans and Swartz [1995] M. Evans and T. Swartz.  Methods for approximating integrals in statistics with special emphasis on Bayesian integration problems.  _Statistical Science_ , pages 254–272, 1995.
  * Evensen [2009] G. Evensen.  _Data Assimilation: the Ensemble Kalman Filter_.  Springer Science and Business Media, 2009.
  * Fisher et al. [2009] M. Fisher, J. Nocedal, Y. Trémolet, and S. Wright.  Data assimilation in weather forecasting: a case study in PDE-constrained optimization.  _Optimization and Engineering_ , 10(3):409–426, 2009.
  * Franklin [1970] J. Franklin.  Well-posed stochastic extensions of ill-posed linear problems.  _Journal of mathematical analysis and applications_ , 31(3):682–716, 1970.
  * Gamerman and Lopes [2006] D. Gamerman and H. Lopes.  Markov chain Monte Carlo: stochastic simulation for Bayesian inference.  _CRC Press_ , 2006.
  * Garcia Trillos and Sanz-Alonso [2018a] N. Garcia Trillos and D. Sanz-Alonso.  Continuum limits of posteriors in graph bayesian inverse problems.  _SIAM Journal on Mathematical Analysis_ , 50(4):4020–4040, 2018a.
  * Garcia Trillos and Sanz-Alonso [2018b] N. Garcia Trillos and D. Sanz-Alonso.  The Bayesian update: variational formulations and gradient flows.  _Bayesian Analysis_ , 2018b.
  * Garcia Trillos et al. [2017] N. Garcia Trillos, Z. Kaplan, T. Samakhoana, and D. Sanz-Alonso.  On the consistency of graph-based Bayesian learning and the scalability of sampling algorithms.  _arXiv preprint arXiv:1710.07702_ , 2017.
  * Garcia Trillos et al. [2019] N. Garcia Trillos, Z. Kaplan, and D. Sanz-Alonso.  Variational characterizations of local entropy and heat regularization in deep learning.  _arXiv preprint arXiv:1901.10082_ , 2019.
  * Ghil et al. [1991] M. Ghil, S. Cohn, J. Tavantzis, and K. Bube.  Application of estimation theory to numerical weather prediction.  _Dynamic Meteorology: Data Assimilation Methods_ , 1991.
  * Gibbs and Su [2002] A. Gibbs and F. Su.  On choosing and bounding probability metrics.  _International Statistical Review_ , 70(3):419–435, 2002.
  * Giles [2015] M. Giles.  Multilevel Monte Carlo methods.  _Acta Numerica_ , 24:259–328, 2015.
  * Gine and Nickl [2015] E. Gine and R. Nickl.  _Mathematical foundations of infinite-dimensional statistical models_.  Cambridge University Press, 2015.
  * Gland et al. [2009] F. Gland, V. Monbet, and V. Tran.  Large sample asymptotics for the ensemble Kalman filter.  _PhD Thesis_ , 2009.
  * Gottwald and Majda [2013] G. Gottwald and A. Majda.  A mechanism for catastrophic filter divergence in data assimilation for sparse observation networks.  _Nonlinear Processes in Geophysic_ , 20(5):705–712, 2013.
  * Hairer et al. [2013] M. Hairer, A. M. Stuart, J. Voss, and P. Wiberg.  Analysis of SPDEs arising in path sampling. Part I: The Gaussian case.  _Communications in Mathematical Sciences_ , 3(4):587–603, 2013.
  * Hammersley and Handscomb [1964] J. Hammersley and D. Handscomb.  Percolation processes.  _Monte Carlo Methods_ , pages 134–141, 1964.
  * Harbey [1964] A. Harbey.  _Forecasting, structural time series models and the Kalman filter_.  Cambridge university press, 1964.
  * Hayden et al. [2011] K. Hayden, E. Olson, and E. Titi.  Discrete data assimilation in the Lorenz and 2D Navier–Stokes equations.  _Physica D: Nonlinear Phenomena_ , 240(18):1416–1425, 2011.
  * Helin and Burger [2015] T. Helin and M. Burger.  Maximum a posteriori probability estimates in infinite-dimensional bayesian inverse problems.  _Inverse Problems_ , 31(8), 2015.
  * Houtekamer and Mitchell [1998] P. Houtekamer and H. Mitchell.  Data assimilation using an ensemble Kalman filter technique.  _Monthly Weather Review_ , 126(3):796–811, 1998\.
  * Iglesias [2015] M. Iglesias.  Iterative regularization for ensemble data assimilation in reservoir models.  _Computational Geosciences_ , 19(1):177–212, 2015.
  * Iglesias et al. [2014a] M. Iglesias, K. Law, and A. M. Stuart.  Ensemble Kalman methods for inverse problems.  _Inverse Problems_ , 29(4):134–141, 2014a.
  * Iglesias et al. [2014b] M. Iglesias, K. Lin, and A. M. Stuart.  Well-posed Bayesian geometric inverse problems arising in subsurface flow.  _Inverse Problems_ , 30(11), 2014b.
  * Iglesias et al. [2016] M. Iglesias, K. Law, and A. M. Stuart.  A regularizing iterative ensemble Kalman method for PDE-constrained inverse problems.  _Inverse Problems_ , 32(2), 2016.
  * Jazwinski [2007] A. Jazwinski.  _Stochastic Processes and Filtering Theory_.  Courier Corporation, 2007.
  * Johansen and Doucet [2008] A. Johansen and A. Doucet.  A note on auxiliary particle filters.  _Statistics and Probability Letters_ , 78(12):1498–1504, 2008.
  * Kaipo and Somersalo [2006] J. Kaipo and E. Somersalo.  Statistical and Computational Inverse Problems.  _Springer Science & Business Media_, 160, 2006.
  * Kalman [1960] R. Kalman.  A new approach to linear filtering and prediction problems.  _Journal of Basic Engineering_ , 82(1):35–45, 1960.
  * Kalman and Bucy [1961] R. Kalman and R. Bucy.  New results in linear filtering and prediction theory.  _Journal of Basic Engineering_ , 83(1):95–108, 1961.
  * Kalnay [2003] E. Kalnay.  _Atmospheric modeling, data assimilation and predictability_.  Cambridge University Press, 2003.
  * Kantas et al. [2014] N. Kantas, A. Beskos, and A. Jasra.  Sequential Monte Carlo methods for high-dimensional inverse problems: a case study for the Navier Stokes equations.  _SIAM Journal on Uncertainty Quantification_ , 2(1):464–489, 2014.
  * Kawai [2017] R. Kawai.  Adaptive importance sampling Monte Carlo simulation for general multivariate probability laws.  _Journal of Computational and Applied Mathematics_ , 319:440–459, 2017.
  * Kelly and Stuart [2014a] D. Kelly and A. M. Stuart.  Well-posedness and accuracy of the ensemble Kalman filter in discrete and continuous time.  _Nonlinearity_ , 27(10), 2014a.
  * Kelly and Stuart [2014b] D. Kelly and A. M. Stuart.  Ergodicity and accuracy of optimal particle filters for Bayesian data assimilation.  _arXiv preprint arXiv:1611.08761_ , 2014b.
  * Knapik et al. [2011] B. Knapik, A. van der Vaart, and J. van Zanten.  Bayesian inverse problems with Gaussian priors.  _Annals of Statistics_ , 39(5):2626–2657, 2011\.
  * Law et al. [2012] K. Law, A. Shukla, and A. M. Stuart.  Analysis of the 3DVAR filter for the partially observed Lorenz’63 model.  _arXiv preprint arXiv:1212.4923_ , 2012.
  * Law et al. [2015] K. Law, A. M. Stuart, and K. Zygalakis.  _Data assimilation_.  Springer, 2015.
  * Law et al. [2016] K. Law, D. Sanz-Alonso, A. Shukla, and A. M. Stuart.  Filter accuracy for the Lorenz 96 model: Fixed versus adaptive observation operators.  _Physica D: Nonlinear Phenomena_ , 325:1–13, 2016.
  * Leeuwen et al. [2015] P. V. Leeuwen, Y. Cheng, and S. Reich.  _Nonlinear data assimilation_.  Springer, 2015.
  * Lieberman et al. [2010] C. Lieberman, C. Willcox, and O. Ghattas.  Parameter and state model reduction for large-scale statistical inverse problems.  _SIAM Journal on Scientific Computing_ , 32(5):2535–2542, 2010.
  * Lindvall [2002] T. Lindvall.  _Lectures on the Coupling Method_.  Springer, 2002.
  * Lorenc [1986] A. Lorenc.  Analysis methods for numerical weather prediction.  _Quarterly Journal of the Royal Meteorological Society_ , 112(474):1177–1194, 1986.
  * Lorentzen et al. [2001] R. Lorentzen, R. Fjelde, J. FrØyen, A. Lage, G. Naevdal, and E. Vefring.  Underbalanced and low-head drilling operations: Real time interpretation of measured data and operational support.  _SPE Annual Technical Conference and Exhibition_ , 2001.
  * Lu et al. [2017] Y. Lu, A. M. Stuart, and H. Weber.  Gaussian Approximations for Probability Measures on $\mathbb{R}^{d}$.  _SIAM/ASA Journal on Uncertainty Quantification_ , 5(1):1136–1165, 2017.
  * MacKay [2003] D. MacKay.  _Information Theory, Inference and Learning Algorithms_.  Cambridge University Press, 2003.
  * Majda and Harlim [2012] A. Majda and J. Harlim.  _Filtering Complex Turbulent Systems_.  Cambridge University Press, 2012.
  * Martin et al. [2012] J. Martin, L. Wilcox, C. Burstedde, and G. Omar.  A stochastic Newton MCMC method for large-scale statistical inverse problems with application to seismic inversion.  _SIAM Journal on Scientific Computing_ , 34(3):A1460–A1487, 2012.
  * Marzouk and Xiu [2009] Y. Marzouk and D. Xiu.  A stochastic collocation approach to Bayesian inference in inverse problems.  _Communications in Computational Physics_ , 6(4):826–847, 2009.
  * Mattingly et al. [2002] J. Mattingly, A. Stuart, and D. Higham.  Ergodicity for PDE’s and approximations: locally Lipschitz vector fields and degenerate noise.  _Stochastic Processes and Their Applications_ , 101(2):185–232, 2002.
  * Meyn and Tweedie [2012] S. Meyn and R. Tweedie.  _Markov Chains and Stochastic Stability_.  Springer Science and Business Media, 2012.
  * Moodey et al. [2013] A. Moodey, A. Lawless, R. Potthast, and P. V. Leeuwen.  Nonlinear error dynamics for cycled data assimilation methods.  _Inverse Problems_ , 29(2), 2013.
  * Naevdal et al. [2002] G. Naevdal, T. Mannseth, and E. Vefring.  Near-well reservoir monitoring through ensemble Kalman filter.  _Proceeding of SPE Improved Oil Recovery Symposium_ , 2002.
  * Nickl [2017] R. Nickl.  Bernstein-von mises theorems for statistical inverse problems : Schrö odinger equation.  _arXiv preprint arXiv:1707.01764)_ , 2017.
  * Nielsen and Garcia [2009] F. Nielsen and V. Garcia.  Statistical exponential families: A digest with flash cards.  _arXiv preprint arXiv:0911.4863)_ , 2009.
  * Oliver et al. [2008] D. Oliver, A. Reynolds, and N. Liu.  _Inverse theory for petroleum reservoir characterization and history matching_.  Cambridge University Press), 2008.
  * Petersen and Pedersen [2008] K. Petersen and M. Pedersen.  _The matrix cookbook_.  Technical University of Denmark), 2008.
  * Pinski et al. [2015a] F. Pinski, F. Simpson, A. M. Stuart, and H. Weber.  Algorithms for Kullback–Leibler approximation of probability measures in infinite dimensions.  _SIAM Journal on Scientific Computing_ , 37(6):A2733–A2757, 2015a.
  * Pinski et al. [2015b] F. Pinski, F. Simpson, A. M. Stuart, and H. Weber.  Kullback–Leibler approximation for probability measures on infinite dimensional spaces.  _SIAM Journal on Mathematical Analysis_ , 47(6):4091–4122, 2015b.
  * Pitt and Shephard [1999] M. Pitt and N. Shephard.  Filtering via simulation: Auxiliary particle filters.  _Journal of the American Statistical Association_ , 94(446):590–599, 1999.
  * Rauch et al. [1965] H. Rauch, C. Striebel, and F. Tung.  Maximum likelihood estimates of linear dynamic systems.  _AIAA journal_ , 3(8):1445–1450, 1965.
  * Rebeschini and Handel [2015] P. Rebeschini and R. V. Handel.  Can local particle filters beat the curse of dimensionality?  _Annals of Applied Probability_ , 25(5):2809–2866, 2015.
  * Reich [2017] S. Reich.  A dynamical systems framework for intermittent data assimilation.  _IBIT Numerical Mathematics_ , 51(1):235–249, 2017.
  * Reich and Cotter [2015] S. Reich and C. Cotter.  _Probabilistic Forecasting and Bayesian Data Assimilation_.  Cambridge University Press), 2015.
  * Sanz-Alonso [2018] D. Sanz-Alonso.  Importance sampling and necessary sample size: An information theory approach.  _SIAM/ASA Journal on Uncertainty Quantification_ , 6(2):867–879, 2018.
  * Sanz-Alonso and Stuart [2015] D. Sanz-Alonso and A. M. Stuart.  Long-time asymptotics of the filtering distribution for partially observed chaotic dynamical systems.  _SIAM/ASA Journal on Uncertainty Quantification_ , 3(1):1200–1220, 2015.
  * Sanz-Alonso and Stuart [2017] D. Sanz-Alonso and A. M. Stuart.  Gaussian approximations of small noise diffusions in Kullback-Leibler divergence.  _Communications in Mathematical Sciences_ , 15(7):2087–2097, 2017.
  * Schillings and Stuart [2017] C. Schillings and A. M. Stuart.  Analysis of the ensemble Kalman filter for inverse problems.  _SIAM Journal on Numerical Analysis_ , 55(3), 2017.
  * Skjervheim et al. [2017] J. Skjervheim, J. Evensen, S. Aanonsen, B. Ruud, and T. Johansen.  Incorporating 4d seismic data in reservoir simulation models using ensemble kalman filter.  _SPE Journal_ , 12(3):282–292, 2017.
  * Snyder [2010] C. Snyder.  Inverse problems: a Bayesian perspective.  _Acta Numerica_ , 19:451–559, 2010.
  * Snyder [2011] C. Snyder.  Particle filters, the optimal proposal and high-dimensional systems.  _Proceedings of the ECMWF Seminar on Data Assimilation for Atmosphere and Ocean_ , pages 1–10, 2011.
  * Snyder et al. [2016] C. Snyder, T. Bengtsson, P. Bickel, and J. Anderson.  Obstacles to high-dimensional particle filtering.  _Monthly Weather Review_ , 136(12):4629–4640, 2016.
  * Tarantola [2015a] A. Tarantola.  Inverse problem theory and methods for model parameter estimation.  _SIAM_ , 2015a.
  * Tarantola [2015b] A. Tarantola.  Towards adjoint-based inversion for rheological parameters in nonlinear viscous mantle flow.  _Physics of the Earth and Planetary Interiors_ , 234:23–34, 2015b.
  * Tokdar et al. [2010] S. Tokdar, S. Kass, and R. Kass.  Importance sampling: a review.  _Wiley Interdisciplinary Reviews: Computational Statistics_ , 2(1):54–60, 2010.
  * Tong et al. [2015] X. Tong, A. Majda, and D. Kelly.  Nonlinear stability of the ensemble Kalman filter with adaptive covariance inflation.  _Nonlinearity_ , 29(2):54–60, 2015.
  * Tong et al. [2016] X. T. Tong, A. J. Majda, and D. Kelly.  Nonlinear stability and ergodicity of ensemble based Kalman filters.  _Nonlinearity_ , 29(2):657, 2016.


# Elliptic flow of light nuclei in Au+Au collisions at $\sqrt{s_{NN}}$ = 14.6,
19.6, 27, and 54.4 GeV using the STAR detector

Rishabh Sharma (for the STAR Collaboration)  Indian Institute of Science
Education and Research, Tirupati 517507, India
<EMAIL_ADDRESS>
###### Abstract

Loosely bound light nuclei are produced in abundance in heavy-ion collisions.
There are two main possible models to explain their production mechanism - the
thermal model and the coalescence model. The thermal model suggests that the
light nuclei are produced from a thermal source, where they are in equilibrium
with other species present in the fireball. However, due to the small binding
energies, the produced nuclei are not likely to survive the high-temperature
conditions of the fireball. The coalescence model tries to explain the
production of light nuclei by assuming that they are formed at later stages by
the coalescence of protons and neutrons near the kinetic freeze-out surface.
The final-state coalescence of nucleons will lead to the mass number scaling
of the elliptic flow ($v_{2}$) of light nuclei. This scaling states that the
$v_{2}$ of light nuclei scaled by their respective mass numbers will follow
very closely the $v_{2}$ of nucleons. Therefore, studying the $v_{2}$ of light
nuclei and comparing it with the $v_{2}$ of protons will help us in
understanding their production mechanism.

In this talk, we will present the transverse momentum ($p_{T}$) and centrality
dependence of $v_{2}$ of $d$, $t$, and ${}^{3}\text{He}$ in Au+Au collisions
at $\sqrt{s_{NN}}$ = 14.6, 19.6, 27, and 54.4 GeV. Mass number scaling of
$v_{2}(p_{T})$ of light (anti-)nuclei will be shown and physics implications
will be discussed.

##  1 Introduction

The study of light nuclei and their interaction in high-energy heavy-ion
collisions has been a subject of active theoretical and experimental
investigations [1]. The production mechanism of light nuclei in heavy-ion
collisions is not very well understood. There are two main models that
describe this mechanism: the thermal model and the coalescence model. The
thermal model suggests that light nuclei are formed near the chemical
freezeout (CFO) surface along with other hadrons [2]. However, the low binding
energies of light nuclei make it unlikely that they will be able to sustain
the high temperature at CFO. The coalescence model, on the other hand,
suggests that light nuclei might be formed by the coalescence of nucleons at
the later stages of evolution of the system [3]. This will result into the
mass number scaling whereby elliptic flow of light nuclei scaled by their
respective mass numbers follows closely to elliptic flow of protons [4].
Therefore, by examining the collective flow of light nuclei, valuable insights
can be gained into how they are produced in heavy-ion collisions.
In the following sections, we will report the elliptic flow ($v_{2}$) of d, t,
and 3He in Au+Au collisions at $\sqrt{s_{NN}}$ = 14.6, 19.6, 27, and 54.4 GeV.
We will also report the results of the centrality dependence study of $v_{2}$
of d in $\sqrt{s_{NN}}$ = 19.6, 27, and 54.4 GeV. Finally, we will show the
results from the mass number scaling study of $v_{2}(p_{T})$ of light nuclei.

##  2 Analysis details

The data presented in these proceedings was collected during the second phase
of the Beam Energy Scan (BES-II) program, conducted by the STAR experiment at
RHIC. Light nuclei identification was done using the Time Projection Chamber
(TPC) [5] and the Time of Flight (TOF) [6] detectors. TPC serves as the main
tracking detector in the STAR experiment and relies on the measurement of
specific ionization energy loss (dE/dx) within a large gas volume to identify
and track various charged particles. The TOF detector, on the other hand,
enables the identification of particles of interest by imposing a constraint
on their mass-square ($m^{2}$).
Elliptic flow, $v_{2}$, is the second order Fourier coefficient of the
azimuthal distribution of the produced nuclei relative to the reaction plane
of the Au+Au collision. Since it is not feasible to directly measure the
reaction plane angle in an experimental setup, we employ the TPC to construct
the second order event plane angle ($\Psi_{2}$) as a substitute for the
reaction plane angle [7]. In the next section, we will discuss the results of
$v_{2}$ of $d$, $t$, and ${}^{3}\text{He}$.

##  3 Results

###  3.1 Elliptic flow of light nuclei

Figure 1 shows $v_{2}$ as a function of $p_{T}$ in 0-80% centrality Au+Au
collisions at $\sqrt{s_{NN}}$ = 14.6, 19.6, 27, and 54.4 GeV. A monotonous
increase with $p_{T}$ in $v_{2}$ of light nuclei is observed across all four
center-of-mass energies.

Figure 1: $v_{2}(p_{T})$ of light nuclei ($d$, $t$, and ${}^{3}\text{He}$) in
0-80% centrality Au+Au collisions at $\sqrt{s_{NN}}$ = 14.6, 19.6, 27, and
54.4 GeV. Vertical lines and shaded bands at each marker represent statistical
and systematic uncertainties, respectively.

###  3.2 Centrality dependence of $v_{2}$

Centrality dependence of $v_{2}$ of $d$ is shown in Fig. 2. The nuclei $v_{2}$
is measured in two centrality ranges 0-30% and 30-80% for Au+Au collisions at
$\sqrt{s_{NN}}$ = 19.6, 27, and 54.4 GeV. It is noted that peripheral
collisions exhibit higher $v_{2}$ values compared to more central collisions.
This observation can be attributed to the greater spatial anisotropy in
peripheral collisions as opposed to central collisions.

Figure 2: Centrality dependence of $v_{2}$ of $d$ as a function of $p_{T}$ in
Au+Au collisions at $\sqrt{s_{NN}}$ = 19.6, 27, and 54.4 GeV. Vertical lines
and shaded bands at each marker represent statistical and systematic
uncertainties, respectively.

###  3.3 Mass number scaling

According to the coalescence model, assuming that protons and neutrons behave
in the same way, for a light nuclei $N$ with mass number $A$, we expect
$v_{2,N}(p_{T})\approx Av_{2,p}(p_{T}/A)$, where $v_{2,p}$ is elliptic flow of
protons [8, 4, 9]. The phenomenon is referred to as mass number scaling.
Figure 3 shows the comparison of $v_{2}/A$ of light nuclei as a function of
$p_{T}/A$ (where $A$ is the mass number of the nuclei) with $v_{2}/A$ of
proton (where $A=1$). Proton $v_{2}$ has been fitted with a third-order
polynomial. The bottom panel in each plot shows the ratio between the
$v_{2}/A$ of light nuclei and the fit to proton $v_{2}$. It is observed that
$v_{2}$ of light nuclei deviates from mass number scaling by 20-30%. However,
additional model studies are required to conclude whether the coalescence
model is the dominant production mechanism of light nuclei.

Figure 3: Mass number scaling of $v_{2}/A$ of light nuclei as a function of
$p_{T}/A$ in minimum bias Au+Au collisions at $\sqrt{s_{NN}}$ = 14.6, 19.6,
27, and 54.4 GeV. Vertical lines and shaded bands at each marker represent
statistical and systematic uncertainties, respectively.

##  4 Conclusion

In summary, we have reported the $v_{2}$ of $d$, $t$, and ${}^{3}\text{He}$ in
Au+Au collisions at $\sqrt{s_{NN}}$ = 14.6, 19.6, 27, and 54.4 GeV. A
monotonic rise of light nuclei $v_{2}$ with $p_{T}$ is observed for all light
nuclei species and studied energies. $v_{2}$ of $d$ is observed to show a
strong centrality dependence being higher for peripheral collisions compared
to central collisions. This behaviour can be attributed to the fact that
peripheral collisions have higher spatial anisotropy compared to the central
ones. In addition, it is also observed that $v_{2}$ of light nuclei deviates
from mass number scaling by 20-30%.

## References

  * [1] Dmytro Oliinychenko  “Overview of light nuclei production in relativistic heavy-ion collisions” The 28th International Conference on Ultra-relativistic Nucleus-Nucleus Collisions: Quark Matter 2019  In _Nuclear Physics A_ 1005, 2021, pp. 121754  DOI: https://doi.org/10.1016/j.nuclphysa.2020.121754
  * [2] Anton Andronic, Peter Braun-Munzinger, Krzysztof Redlich and Johanna Stachel  “Decoding the phase structure of QCD via particle production at high energy”  In _Nature_ 561.7723 Springer ScienceBusiness Media LLC, 2018, pp. 321–330  DOI: 10.1038/s41586-018-0491-6
  * [3] S.. Butler and C.. Pearson  “Deuterons from High-Energy Proton Bombardment of Matter”  In _Phys. Rev._ 129 American Physical Society, 1963, pp. 836–842  DOI: 10.1103/PhysRev.129.836
  * [4] T.Z. Yan et al.  “Scaling of anisotropic flow and momentum-space densities for light particles in intermediate energy heavy ion collisions”  In _Physics Letters B_ 638.1, 2006, pp. 50–54  DOI: https://doi.org/10.1016/j.physletb.2006.05.018
  * [5] K.H. Ackermann et al.  “STAR detector overview” The Relativistic Heavy Ion Collider Project: RHIC and its Detectors  In _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 499.2, 2003, pp. 624–632  DOI: https://doi.org/10.1016/S0168-9002(02)01960-5
  * [6] W.J. Llope  “The large-area time-of-flight upgrade for STAR”  In _Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms_ 241.1, 2005, pp. 306–310  DOI: https://doi.org/10.1016/j.nimb.2005.07.089
  * [7] A.. Poskanzer and S.. Voloshin  “Methods for analyzing anisotropic flow in relativistic nuclear collisions”  In _Phys. Rev. C_ 58 American Physical Society, 1998, pp. 1671–1678  DOI: 10.1103/PhysRevC.58.1671
  * [8] L. Adamczyk et al.  “Measurement of elliptic flow of light nuclei at $\sqrt{{s}_{NN}}=200$, 62.4, 39, 27, 19.6, 11.5, and 7.7 GeV at the BNL Relativistic Heavy Ion Collider”  In _Phys. Rev. C_ 94 American Physical Society, 2016, pp. 034908  DOI: 10.1103/PhysRevC.94.034908
  * [9] Yongseok Oh and Che Ming Ko  “Elliptic flow of deuterons in relativistic heavy-ion collisions”  In _Phys. Rev. C_ 76 American Physical Society, 2007, pp. 054910  DOI: 10.1103/PhysRevC.76.054910