# Gradient flow structure and convergence analysis of the ensemble Kalman
inversion for nonlinear forward models

Simon Weissmann Universität Heidelberg, Interdisziplinäres Zentrum für
Wissenschaftliches Rechnen, D-69120 Heidelberg, Germany (simon.weissmann@uni-
heidelberg.de)

###### ?abstractname?

The ensemble Kalman inversion (EKI) is a particle based method which has been
introduced as the application of the ensemble Kalman filter to inverse
problems. In practice it has been widely used as derivative-free optimization
method in order to estimate unknown parameters from noisy measurement data.
For linear forward models the EKI can be viewed as gradient flow
preconditioned by a certain sample covariance matrix. Through the
preconditioning the resulting scheme remains in a finite dimensional subspace
of the original high-dimensional (or even infinite dimensional) parameter
space and can be viewed as optimizer restricted to this subspace. For general
nonlinear forward models the resulting EKI flow can only be viewed as gradient
flow in approximation. In this paper we discuss the effect of applying a
sample covariance as preconditioning matrix and quantify the gradient flow
structure of the EKI by controlling the approximation error through the spread
in the particle system. The ensemble collapse on the one side leads to an
accurate gradient approximation, but on the other side to degeneration in the
preconditioning sample covariance matrix. In order to ensure convergence as
optimization method we derive lower as well as upper bounds on the ensemble
collapse. Furthermore, we introduce covariance inflation without breaking the
subspace property intending to reduce the collapse rate of the ensemble such
that the convergence rate improves. In a numerical experiment we apply EKI to
a nonlinear elliptic boundary-value problem and illustrate the dependence of
EKI as derivative-free optimizer on the choice of the initial ensemble.

Keywords. Ensemble Kalman inversion, Tikhonov regularization, derivative-free
optimization, subspace property, covariance inflation
AMS(MOS) subject classifications. 37C10, 65N21, 65N75, 93D05

######  ?contentsname?

  1. 1 Introduction
    1. 1.1 Literature review
    2. 1.2 Main contributions
  2. 2 Mathematical setup
    1. 2.1 The ensemble Kalman filter applied to inverse problems
  3. 3 Main result: Convergence of TEKI for nonlinear forward models
  4. 4 Preliminaries
    1. 4.1 Existence of solutions and ensemble collapse
    2. 4.2 Gradient approximation via Taylor expansion
    3. 4.3 TEKI as constrained optimization method
  5. 5 Proof of Theorem 3.2
  6. 6 Stabilization via covariance inflation: A generalized ensemble square root filter for inverse problems
    1. 6.1 Incorporation of covariance inflation
  7. 7 Numerical experiment
    1. 7.1 Fixed regularization
    2. 7.2 Adaptive regularization
  8. 8 Conclusion

##  1 Introduction

In inverse problems we are concerned with the task of recovering some unknown
quantity of interest, which can not be observed directly. Its research area
has a wide range of applications in science and engineering, including among
others medical imaging and geophysics. Since inverse problems are typically
ill-posed, much research focuses on theoretical and practical analysis of
regularization methods. Regularization techniques help to overcome instability
issues arising for example through noise in the data. Beside deterministic
regularization methods, the Bayesian approach formulates inverse problems in a
statistical framework. By formulating a probabilistic approach, it
incorporates uncertainty of the underlying model. The regularization can then
be viewed as incorporating prior information through a probability
distribution on the unknown parameter. The solution of the Bayesian inverse
problem is given by the posterior distribution - the distribution of the
unknown parameter conditioned on the realization of the observation. The
resulting posterior distribution is often not accessible directly, such that
sampling or suitable integration methods are needed.

In this document, we focus on a particle based method which is commonly used
for data assimilation problems - the ensemble Kalman filter (EnKF). The EnKF
has been introduced by Evensen [27, 28] and more recently, formulated for
solving inverse problems [38]. The application of the EnKF to inverse problem
has been established as widely used tool and is known as ensemble Kalman
inversion (EKI). The aim of this manuscript is to analyse EKI as derivative-
free optimization method for general nonlinear forward models. Our analysis
will be based on the continuous-time formulation of the scheme, which can be
viewed as (stochastic) ordinary differential equation ((S)ODE) in time, with
focus on the deterministic setting, i.e. without perturbed observations. We
will quantify the approximation of derivatives through the ensemble spread and
derive sufficient lower and upper bounds on the collapse in order to verify
convergence of the scheme.

###  1.1 Literature review

The EnKF has been introduced by Evensen [27] as particle based approximation
of the filtering distribution arising in data assimilation problems. The
method has been applied in the context of Bayesian inverse problems [16, 23].
An important focus lies in the analysis of the large ensemble size limit,
which has been studied for linear and Gaussian models [48, 42] as well as in a
nonlinear setting [47] and for the ensemble square root filter [44]. Viewing
the particle system as Monte Carlo method had led to much focus on the
formulation of multilevel variants [33, 17, 34, 11]. Moreover, ensemble Kalman
methods have been studied in the long time behavior [39, 59, 40]. The accuracy
for a fixed ensemble size has been studied in [60, 49] for the EnKF and in
[19, 18] for the ensemble Kalman-Bucy filter. Furthermore, a high focus of
research for ensemble Kalman methods lies in the continuous-time formulation
[4, 5, 51, 46, 45, 43].

In [38] the authors propose the application of the EnKF to inverse problems.
In the literature EKI has been analysed as sequential Monte Carlo type methods
as well as derivative-free optimization method. In the setting of linear
forward maps and Gaussian prior assumptions it is well-known that EKI
approximates the posterior distribution in its mean field limit. However, in
[26] it has been demonstrated that in general nonlinear settings the resulting
estimator is not consistent with respect to the posterior distribution. In the
mean field limit EKI has been analysed based on the Fokker–Planck equation
[20, 32]. The continuous-time limit of EKI has been formally derived in [56]
and theoretically analysed in [8, 7]. In the continuous-time formulation EKI
can be viewed as derivative-free optimization method [56] and has shown
promising results for the training task in different machine learning
applications [41, 31]. Building up on the continuous-time formulation there
has been much analysis on deterministic EKI, which ignores the diffusion of
the underlying SDE, and stochastic EKI including the perturbed observations.
For linear models well-posedness and convergence results have been derived for
the deterministic formulation [56] and the stochastic formulation [9]. The
nonlinear setting has been analysed in [14] where the authors consider the
discrete time setting with a non constant step size scheme. The authors
include variance inflation breaking the subspace property in order to verify
convergence of the method. Moreover, EKI has been studied for nonlinear models
in the mean field limit [22], where weights have been incorporated in order to
correct the resulting posterior estimate. In our presented analysis it turns
out that the gradient flow structure highly depends on the behavior of the
sample covariance which is used as preconditioner. The dynamical behavior of
EKI and its sample covariance has been described and analysed by a spectral
decomposition [10]. Furthermore, in [61] a localization of the sample
covariance has been introduced and analysed based on the deterministic
continuous-time formulation for nonlinear forward models. Since EKI is applied
as solver for inverse problems, one has to handle noise in the data. In this
context, an early stopping criterion has been proposed in [57] and discrete
regularization has been analysed in [36, 37]. The focus in this document is on
the Tikhonov regularized modification of EKI (TEKI) which has been introduced
in [13] and further analysed in [62]. Further, adaptive regularization methods
within EKI have been studied in [50, 35].

From an alternative point of view, EKI has been modified to a particle based
sampling method [29, 30, 21, 53]. Shifting the noise in the observations to
the parameter space itself the resulting SDE can be treated as Langevin
dynamic, where the ergodicity can be used to build a sampling method.

###  1.2 Main contributions

In this article we are going to quantify EKI as derivative-free optimization
method in the context of general nonlinear forward maps. We employ a gradient
flow structure and present a convergence analysis based on Lyapunov functions.
With this document, we extend the convergence analysis for TEKI presented in
[13] where the convergence result was based on linear forward models. We make
the following contributions:

  * •

We present the well-posedness of the TEKI flow by proving unique existence of
solutions for (10). Furthermore, we view the TEKI flow as approximate gradient
flow by decomposing the flow into the preconditioned gradient direction and
the approximation error resulting through Taylor expansion.

  * •

We quantify the approximation error through the spread in the particle system.
While the ensemble collapses in time with rate $1/t$, we will see that the
approximation error degenerates with rate $(1/t)^{3/2}$. In contrast, we prove
that the sample covariance remains strictly-positive definite as operator
acting on the subspace spanned by the initial ensemble spread
$\mathcal{B}=u_{0}^{\perp}+\operatorname{span}\\{u_{0}^{(j)}-\bar{u}_{0}\mid
j=1,\dots,J\\}$ for a specific $u_{0}^{\perp}\in\mathcal{X}$.

  * •

We describe the effect of preconditioning through the sample covariance.
Therefore, we view the minimization task of $\Phi_{R}$ over the subspace
$\mathcal{B}$ as equality constrained optimization method and derive a Polyak-
Łojasiewicz (PL)-type inequality restricted to the subspace $\mathcal{B}$.
Under strong convexity and smoothness of $\Phi_{R}$ we are then able to prove
convergence with rate $(1/t)^{\frac{1}{\alpha}}$, $\alpha\geq 1$.

  * •

We incorporate covariance inflation without breaking the subspace property.
The inflation sufficiently slows down the ensemble collapse such that the
convergence rate $\alpha$ can be improved to $(1-\rho)\alpha$, where
$\rho\in[0,1)$ is the inflation factor. Furthermore, we view our proposed
covariance inflation scheme as generalization of the ensemble square root
filter (ESRF) applied to inverse problems.

#### Outline

The remainder of this manuscript is structured as follows. In Section 2 we
provide the mathematical background for (T)EKI, and we formulate our main
convergence result with the corresponding assumptions in Section 3. We provide
a row of auxiliary properties of the TEKI flow in Section 4 which are then
applied to present the proof of our main Theorem 3.2 in Section 5. The
incorporation of the covariance inflation and the improved convergence
analysis are presented in Section 6. In Section 7 we illustrate the presented
convergence results in a one-dimensional elliptic boundary-value problem. The
document will be closed with a summary and outlook for open research
directions in Section 8.

##  2 Mathematical setup

In the following we will introduce the mathematical background. We consider
the inverse problem of recovering $u\in\mathcal{X}$ given noisy observations

(1) $y=G(u)+\eta,$

where $\mathcal{X}$ is some Hilbert space, $G:\mathcal{X}\to\mathbb{R}^{K}$ is
the possibly nonlinear forward map and $\eta$ denotes additive Gaussian noise,
i.e. $\eta\sim\mathcal{N}(0,\Gamma)$ for some symmetric positive definite
noise covariance matrix $\Gamma\in\mathbb{R}^{K\times K}$. Due to ill-
posedness of inverse problems, we typically consider the task of minimization
of a regularized loss functional of the form

$\min_{x\in\mathcal{X}}\
\mathcal{L}_{\mathbb{R}^{K}}(G(x),y)+\mathcal{R}_{\mathcal{X}}(x),$

where we describe the discrepancy to the data through
$\mathcal{L}_{\mathbb{R}^{K}}:\mathbb{R}^{K}\times\mathbb{R}^{K}\to\mathbb{R}_{+}$
and incorporate prior information through the regularization function
$\mathcal{R}_{\mathcal{X}}:\mathcal{X}\to\mathbb{R}_{+}$. We refer to [24, 3]
for an overview of various types of regularization methods. Classical choices
of regularization include Tikhonov regularization [25] and total variation
regularization [15, 54]. Throughout this document, our focus will be on the
finite dimensional case $\mathcal{X}=\mathbb{R}^{n_{x}}$ and the Tikhonov
regularized lossfunction of the form

(2) $\Phi_{R}(x)=\Phi(x;y)+R(x),\quad R(x)=\frac{1}{2}\|x\|_{C_{0}}^{2},$

where we define the least-square data misfit by

(3) $\Phi(x;y)=\frac{1}{2}\|G(x)-y\|_{\Gamma}^{2}$

and $R$ denotes the regularization function with symmetric positive definite
regularization matrix $C_{0}\in\mathbb{R}^{n_{x}\times n_{x}}$. Here, we have
introduced the notation $\|\cdot\|_{\Gamma}:=\|\Gamma^{-1/2}\cdot\|$,
$\|\cdot\|_{C_{0}}=\|C_{0}^{-1/2}\cdot\|$, where $\|\cdot\|$ denotes the
euclidean norm in $\mathbb{R}^{K}$ and $\mathbb{R}^{n_{x}}$ respectively.
Moreover, we will make use of $\|\cdot\|_{\mathrm{HS}}$ as notation for the
Hilbert-Schmidt operator norm.

The Tikhonov regularization is closely connected to the Bayesian approach for
inverse problems [58] which incorporates regularization from a statistical
point of view. In a probabilistic model we let $u$ be an $\mathcal{X}$-valued
random variable with prior distribution $\mu_{0}$ stochastically independent
of the noise $\eta$. When viewing $(u,y)$ as a jointly distributed random
variable on $\mathbb{R}^{K}\times\mathcal{X}$, the solution of the Bayesian
inverse problem is given by the posterior distribution of $u\mid y$

(4) $\mu(\mathrm{d}x)=\frac{1}{Z}\exp(-\Phi(x;y))\mu_{0}({\mathrm{d}}x),$

with normalization constant

$Z:=\int_{\mathbb{R}^{n_{x}}}\exp(-\Phi(x;y))\mu_{0}({\mathrm{d}}x).$

We note for a Gaussian prior assumption $\mu_{0}=\mathcal{N}(0,C_{0})$ the
maximum a-posteriori estimate computes as

$\min_{x\in\mathcal{X}}\ \Phi(x;y)+\frac{1}{2}\|x\|_{C_{0}}^{2},$

which relates the Bayesian approach for inverse problems to the Tikhonov
regularization through (2). Throughout this manuscript we assume that
$y\in\mathbb{R}^{n_{y}}$ is fixed and suppress the dependence of $\Phi$ on $y$
by writing $\Phi(x)=\Phi(x;y)$.

###  2.1 The ensemble Kalman filter applied to inverse problems

The EKI has been introduced as the application of the EnKF to inverse problems
[38]. We will follow the derivation of EKI as sequential Monte Carlo method
aiming to approximate the posterior distribution. Therefore, we introduce the
tempered distribution

(5)
$\mu_{n+1}({\mathrm{d}}u)=\frac{1}{Z_{n}}\exp(-h\Phi(u;y))\mu_{n}({\mathrm{d}}u),{\quad
n=1,\dots,N,\quad N\in\mathbb{N},}$

with $h=1/N$ and normalizing constants $Z_{n}$, where $\mu_{0}$ is the prior
distribution and $\mu_{N}$ corresponds to the posterior distribution. This
step incorporates an artificial discrete time system shifting weight from
prior to posterior distribution. Given a sample from the prior distribution,
the EKI evolves the sample as particle system through Gaussian approximation
steps. Let $(u_{0}^{(j)})_{j\in\\{1,\dots,J\\}}$ be the initial ensemble of
size $J\geq 2$ and consider the empirical approximation of the tempering
distribution defined in (5)

$\mu_{n}(\mathrm{d}u)\approx\frac{1}{J}\sum_{j=1}^{J}\delta_{u_{n}^{(j)}}(\mathrm{d}u).$

Defining the following empirical means and covariances

$\displaystyle C(u)$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}(u^{(j)}-\bar{u})\otimes(u^{(j)}-\bar{u}),\quad\bar{u}=\frac{1}{J}\sum_{j=1}^{J}u^{(j)},$
$\displaystyle C^{G,G}(u)$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}(G(u^{(j)})-\bar{G})\otimes(G(u^{(j)})-\bar{G}),\quad\bar{G}=\frac{1}{J}\sum_{j=1}^{J}G(u^{(j)}),$
$\displaystyle C^{u,G}(u)$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}(u^{(j)}-\bar{u})\otimes(G(u^{(j)})-\bar{G}),$
$\displaystyle C^{G,u}(u)$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}(G(u^{(j)})-\bar{G})\otimes(u^{(j)}-\bar{u}),$

the ensemble Kalman iteration for the current particle system
$(u_{n}^{(j)})_{j\in\\{1,\dots,J\\}}$ is given by

(6)
$u_{n+1}^{(j)}=u_{n}^{(j)}+C^{u,G}(u_{n})(C^{G,G}(u_{n})+h^{-1}\Gamma)^{-1}(y_{n+1}^{(j)}-G(u_{n}^{(j)})),\quad
j=,1\dots,J,$

where $h>0$ is a artificial step size and $y_{n+1}^{(j)}$ are perturbed
observations

$y_{n+1}^{(j)}=y+\xi_{n+1}^{(j)},\quad\xi_{n+1}^{(j)}\sim\mathcal{N}(0,h^{-1}\Gamma)$

with $\xi_{n+1}^{(j)}$ being i.i.d. with respect to both $j$ and $n$.
Considering the limit $h\to 0$ the discrete EKI (6) can be viewed as time
discretization of the system of coupled SDEs

(7)
${\mathrm{d}}u_{t}^{(j)}=C^{u,G}(u_{t})\Gamma^{-1}(y-G(u_{t}^{(j)}))\,{\mathrm{d}}t+C^{u,G}(u_{t})\Gamma^{-1/2}\,{\mathrm{d}}W_{t}^{(j)},$

where $W^{(j)}$ are independent Brownian motions on $\mathbb{R}^{K}$. The
continuous-time limit has been formulated in [56] and theoretical verified in
[8, 7].

Since the EKI is not consistent with the posterior distribution for general
nonlinear forward maps [26], we will focus on EKI as derivative-free
optimization method. Instead of considering perturbed observations, sometimes
also referred to as stochastic EKI, the EKI is often analysed in its
deterministic formulation represented by a system of coupled ODEs

(8)
$\frac{{\mathrm{d}}u_{t}^{(j)}}{{\mathrm{d}}t}=C^{u,G}(u_{t})\Gamma^{-1}(y-G(u_{t}^{(j)})),$

which describes the drift term of the SDE formulation [56, 57]. Considering a
linear forward map $G(\cdot)=A\cdot$ with
$A\in\mathcal{L}(\mathcal{X},\mathbb{R}^{K})$ the deterministic EKI can be
viewed as preconditioned gradient flow written through

(9)
$\frac{{\mathrm{d}}u_{t}^{(j)}}{{\mathrm{d}}t}=-C(u_{t})\nabla\Phi(u_{t}^{(j)}).$

For general nonlinear forward maps the representation as preconditioned
gradient flow holds only approximatively. In particular, using Taylor
approximations we can justify that

$C^{u,G}(u)\approx C(u){\mathrm{D}}G^{\ast}(u),$

where ${\mathrm{D}}G$ denotes the (Fréchet) derivative of $G$, and we are
going to quantify the error between the nonlinear EKI and a preconditioned
gradient flow. Since our aim is to view the EKI as optimization method for
minimizing the data misfit functional $\Phi$, the resulting scheme is ill-
posed due to the ill-posed inverse and it is natural to consider a regularized
modification of the EKI. Therefore, the authors in [13] proposed to
incorporate Tikhonov regularization within EKI resulting in the deterministic
continuous-time formulation as coupled system of ODEs

(10)
$\frac{{\mathrm{d}}u_{t}^{(j)}}{{\mathrm{d}}t}=C^{u,G}(u_{t})\Gamma^{-1}(y-G(u_{t}^{(j)}))-C(u_{t})C_{0}^{-1}u_{t}^{(j)},$

which can then be viewed as derivative-free approximation to the
preconditioned gradient flow

(11)
$\frac{{\mathrm{d}}u_{t}^{(j)}}{{\mathrm{d}}t}=-C(u_{t})\nabla\Phi_{R}(u_{t}^{(j)}).$

We will see that this approximation is getting more accurate for a decreasing
spread of the particle system. To do so, we are going to analyse the long-time
behavior of (10) for a fixed number of particles. Therefore, we will extend
the convergence analysis presented in [13] from the linear to the nonlinear
setting. Here, it has been shown that assuming a linear forward model
$G(\cdot)=A\cdot,$ $A\in\mathcal{L}(\mathcal{X},\mathbb{R}^{K})$ the resulting
solution of the TEKI flow (10) minimizes the objective function $\Phi_{R}$ in
the subspace spanned by initial ensemble spread
$\mathcal{B}=u_{0}^{\perp}+\operatorname{span}\\{u_{0}^{(j)}-\bar{u}_{0}\mid
j=1,\dots,J\\}$ for a specific $u_{0}^{\perp}\in\mathcal{X}$.

##  3 Main result: Convergence of TEKI for nonlinear forward models

In order to quantify the nonlinear EKI as derivative-free approximation we
make the following assumptions.

######  Assumption .

We assume that the objective function $\Phi_{R}$ is
$C^{2}(\mathcal{X};\mathbb{R}_{+})$ and $\mu$-strongly convex, i.e.

$\Phi_{R}(x_{1})-\Phi_{R}(x_{2})\geq\langle\nabla\Phi_{R}(x_{2}),x_{1}-x_{2}\rangle+\frac{\mu}{2}\|x_{1}-x_{2}\|^{2}$

for some $\mu>0$ and all $x_{1},x_{2}\in\mathcal{X}$. Moreover, we assume that
$\Phi_{R}$ is $L$-smooth, i.e.

$\Phi_{R}(x_{1})-\Phi_{R}(x_{2})\leq\langle\nabla\Phi_{R}(x_{2}),x_{1}-x_{2}\rangle+\frac{L}{2}\|x_{1}-x_{2}\|^{2}$

for some $L>0$ and all $x_{1},x_{2}\in\mathcal{X}$ and therefore satisfies the
PL-inequality

$\nu\|\nabla\Phi(x)\|^{2}\geq\Phi_{R}(x)-\Phi(x_{\ast}),\quad x\in\mathcal{X}$

for some $\nu>0$, where $x_{\ast}$ is a stationary point, i.e.
$\nabla\Phi(x_{\ast})=0$.

In the above assumption, we suppose global $\mu$-strong convexity and
$L$-smoothness, which is needed to prove convergence results as global
optimizer. In the case, where those assumptions are only satisfied locally, we
can only expect the results to hold in local neighborhoods around critical
points. The necessary convergence analysis as local optimizer is work in
progress.

######  Assumption .

We assume that the forward map $G$ is $C^{2}(\mathcal{X};\mathbb{R}^{n_{y}})$,
locally Lipschitz with constant $c_{\mathrm{lip}}>0$ and can be approximated
linearly by

(12)
$G(x_{1})=G(x_{2})+{\mathrm{D}}G(x_{2})(x_{1}-x_{2})+{\mathrm{Res}}(x_{1},x_{2}),$

where the approximation error ${\mathrm{Res}}$ is bounded by

(13) $\|{\mathrm{Res}}(x_{1},x_{2})\|_{\mathbb{R}^{n_{y}}}\leq
b_{\mathrm{res}}\|x_{1}-x_{2}\|_{\mathcal{X}}^{2}$

for some $b_{\mathrm{res}}>0$ independent from $x_{1},x_{2}\in\mathcal{X}$.

We note that Assumption 3 is satisfied for linear forward maps
$G(\cdot)=A\cdot$, $A\in\mathcal{L}(\mathcal{X},\mathbb{R}^{K})$ with
${\mathrm{Res}}(x_{1},x_{2})=0$ for all $x_{1},x_{2}\in\mathcal{X}$.
Furthermore, it is sufficient to assume that the Hessian of $G$ is uniformly
bounded in the sense that $\langle z,{\mathrm{D}}^{2}G(x)z\rangle\leq
c\|z\|^{2}$ for all $z\in\mathcal{X}$. Alternatively, the tangential cone
condition

$\|G(x_{1})-G(x_{2})-{\mathrm{D}}G(x_{2})(x_{1}-x_{2})\|\leq\tilde{b}\|x_{1}-x_{2}\|\|G(x_{1})-G(x_{2})\|$

for some uniform $\tilde{b}>0$ together with Lipschitz continuity implies
Assumption 3. This condition is often assumed (locally) when studying
iterative methods for nonlinear inverse problems [55]. As soon as we are able
to proof uniform boundedness of the solution of the TEKI flow (10), it is
sufficient to consider Assumption 3 locally, i.e. to assume that for $B>0$
there exists $b_{\mathrm{res}}(B)>0$ such that (13) holds for all
$x_{1},x_{2}\in\mathcal{X}$ with
$\|x_{1}\|_{\mathcal{X}},\|x_{2}\|_{\mathcal{X}}^{2}\leq B$. While the proof
of boundedness of the TEKI flow (10) is typically challenging, one can enforce
boundedness through a smooth modification of the forward map $G$, such that
$G(x)=0$ for $\|x\|>B$, see for example [13, 14]. Considering this type of
modification might be reasonable for practical implementation if it is known
that solutions of the corresponding inverse problem should be bounded.

Furthermore, we note that the above Assumption 3 can be weaken through the
assumption $\|{\mathrm{Res}}(x_{1},x_{2})\|\leq
b_{\mathrm{res}}\|x_{1}-x_{2}\|^{p}$, for $p>2$. However, to avoid technical
details we will keep the assumption with $p=2$ for the main proof and state
more details on the extension to $p>2$ in Remark 5.2 after the proof of our
main result.

######  Proposition 3.1.

Suppose Assumption 3 is satisfied, then it holds true that

$\displaystyle\frac{{\mathrm{d}}\frac{1}{J}\sum_{j=1}^{J}\left(\Phi_{R}(u_{t}^{(j)})\right)}{{\mathrm{d}}t}$
$\displaystyle\leq-\frac{1}{J}\sum_{j=1}^{J}\langle\nabla\Phi_{R}(u_{t}^{(j)}),C(u_{t})\nabla\Phi_{R}(u_{t}^{(j)})\rangle$
$\displaystyle\quad+b_{1}JV_{e}(t)^{3/2}\frac{1}{J}\sum_{j=1}^{J}\left(\frac{1}{2}\|\nabla\Phi_{R}(u_{t}^{(j)})\|^{2}+\Phi_{R}(u_{t}^{(j)})\right),$

where $b_{1}>0$ is independent from $J$ and
$V_{e}(t):=\frac{1}{J}\sum_{k=1}^{J}\|u_{t}^{(k)}-\bar{u}_{t}\|^{2}$.

The above result quantifies the (approximate) gradient flow structure of TEKI.
We will see, that $V_{e}(t)$ is of the order $1/t$, which means that the
approximation error of the gradient through the covariance degenerates in
time. However, the degeneration of $C(u_{t})$ influences the preconditioning
effect for the gradient as well. Although we can show, that the degeneration
of $C(u_{t})$ is bounded from below by the order $1/t$, it is not possible to
imply lower bounds on the eigenvalues $C(u_{t})$. This issue comes from the
fact that $C(u_{t})$ is a sample covariance and has at most rank
$\min(n_{x},J-1)$. To alleviate this problem, we will show that $C(u_{t})$
will remain positive definite in the subspace spanned through the initial
ensemble spread. Along the solution of the TEKI flow the particle system
$(u_{t}^{(j)})_{j=1,\dots,J}$ remains in the subspace

$\mathcal{B}:=u_{0}^{\perp}+\operatorname{span}\\{u_{0}^{(j)}-\bar{u}_{0}\mid
j=1,\dots,J\\}\subset\operatorname{span}\\{u_{0}^{j}\mid
j=1,\dots,J\\}=:\mathcal{S},$

where $u_{0}^{\perp}=\bar{u}_{0}-\mathcal{P}_{E}\bar{u}_{0}$, where
$\mathcal{P}_{E}:=E(E^{\top}E)^{-1}E^{T}$ with
$E=((u_{0}^{(1)}-\bar{u}_{0})^{\top},\dots,(u_{0}^{(J)}-\bar{u}_{0})^{\top})\in\mathbb{R}^{n_{x}\times
J}$ is the orthogonal projection on the subspace
$\operatorname{span}\\{u_{0}^{(j)}-\bar{u}_{0}\mid j=1,\dots,J\\}$. This
subspace property also holds for general nonlinear forward maps $G$, see Lemma
3.7 and Corollary 3.8 in [13]. Hence, our aim is to quantify TEKI as gradient
flow w.r.t. the constrained optimization problem

(14) $\min_{x\in{\mathcal{B}}}\ \Phi_{R}(x),$

where $(u_{0}^{(j)})_{j=1,\dots,J}$ is the initial ensemble for TEKI. In fact,
a similar result has been shown in [13] for the linear setting, where the
gradient approximation is exact. We are going to extend the theory to
nonlinear forward maps satisfying Assumptions 3 and 3. Note that we will
reformulate the optimization problem (14) as equality constrained optimization
problem, implying that there exists a unique minimizer
$u_{\ast}\in\mathcal{B}$ of (14). In the following we formulate our main
convergence result.

######  Theorem 3.2.

Suppose Assumption 3 and Assumption 3 are satisfied, let $(u_{t}^{(j)},t\geq
0)_{j=1,\dots,J}$ solve (10) with linearly independent initial ensemble
$(u_{0}^{(j)})_{j=1,\dots,J}$ and let $u_{\ast}\in\mathcal{B}$ be the unique
global minimizer of (14). Then there exist $c_{1},c_{2}>0$ such that

(15)
$\frac{1}{J}\sum_{j=1}^{J}\Phi_{R}(u_{t}^{(j)})-\Phi_{R}(u_{\ast})\leq\left(\frac{c_{1}}{t+c_{2}}\right)^{\frac{1}{\alpha}},$

where
$0<\alpha<\frac{L}{\mu}(\sigma_{\max}+c_{\mathrm{lip}}\lambda_{\max}\|C(u_{0})\|_{\mathrm{HS}})$.

We note that the constants $c_{1},c_{2}>0$ depend on various constants arising
in Assumption 3-3 and may be increasing for less well-behaved scenarios.
However, these constants are independent of the time $t>0$ and the convergence
still holds asymptotically.

##  4 Preliminaries

In the following we present a row of auxiliary results, which are needed to
prove the main convergence result in Theorem 3.2.

###  4.1 Existence of solutions and ensemble collapse

Recall that we assume a finite dimensional parameter space
$\mathcal{X}=\mathbb{R}^{n_{x}}$. We start this section with verifying the
well-posedness of the scheme. Therefore, we state that there exists a unique
solution for the considered TEKI flow.

######  Theorem 4.1.

Let $(u_{0}^{(j)})_{j=1,\dots,J}$ be the linearly independent initial ensemble
and define its linear span $S=\operatorname{span}\\{u_{0}^{(j)}\mid
j=1,\dots,J\\}$. Then for all $T>0$ the equation (10) has a unique solution
$(u_{t}^{(j)},t\in[0,T])\in C([0,T],\mathcal{S})$. In particular, it holds
true that $\langle z,u_{t}^{(j)}\rangle=0$ for all $z\in S^{\perp}$ and all
$t\in[0,T]$, $j=1,\dots,J$.

######  Proof 4.2.

Due to local Lipschitz continuity local solutions follow from standard ODE
theory, see e.g. [13, Theorem 3.2] for details. Moreover, we refer to [13,
Corollary 3.8] for the subspace property, i.e. while a solution exists it
holds true that $\langle z,u_{t}^{(j)}\rangle=0$ for all $z\in S^{\perp}$ and
all $j=1,\dots,J$. In the following, we extend the existence of solutions
globally. Therefore, we consider the Lyapunov function

$V_{\bar{u}}(u_{t})=\varphi(\|\bar{u}_{t}\|^{2})+\frac{B_{1}}{J}\sum_{j=1}^{J}\|u_{t}^{(j)}-\bar{u}_{t}\|^{2}+B_{2}$

for $\varphi:\mathbb{R}\to\mathbb{R}$ with $\varphi(z)=\log(1+z)$ and
sufficiently large $B_{1},B_{2}>0$. We firstly observe that
$z\varphi^{\prime}(z)=\frac{z}{1+z}\leq 1$ for $z\geq 0$. The time evolution
of $\varphi(\|\bar{u}_{t}\|^{2})$ is given by

$\displaystyle\frac{\mathrm{d}}{{\mathrm{d}}t}\varphi(\|\bar{u}_{t}\|^{2})$
$\displaystyle=2\varphi^{\prime}(\|\bar{u}_{t}\|^{2})\left(-\langle\bar{u}_{t},C^{u,G}(u_{t})\Gamma^{-1}(\bar{G}_{t}-y)\rangle-\langle\bar{u}_{t},C(u_{t})C_{0}^{-1}\bar{u}_{t}\rangle\right)$
$\displaystyle\leq
2\varphi^{\prime}(\|\bar{u}_{t}\|^{2})\left(\|\bar{u}_{t}\|\|\bar{G}_{t}-y\|_{\Gamma}\|C^{u,G}(u_{t})\Gamma^{-1/2}\|_{\mathrm{HS}}+\|\bar{u}_{t}\|\|\bar{u}_{t}\|_{C_{0}}\|C(u_{t})C_{0}^{-1/2}\|_{\mathrm{HS}}\right)$
$\displaystyle\leq\varphi^{\prime}(\|\bar{u}_{t}\|^{2})^{2}\|\bar{u}_{t}\|^{2}(2\|\bar{G}_{t}-G(\bar{u}_{t})\|_{\Gamma}^{2}+2\|G(\bar{u}_{t})-y\|_{\Gamma}^{2})+\|C^{u,G}(u_{t})\Gamma^{-1/2}\|_{\mathrm{HS}}^{2}$
$\displaystyle\quad+\varphi^{\prime}(\|\bar{u}_{t}\|^{2})^{2}\|\bar{u}_{t}\|^{2}\|\bar{u}_{t}\|_{C_{0}}^{2}+\|C(u_{t})C_{0}^{-1/2}\|_{\mathrm{HS}}^{2},$

where we have applied Young’s inequality and

$\|\bar{G}_{t}-y\|_{\Gamma}^{2}\leq(\|\bar{G}_{t}-G(\bar{u}_{t})\|_{\Gamma}+\|G(\bar{u}_{t})-y\|_{\Gamma})^{2}\leq
2\|\bar{G}_{t}-G(\bar{u}_{t})\|_{\Gamma}^{2}+2\|G(\bar{u}_{t})-y\|_{\Gamma})^{2}$

for the last inequality. With Lipschitz continuity of $G$, Cauchy-Schwarz
inequality and Jensen’s inequality we have that

$\displaystyle\|\bar{G}-G(\bar{u})\|^{2}=\frac{1}{J^{2}}\sum_{j,l}^{2}\langle
G(u^{(j)})-G(\bar{u}),G(u^{(l)})-G(\bar{u})\rangle$ $\displaystyle\leq
c_{\mathrm{lip}}^{2}\left(\frac{1}{J}\sum_{j=1}^{J}\|u^{(j)}-\bar{u}\|\right)^{2}$
$\displaystyle\leq
c_{\mathrm{lip}}^{2}\frac{1}{J}\sum_{j=1}^{J}\|u^{(j)}-\bar{u}\|^{2}.$

Hence, using $z\varphi^{\prime}(z)\leq 1$ we obtain

$\displaystyle\frac{\mathrm{d}}{{\mathrm{d}}t}\varphi(\|\bar{u}_{t}\|^{2})$
$\displaystyle\leq
c_{1}\frac{1}{J}\sum_{j=1}^{J}\|u_{t}^{(j)}-\bar{u}_{t}\|^{2}+\tilde{c}_{2}\varphi^{\prime}(\|\bar{u}\|^{2})\Phi_{R}(\bar{u}_{t})$
$\displaystyle\quad+\|C^{u,G}(u_{t})\Gamma^{-1/2}\|_{\mathrm{HS}}^{2}+\|C(u_{t})C_{0}^{-1/2}\|_{\mathrm{HS}}^{2}$
$\displaystyle\leq
c_{1}\frac{1}{J}\sum_{j=1}^{J}\|u_{t}^{(j)}-\bar{u}_{t}\|^{2}+c_{2}+\|C^{u,G}(u_{t})\Gamma^{-1/2}\|_{\mathrm{HS}}^{2}+\|C(u_{t})C_{0}^{-1/2}\|_{\mathrm{HS}}^{2}$

for some constants $c_{1},c_{2}>0$. Here we have used, that $\Phi_{R}(u)\leq
c\|u\|^{2}$ for local Lipschitz continuous $G$. On the other side, we have

$\displaystyle\frac{{\mathrm{d}}\frac{1}{J}\sum_{j=1}^{J}\|u_{t}^{(j)}-\bar{u}_{t}\|^{2}}{{\mathrm{d}}t}$
$\displaystyle=-\frac{2}{J^{2}}\sum_{j,k=1}^{J}\langle
u_{t}^{(j)}-\bar{u}_{t},u_{t}^{(k)}-\bar{u}_{t}\rangle\langle
G(u_{t}^{(k)})-\bar{G}_{t},G(u_{t}^{(j)})-\bar{G}_{t}\rangle_{\Gamma}$
$\displaystyle\quad-\frac{2}{J^{2}}\sum_{j,k=1}^{J}\langle
u_{t}^{(j)}-\bar{u}_{t},u_{t}^{(k)}-\bar{u}_{t}\rangle\langle
u_{t}^{(k)}-\bar{u}_{t},C_{0}^{-1}(u_{t}^{(j)}-\bar{u}_{t})\rangle$
$\displaystyle=-2\|C^{u,G}(u_{t})\Gamma^{-1/2}\|_{\mathrm{HS}}^{2}-2\|C(u_{t})C_{0}^{-1/2}\|_{\mathrm{HS}}^{2}.$

It follows for $B_{1}>c_{1}$ and $B_{2}>c_{2}$ that

$\frac{\mathrm{d}}{{\mathrm{d}}t}V_{\bar{u}}(u_{t})\leq cV_{\bar{u}}(u_{t})$

and by application of Gronwall’s inequality for all $T>0$ the solution of (10)
remains bounded.

We note that a similar type of Lyapunov function $\varphi$ has been used in
[7] in order to show finite logarithmic moments of the stochastic formulation.
While it is not clear whether $L_{2}$ bounds hold for (T)EKI, one is at least
able to derive bounds on logarithmic moments.

In the following we define the spread of the particle system by

(16) $e_{t}^{(j)}=u_{t}^{(j)}-\bar{u}_{t}.$

From Lemma 3.7 and Corollary 3.8 in [13] it is known that the TEKI flow
satisfies the following extended subspace property.

######  Proposition 4.3 (Corollary 3.8, [13]).

Let $(u_{0}^{(j)})_{j=1,\dots,J}$ be the linearly independent initial ensemble
and define the subspace

$\mathcal{B}=u_{0}^{\perp}+\operatorname{span}\\{e_{0}^{(j)}\mid
j=1,\dots,J\\},$

where $u_{0}^{\perp}=\bar{u}_{0}-\mathcal{P}_{E}\bar{u}_{0}$, where
$\mathcal{P}_{E}=E(E^{\top}E)^{-1}E^{T}$ with
$E=((e_{0}^{(1)})^{\top},\dots,(e_{0}^{(J)})^{\top})\in\mathbb{R}^{n_{x}\times
J}$ is the orthogonal projection on the subspace
$\operatorname{span}\\{e_{0}^{(j)}\mid j=1,\dots,J\\}$. Then for $T\geq 0$ the
unique solution $(u_{t}^{(j)},t\in[0,T])$ of the TEKI flow (10) remains in
$\mathcal{B}$, i.e. $u_{t}^{(j)}\in\mathcal{B}$ for all $t\in[0,T]$ and all
$j=1,\dots,J$.

On the one side the spread of the particle system describes the accuracy of
the gradient approximation, which gets more accurate for a less spread
particle system. On the other side, the spread of the particle system
determines the preconditioning effect through the sample covariance. For this
reason we aim for a collapse of the ensemble but not too fast. In our first
result, we prove that the ensemble collapses with rate $1/t$. Recall that
$C_{0}$ is assumed to be a fixed symmetric positive definite matrix, and hence
all of its eigenvalues are strictly positive.

######  Lemma 4.4.

Let $(u_{t}^{(j)},t\geq 0)_{j=1,\dots,J}$ be the solution of (8) initialized
by a linearly independent ensemble $(u_{0}^{(j)})_{j=1,\dots,J}$. Then the
mapping $V_{e}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ with
$t\mapsto\frac{1}{J}\sum_{j=1}^{J}\|e_{t}\|^{2}$ is bounded by

$V_{e}(t)\leq\frac{1}{\frac{2\sigma_{\min}}{J}t+V_{e}(0)^{-1}},$

where $\sigma_{\min}$ denotes the smallest eigenvalue of $C_{0}^{-1}$.

######  Proof 4.5.

The evolution of $V_{e}(t)$ is given by

$\displaystyle\frac{{\mathrm{d}}\frac{1}{J}\sum_{j=1}^{J}\|u_{t}^{(j)}-\bar{u}_{t}\|^{2}}{{\mathrm{d}}t}$
$\displaystyle=-\frac{2}{J^{2}}\sum_{j,k=1}^{J}\langle
u_{t}^{(j)}-\bar{u}_{t},u_{t}^{(k)}-\bar{u}_{t}\rangle\langle
G(u_{t}^{(k)})-\bar{G}_{t},G(u_{t}^{(j)})-\bar{G}_{t}\rangle_{\Gamma}$
$\displaystyle\quad-\frac{2}{J^{2}}\sum_{j,k=1}^{J}\langle
u_{t}^{(j)}-\bar{u}_{t},u_{t}^{(k)}-\bar{u}_{t}\rangle\langle
u_{t}^{(k)}-\bar{u}_{t},C_{0}^{-1}(u_{t}^{(j)}-\bar{u}_{t})\rangle$
$\displaystyle=-2\|C^{u,G}(u_{t})\Gamma^{-1/2}\|_{\mathrm{HS}}^{2}$
$\displaystyle\quad-\frac{2}{J^{2}}\sum_{j,k=1}^{J}\langle
u_{t}^{(j)}-\bar{u}_{t},u_{t}^{(k)}-\bar{u}_{t}\rangle\langle
u_{t}^{(k)}-\bar{u}_{t},C_{0}^{-1}(u_{t}^{(j)}-\bar{u}_{t})\rangle$
$\displaystyle\leq-\frac{2\sigma_{\min}}{J}\left(\frac{1}{J}\sum_{j=1}^{J}\|u_{t}^{(j)}-\bar{u}_{t}\|^{2}\right)^{2}.$

Hence, we have shown

$\frac{{\mathrm{d}}V_{e}(t)}{{\mathrm{d}}t}\leq-\frac{2\sigma_{\min}}{J}V_{e}(t)^{2}$

and the assertion follows by Lyapunov-type argument.

In the following result, we state that, although the sample covariance
degenerates along the particle evolution, it remains a strictly positive
definite operator acting on the subspace $\mathcal{B}$.

######  Lemma 4.6.

Let $(u_{t}^{(j)},t\geq 0)_{j=1,\dots,J}$ be the solution of (8) initialized
by a linearly independent ensemble $(u_{0}^{(j)})_{j=1,\dots,J}$ such that

(17) $\zeta_{0}=\min_{z\in{\mathcal{B}},\ \|z\|=1}\ \langle
z,C(u_{0})z\rangle>0.$

For each $z\in{\mathcal{B}}$ with $\|z\|=1$ it holds true that

$\langle z,C(u_{t})z\rangle\geq\frac{1}{mt+\zeta_{0}^{-1}},$

where $m=2(c_{\mathrm{lip}}^{2}\lambda_{\max}V_{e}(0)+\sigma_{\max})>0$
depends on the smallest and largest eigenvalue $\sigma_{\max}$ and
$\sigma_{\min}$ of $C_{0}^{-1}$, the largest eigenvalue $\lambda_{\max}$ of
$\Gamma^{-1}$ and the Lipschitz constant $c_{\mathrm{lip}}$ of $G$.

######  Proof 4.7.

We note that the following proof is closely related to the proof of Theorem
3.5 in [13], where the authors derived an upper bound on the sample covariance
matrix. The time evolution of the sample covariance is given by

$\displaystyle\frac{{\mathrm{d}}C(u_{t})}{{\mathrm{d}}t}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\frac{{\mathrm{d}}e_{t}^{(j)}}{{\mathrm{d}}t}\left(e_{t}^{(j)}\right)^{\top}+\frac{1}{J}\sum_{j=1}^{J}e_{t}^{(j)}\left(\frac{{\mathrm{d}}e_{t}^{(j)}}{{\mathrm{d}}t}\right)^{\top}$
$\displaystyle=-2C^{u,G}(u_{t})\Gamma^{-1}C^{G,u}(u_{t})-2C(u_{t})C_{0}^{-1}C(u_{t}).$

Let $z\in{\mathcal{B}}$ with $\|z\|=1$, then we have that

$\displaystyle\langle z,\frac{{\mathrm{d}}C(u_{t})}{{\mathrm{d}}t}z\rangle$
$\displaystyle=-2\|C^{G,u}(u_{t})z\|_{\Gamma}^{2}-2\|C(u_{t})z\|_{C_{0}}^{2}.$

Since $G$ is locally Lipschitz, we obtain with Cauchy-Schwarz

$\displaystyle\|C^{G,G}(u_{t})\|_{\mathrm{HS}}=\operatorname{Tr}(C^{G,G}(u_{t})C^{G,G}(u_{t}))^{1/2}$
$\displaystyle\leq\frac{1}{J}\sum_{j=1}^{J}\|G(u_{t}^{(j)})-\bar{G}\|^{2}$
$\displaystyle\leq c_{\mathrm{lip}}^{2}V_{e}(t)\leq
c_{\mathrm{lip}}^{2}V_{e}(0)$

and therefore it follows with
$\|C^{G,u}(u_{t})z\|^{2}\leq\|C^{G,G}\|_{\mathrm{HS}}\|C(u_{t})z\|_{\mathrm{HS}}$
that

$\|C^{G,u}(u_{t})z\|_{\Gamma}^{2}\leq
c_{\mathrm{lip}}^{2}\lambda_{\max}V_{e}(0)\|C(u_{t})z\|_{C_{0}}^{2}.$

Hence, we have derived

$\frac{{\mathrm{d}}\langle
z,C(u_{t})z\rangle}{{\mathrm{d}}t}\geq-2(c_{\mathrm{lip}}^{2}\lambda_{\max}V_{e}(0)+\sigma_{\max})\|C(u_{t})z\|^{2}.$

Next, we consider $C(u_{t})$ as operator acting on $\mathcal{B}$ and denote
its smallest eigenvalue $\zeta_{t}\geq 0$ with unit-norm eigenvector
$\varphi(t)\in{\mathcal{B}}$. In the following we will prove that $\zeta_{t}$
is indeed strictly positive. Since $\|\varphi(t)\|=1$ for all $t\geq 0$, we
observe that

$0=\frac{{\mathrm{d}}\|\varphi(t)\|^{2}}{{\mathrm{d}}t}=2\langle\varphi(t),\frac{{\mathrm{d}}\varphi(t)}{{\mathrm{d}}t}\rangle.$

Hence, we can describe the time-evolution of the smallest eigenvalue
$\zeta_{t}$ of $C(u_{t})$ (initialized with $\zeta_{0}>0$) through

$\displaystyle\frac{{\mathrm{d}}}{{\mathrm{d}}t}\zeta_{t}=\frac{{\mathrm{d}}}{{\mathrm{d}}t}\langle\varphi(t),C(u_{t})\varphi(t)\rangle$
$\displaystyle=\langle\varphi(t),\frac{{\mathrm{d}}C(u_{t})}{{\mathrm{d}}t}\varphi(t)\rangle+2\langle\frac{{\mathrm{d}}\varphi(t)}{{\mathrm{d}}t},C(u_{t})\varphi(t)\rangle$
$\displaystyle=\langle\varphi(t),\frac{{\mathrm{d}}C(u_{t})}{{\mathrm{d}}t}\varphi(t)\rangle+2\zeta_{t}\langle\frac{{\mathrm{d}}\varphi(t)}{{\mathrm{d}}t},\varphi(t)\rangle$
$\displaystyle=\langle\varphi(t),\frac{{\mathrm{d}}C(u_{t})}{{\mathrm{d}}t}\varphi(t)\rangle$
$\displaystyle\geq-2(c_{\mathrm{lip}}^{2}\lambda_{\max}V_{e}(0)+\sigma_{\max})\|C(u_{t})\varphi(t)\|^{2}$
$\displaystyle\geq-2(c_{\mathrm{lip}}^{2}\lambda_{\max}V_{e}(0)+\sigma_{\max})\zeta_{t}^{2}.$

With the above computation we obtain the lower bound on the smallest
eigenvalue of $C(u_{t})$ as operator acting on $\mathcal{B}$ by

$\zeta_{t}\geq\frac{1}{2(c_{\mathrm{lip}}^{2}\lambda_{\max}V_{e}(0)+\sigma_{\max})t+\zeta_{0}^{-1}}>0,$

and the assertion follows with

(18) $m:=2(c_{\mathrm{lip}}^{2}\lambda_{\max}V_{e}(0)+\sigma_{\max})>0.$

###  4.2 Gradient approximation via Taylor expansion

In the following we will quantify how the TEKI flow (10) approximates
gradients and hence can be viewed as gradient flow.

######  Lemma 4.8.

Suppose that Assumption 3 is satisfied. Given a particle system
$(u^{(j)})_{j=1,\dots,J}$ it holds true that

$\|C^{u,G}(u)\Gamma^{-1}(G(u^{(j)})-y)-C(u)\nabla\Phi(u^{(j)})\|\leq
b_{1}J\sqrt{\Phi(u^{(j)})}\left(\frac{1}{J}\sum_{k=1}^{J}\|u^{(k)}-\bar{u}\|^{2}\right)^{3/2},$

for some $b_{1}$ independent of $J$.

######  Proof 4.9.

We apply (12) and using bi-linearity of the inner product we obtain

$\displaystyle
C^{u,G}(u)\Gamma^{-1}(G(u^{(j)})-y)=\frac{1}{J}\sum_{k=1}^{J}\langle
G(u^{(k)})-\bar{G},G(u^{(j)})-y\rangle_{\Gamma}(u^{k}-\bar{u})$
$\displaystyle=\frac{1}{J}\sum_{k=1}^{J}\langle
G(u^{(j)})+{\mathrm{D}}G(u^{(j)})(u^{(k)}-u^{(j)})+{\mathrm{Res}}(u^{(k)},u^{(j)}),G(u^{(j)})-y\rangle_{\Gamma}(u^{(k)}-\bar{u})$
$\displaystyle\quad-\frac{1}{J}\sum_{l=1}^{J}\langle
G(u^{(j)})+{\mathrm{D}}G(u^{(j)})(u^{(l)}-u^{(j)})+{\mathrm{Res}}(u^{(l)},u^{(j)}),G(u^{(j)})-y\rangle_{\Gamma}(u^{(k)}-\bar{u})$
$\displaystyle=\frac{1}{J}\sum_{k=1}^{J}\langle{\mathrm{D}}G(u^{(j)})(u^{(k)}-\bar{u}),G(u^{(j)})-y\rangle_{\Gamma}(u^{(k)}-\bar{u})$
$\displaystyle\quad+\frac{1}{J}\sum_{k=1}^{J}\langle{\mathrm{Res}}(u^{(k)},u^{(j)})-\overline{\mathrm{Res}}^{(j)},G(u^{(j)})-y\rangle_{\Gamma}(u^{(k)}-\bar{u})$
$\displaystyle=C(u)\nabla\Phi(u^{(j)})+\frac{1}{J}\sum_{k=1}^{J}\langle{\mathrm{Res}}(u^{(k)},u^{(j)})-\overline{\mathrm{Res}}^{(j)},G(u^{(j)})-y\rangle_{\Gamma}(u^{(k)}-\bar{u}),$

where we have defined

$\overline{\mathrm{Res}}^{(j)}=\frac{1}{J}\sum_{l=1}^{J}{\mathrm{Res}}(u^{(l)},u^{(j)}).$

It follows with Assumption 3 and Cauchy-Schwarz that

$\displaystyle\|C^{u,G}(u)\Gamma^{-1}(G(u^{(j)})-y)-C(u)\nabla\Phi(u^{(j)})\|$
$\displaystyle\leq
b_{\mathrm{res}}\|\Gamma^{-1/2}\|\frac{1}{J}\sum_{k=1}^{J}4\left(\|u^{(k)}-u^{(j)}\|^{2}+\frac{1}{J}\sum_{l=1}^{J}\|u^{(l)}-u^{(j)}\|^{2}\right)\|u^{(k)}-\bar{u}\|\|G(u^{(j)})-y\|_{\Gamma}$
$\displaystyle\leq
b_{\mathrm{res}}\|\Gamma^{-1/2}\|\frac{1}{J}\sum_{k=1}^{J}4\left((1+\frac{1}{J})\sum_{l=1}^{J}\|u^{(l)}-u^{(j)}\|^{2}\right)\|u^{(k)}-\bar{u}\|\|G(u^{(j)})-y\|_{\Gamma}$
$\displaystyle\leq
b_{\mathrm{res}}\|\Gamma^{-1/2}\|(8+\frac{2}{J})J\|G(u^{(j)})-y\|_{\Gamma}\left(\frac{1}{J}\sum_{k=1}^{J}\|u^{(k)}-\bar{u}\|^{2}\right)^{3/2},$

where we have used that

$\|u^{(l)}-u^{(j)}\|^{2}\leq
2\|u^{(l)}-\bar{u}\|^{2}+2\|u^{(j)}-\bar{u}\|^{2}.$

We are now ready to prove Proposition 3.1 in order to quantify the gradient
flow structure of TEKI.

######  Proof 4.10 (Proof of Proposition 3.1).

We apply chain rule for deriving the time-evolution of
$V_{\Phi_{R}}(t):=\frac{1}{J}\sum_{j=1}^{J}\Phi_{R}(u_{t}^{(j)})$ as

$\displaystyle\frac{{\mathrm{d}}V_{\Phi_{R}}(t)}{{\mathrm{d}}t}$
$\displaystyle=\frac{{\mathrm{d}}V_{\Phi_{R}}(t)}{{\mathrm{d}}t}=\frac{1}{J}\sum_{j=1}^{J}\langle\nabla\Phi_{R}(u_{t}^{(j)}),\frac{{\mathrm{d}}u_{t}^{(j)}}{{\mathrm{d}}t}\rangle$
$\displaystyle=-\frac{1}{J}\sum_{j=1}^{J}\langle\nabla\Phi_{R}(u_{t}^{(j)}),C(u_{t})\nabla\Phi_{R}(u_{t}^{(j)})\rangle$
$\displaystyle\quad+\frac{1}{J}\sum_{j=1}^{J}\langle\nabla\Phi_{R}(u_{t}^{(j)}),C^{u,G}(u_{t})\Gamma^{-1}(G(u_{t}^{(j)})-y)-C(u_{t})\nabla\Phi(u_{t}^{(j)})\rangle.$

We apply Cauchy-Schwarz inequality and Lemma 4.8 such that

$\displaystyle\frac{{\mathrm{d}}V_{\Phi_{R}}(t)}{{\mathrm{d}}t}$
$\displaystyle\leq-\frac{1}{J}\sum_{j=1}^{J}\langle\nabla\Phi_{R}(u_{t}^{(j)}),C(u_{t})\nabla\Phi_{R}(u_{t}^{(j)})\rangle$
$\displaystyle\quad+b_{1}JV_{e}(t)^{3/2}\frac{1}{J}\sum_{j=1}^{J}\|\nabla\Phi_{R}(u_{t}^{(j)})\|\sqrt{\Phi(u_{t}^{(j)})}$
$\displaystyle\leq-\frac{1}{J}\sum_{j=1}^{J}\langle\nabla\Phi_{R}(u_{t}^{(j)}),C(u_{t})\nabla\Phi_{R}(u_{t}^{(j)})\rangle$
$\displaystyle\quad+b_{1}JV_{e}(t)^{3/2}\frac{1}{J}\sum_{j=1}^{J}\left(\frac{1}{2}\|\nabla\Phi_{R}(u_{t}^{(j)})\|^{2}+\Phi(u_{t}^{(j)})\right).$

The claim follows with $\Phi(x)\leq\Phi_{R}(x)$.

Considering the upper bound on the derived gradient approximation in Lemma
4.8, we observe that it is increasing in the ensemble size $J$. Although we do
not know whether this bound is sharp, it suggests that the gradient
approximation might be less well-behaved in the large ensemble size regime. To
alleviate this issue, we expect improvement when considering localised
preconditioning covariance matrices $C^{u,G}(u^{(j)})$ for each particle
$u^{(j)}$. For the localised covariance we include more weights on particles
in the neighborhood around $u^{(j)}$ and less weight for particles far away.
This type of localisation has been considered for the ensemble Kalman sampler,
where we do not expect the particle system to collapse [53]. The convergence
analysis of (T)EKI under localisation is work in progress.

###  4.3 TEKI as constrained optimization method

Before proving our main convergence result, we need to quantify the behavior
of the preconditioning effect through the sample covariance, i.e. what does it
mean to consider the direction $C(u_{t})\nabla\Phi_{R}(u_{t}^{(j)})$ instead
of the steepest descent direction given through $\nabla\Phi_{R}(u_{t}^{(j)})$.
We first note that once $C(u_{t})$ can be identified as strictly positive
definite operator on the whole parameter space $\mathcal{X}$, the scheme can
be viewed as global optimization method. However, in general it is not true
that $C(u_{t})$ describes a strictly positive definite operator on the
parameter space $\mathcal{X}$ and indeed, $C(u_{t})$ has at most rank
$\min(n_{x},J-1)$ and hence is positive semi-definite in case $J-1<n_{x}$.
While $C(u_{t})$ is only positive semi-definite on the whole space
$\mathcal{X}$, we have verified in Lemma 4.6 that $C(u_{t})$ is a strictly
positive definite operator acting on the subspace $\mathcal{B}$ spanned
through the initial ensemble. We will make use of this property in order to
quantify the TEKI flow (10) as optimizer for the constrained optimization
problem

$\min_{x\in\mathcal{X}}\ \Phi_{R}(x),\quad
x\in\mathcal{B}=u_{0}^{\perp}+\operatorname{span}\\{e_{0}^{(j)}\mid
j=1,\dots,J\\}.$

Since we have assumed the finite dimensional setting
$\mathcal{X}=\mathbb{R}^{n_{x}}$, we reformulate this problem as equality
constrained optimization problem

(19) $\min_{x\in\mathbb{R}^{n_{x}}}\ \Phi_{R}(x),\quad\text{s.t.}\quad
h_{i}(x)=\langle q^{(i)},x\rangle=0,\ i=J+1,\dots,n_{x},$

where $q^{(i)}\perp e_{0}^{(j)}$ and $q^{(i)}\perp u_{0}^{\perp}$ for all
$j=1,\dots,J$ and $i=J+1,\dots,n_{x}$. Since under Assumption 3 $\Phi_{R}$ is
assumed to be strongly convex, there exists a unique global solution
$u_{\ast}$ of (19) [6, Proposition 2.1.1]. By the Lagrange multiplier theorem
[6, Proposition 3.1.1.] there exist
$\kappa_{J+1}^{\ast},\dots,\kappa_{n_{x}}^{\ast}\in\mathbb{R}$

(20) $\nabla\Phi_{R}(u_{\ast})+\sum_{i=J+1}^{n_{x}}\kappa_{i}\nabla
h_{i}(x_{\ast})=0.$

We will make use of this necessary condition in order to derive the following
constrained PL-type inequality.

######  Lemma 4.11.

Suppose Assumption 3 is satisfied and let $x_{\ast}$ be the unique global
minimizer of (19). Then there exists $\nu>0$ such that for all
$x\in{\mathcal{B}}$ it holds true that

$\nu\|\nabla\Phi_{R}(x)\|^{2}\geq\Phi_{R}(x)-\Phi_{R}(x_{\ast}).$

Moreover, let $\mathcal{P}_{\mathcal{B}}$ be the orthogonal projection onto
$\mathcal{B}$. Then it holds true that

${\nu\|\mathcal{P}_{\mathcal{B}}\nabla\Phi_{R}(\mathcal{P}_{\mathcal{B}}x)\|^{2}\geq\Phi_{R}(\mathcal{P}_{\mathcal{B}}x)-\Phi_{R}(x_{\ast})}$

for all $x\in\mathbb{R}^{n_{x}}$.

######  Proof 4.12.

Consider some arbitrary $x\in\mathcal{B}$. By $L$-smoothness and the necessary
optimality condition (20) we have that

$\displaystyle\Phi_{R}(x)$
$\displaystyle\leq\Phi_{R}(x_{\ast})+\langle\nabla\Phi_{R}(x_{\ast}),x-x_{\ast}\rangle+\frac{L}{2}\|x-x_{\ast}\|^{2}$
$\displaystyle=\Phi_{R}(x_{\ast})-\sum_{i=J+1}^{n_{x}}\kappa_{i}\langle\nabla
h_{i}(x_{\ast}),x-x_{\ast}\rangle+\frac{L}{2}\|x-x_{\ast}\|^{2}$
$\displaystyle=\Phi_{R}(x_{\ast})-\sum_{i=J+1}^{n_{x}}\kappa_{i}\langle
q^{(i)},x-x_{\ast}\rangle+\frac{L}{2}\|x-x_{\ast}\|^{2}$
$\displaystyle=\Phi_{R}(x_{\ast})+\frac{L}{2}\|x-x_{\ast}\|^{2},$

where we have used $x,x_{\ast}\perp q^{(i)}$ for all $i=J+1,\dots,n_{x}$. On
the other side with $\mu$-strong convexity and Cauchy-Schwarz we obtain

$\frac{\mu}{2}\|x_{\ast}-x\|\leq\|\nabla\Phi(x)\|$

leading to

$\Phi_{R}(x)-\Phi_{R}(x_{\ast})\leq\frac{L}{2}\|x-x_{\ast}\|^{2}\leq\frac{L}{\mu}\|\nabla\Phi_{R}(x)\|^{2}.$

For the second claim, consider $x\in\mathbb{R}^{n_{x}}$. By $\mu$-strong
convexity it follows

$\langle{\mathcal{P}_{\mathcal{B}}}\nabla\Phi_{R}(x),{\mathcal{P}_{\mathcal{B}}}(x-x_{\ast})\rangle=\langle\nabla\Phi_{R}(x),{\mathcal{P}_{\mathcal{B}}}x-x_{\ast})\rangle\geq\Phi_{R}({\mathcal{P}_{\mathcal{B}}}x)-\Phi_{R}(x_{\ast})+\frac{\mu}{2}\|{\mathcal{P}_{\mathcal{B}}}x-x_{\ast}\|^{2}.$

The assertion follows again by applying Cauchy-Schwarz inequality and
$L$-smoothness.

##  5 Proof of Theorem 3.2

We are now ready to present the proof of our main convergence result.

######  Proof 5.1 (Proof of Theorem 3.2).

Let $\mathcal{B}=u_{0}^{\perp}+\operatorname{span}\\{e_{0}^{(j)}\mid
j=1,\dots,J\\}$ and let $\mathcal{P}_{\mathcal{B}}$ be the orthogonal
projection onto $\mathcal{B}$. Recall that from Proposition 3.1 we have that

$\displaystyle\frac{{\mathrm{d}}\frac{1}{J}\sum_{j=1}^{J}\left(\Phi_{R}(u_{t}^{(j)})\right)}{{\mathrm{d}}t}$
$\displaystyle\leq-\frac{1}{J}\sum_{j=1}^{J}\langle\nabla\Phi_{R}(u_{t}^{(j)}),C(u_{t})\nabla\Phi_{R}(u_{t}^{(j)})\rangle$
$\displaystyle\quad+b_{1}JV_{e}(t)^{3/2}\frac{1}{J}\sum_{j=1}^{J}\left(\frac{1}{2}\|\nabla\Phi_{R}(u_{t}^{(j)})\|^{2}+\Phi_{R}(u_{t}^{(j)})\right)$
$\displaystyle\leq-\frac{1}{J}\sum_{j=1}^{J}\langle\nabla\Phi_{R}(u_{t}^{(j)}),C(u_{t})\nabla\Phi_{R}(u_{t}^{(j)})\rangle$
$\displaystyle\quad+\frac{3}{2}b_{1}J\left(\frac{1}{\frac{2\sigma_{\min}}{J}t+V_{e}(0)^{-1}}\right)^{3/2}\frac{1}{J}\sum_{j=1}^{J}(\Phi_{R}(u_{t}^{(j)})-\Phi_{R}(u_{\ast}))$
$\displaystyle\quad+\tilde{b}_{1}J\left(\frac{1}{\frac{2\sigma_{\min}}{J}t+V_{e}(0)^{-1}}\right)^{3/2}\Phi_{R}(u_{\ast})),$

where we have used strong convexity. Due to the subspace space property in
Proposition 4.3 we obtain
$C(u_{t})=\mathcal{P}_{\mathcal{B}}C(u_{t})\mathcal{P}_{\mathcal{B}}$ and with
Lemma 4.6 it holds true that

$\displaystyle\langle\nabla\Phi_{R}(u_{t}^{(j)}),C(u_{t})\nabla\Phi_{R}(u_{t}^{(j)})\rangle$
$\displaystyle=\langle{\mathcal{P}_{\mathcal{B}}}\nabla\Phi_{R}(u_{t}^{(j)}),C(u_{t}){\mathcal{P}_{\mathcal{B}}}\nabla\Phi_{R}(u_{t}^{(j)})\rangle$
$\displaystyle\geq\frac{1}{mt+\zeta_{0}^{-1}}\|{\mathcal{P}_{\mathcal{B}}}\nabla\Phi_{R}(u_{t}^{(j)})\|^{2},$

where $m>0$ is defined in (18) and $\zeta_{0}$ is defined in (17). Hence, we
define
$V_{\Phi_{R}}(t):=\frac{1}{J}\sum_{j=1}^{J}\Phi_{R}(u_{t}^{(j)})-\Phi_{R}(u_{\ast})$
and apply Lemma 4.11 in order to obtain

$\displaystyle\frac{{\mathrm{d}}V_{\Phi_{R}}(t)}{{\mathrm{d}}t}=\frac{{\mathrm{d}}\frac{1}{J}\sum_{j=1}^{J}\left(\Phi_{R}(u_{t}^{(j)})\right)}{{\mathrm{d}}t}$
$\displaystyle\leq-\frac{1}{\frac{L}{\mu}mt+\zeta_{0}^{-1}}V_{\Phi_{R}}(t)$
$\displaystyle\quad+\tilde{b}_{1}J\left(\frac{1}{\frac{2\sigma_{\min}}{J}t+V_{e}(0)^{-1}}\right)^{3/2}V_{\Phi_{R}}(t)$
$\displaystyle\quad+\tilde{b}_{1}J\left(\frac{1}{\frac{2\sigma_{\min}}{J}t+V_{e}(0)^{-1}}\right)^{3/2}\Phi_{R}(u_{\ast}).$

It follows with $\alpha=\frac{L}{\mu}m>0$ that

$\displaystyle\frac{{\mathrm{d}}V(t)}{{\mathrm{d}}t}\leq-\frac{1}{\alpha
t+a}V(t)+\left(\frac{z_{1}}{z_{2}+z_{3}t}\right)^{3/2}(V(t)+c)$

for some constants $a,z_{1},z_{2},z_{3},c>0$. By Gronwall’s inequality we
obtain

$\displaystyle V(t)$ $\displaystyle\leq
V(0)\left(\int_{0}^{t}\left(\frac{z_{1}}{z_{2}+z_{3}s}\right)^{3/2}\,{\mathrm{d}}s\right)\exp\left(\int_{0}^{t}\left(\frac{z_{1}}{z_{2}+z_{3}s}\right)^{3/2}\,{\mathrm{d}}s-\int_{0}^{t}\frac{1}{\alpha
s+a}\,{\mathrm{d}}s\right)$
$\displaystyle\leq\left(\frac{c_{1}}{t+c_{2}}\right)^{\frac{1}{\alpha}},$

for some constants $c_{1},c_{2}>0$.

######  Remark 5.2.

We note that under the relaxed assumption $\|{\mathrm{Res}}(x,y)\|\leq
b_{\mathrm{res}}\|x-y\|^{p}$, for $p\geq 2$ such that $p/2\in\mathbb{N}$, we
obtain a similar bound for the gradient approximation as in Lemma 4.8 given by

$\displaystyle\|C^{u,G}(u)\Gamma^{-1}(G(u^{(j)})-y)-C(u)\nabla\Phi(u^{(j)})\|$
$\displaystyle\leq
b_{\mathrm{res}}\|\Gamma^{-1/2}\|(8+\frac{2}{J})J(4J)^{p/2}\|G(u^{(j)})-y\|_{\Gamma}\left(\frac{1}{J}\sum_{k=1}^{J}\|u^{(k)}-\bar{u}\|^{2}\right)^{1/2+p/2}.$

Here, we have applied Binomial theorem to imply the inequality

$\sum_{j=1}^{J}a_{j}^{p}\leq\left(\sum_{j=1}^{J}a_{j}^{2}\right)^{p/2},\quad
a_{j}\geq 0.$

Since $p\geq 2$ the integral
$\int_{0}^{t}\left(\frac{z_{1}}{z_{2}+z_{3}s}\right)^{1/2+p/2}\,{\mathrm{d}}s$
in the above proof remains finite and hence the asymptotic convergence result
in Theorem 3.2 still holds.

##  6 Stabilization via covariance inflation: A generalized ensemble square
root filter for inverse problems

The incorporation of covariance inflation is a data assimilation tool in order
to stabilize convergence and often used in practice [2, 1, 59]. In the
following we present a form of covariance inflation for the TEKI flow for
which it turns out to be a generalized formulation of the ensemble square root
filter (ESRF) [52] applied to inverse problems.

Through our incorporation of covariance inflation we obtain in approximation a
weighted gradient flow structure between the current particle position and the
position of the ensemble mean. Roughly speaking, the TEKI flow can then be
viewed as approximation of the equation:

(21)
$\frac{\mathrm{d}u_{t}^{(j)}}{{\mathrm{d}}t}=-C(u_{t})\left((1-\rho)\nabla\Phi_{R}(u_{t}^{(j)})+\rho\nabla\Phi_{R}(\bar{u}_{t})\right),$

where $\rho\in[0,1)$ describes the weight of the gradients of the particles
and the mean. Based on [5, 4] the continuous time formulation of the ESRF
applied to inverse problems has been analysed in [12] for linear forward
models and the special case $\rho=1/2$. In the following, we will incorporate
the covariance inflation and theoretical quantify the approximate gradient
flow structure. Finally, we will observe an acceleration in convergence
towards the optimal solution restricted to the subspace spanned through the
initial ensemble spread.

?figurename? 1: Illustration of the effect through the incorporation of the
covariance inflation. Each particle is getting moved into its own direction as
well as into a joint direction.

###  6.1 Incorporation of covariance inflation

The covariance inflation is incorporated through the following modified TEKI
flow

(22)
$\begin{split}\frac{{\mathrm{d}}u_{t}^{(j)}}{{\mathrm{d}}t}&=C^{u,G}(u_{t})\Gamma^{-1}(y-G(u_{t}^{(j)}))-C(u_{t})C_{0}^{-1}u_{t}^{(j)}\\\
&\quad+\rho C^{u,G}(u_{t})\Gamma^{-1}(G(u_{t}^{(j)})-\bar{G})+\rho
C(u_{t})C_{0}^{-1}(u_{t}^{(j)}-\bar{u}).\end{split}$

The idea behind this incorporation comes from the fact that the inflating term
does not directly change the dynamical behavior of the particles mean, since
the evolution

(23)
$\frac{{\mathrm{d}}\bar{u}_{t}}{{\mathrm{d}}t}=C^{u,G}(u_{t})\Gamma^{-1}(y-\bar{G}_{t})-C(u_{t})C_{0}^{-1}\bar{u}_{t}$

coincides with the unmodified dynamic (i.e. $\rho=0$). However the inflating
term decreases the force of collapsing the particle system and hence leads to
a slower ensemble collapse without changing the collapse rate $1/t$ but
decreasing the constant. For the interpretation from the optimization point of
view we rewrite (22) as

(24)
$\begin{split}\frac{{\mathrm{d}}u_{t}^{(j)}}{{\mathrm{d}}t}&=(1-\rho)\left\\{C^{u,G}(u_{t})\Gamma^{-1}(y-G(u_{t}^{(j)}))-C(u_{t})C_{0}^{-1}u_{t}^{(j)}\right\\}\\\
&\quad+\rho\left\\{C^{u,G}(u_{t})\Gamma^{-1}(y-\bar{G}_{t})-C(u_{t})C_{0}^{-1}\bar{u}_{t}\right\\}.\end{split}$

motivating the modified TEKI flow as approximation of (21). The movement of
the particle system through (24) is described through two forces. The first
force is given by

$C^{u,G}(u_{t})\Gamma^{-1}(y-G(u_{t}^{(j)}))-C(u_{t})C_{0}^{-1}u_{t}^{(j)}$

moving each particle in its own direction and leading to a collapse of the
ensemble and hence to a better gradient approximation as we have already seen
in Lemma 4.8. The second force

$C^{u,G}(u_{t})\Gamma^{-1}(y-\bar{G}_{t})-C(u_{t})C_{0}^{-1}\bar{u}_{t}$

moves the whole particle system into the same direction without changing the
spread and hence does not improve the gradient approximation. However, it
accelerates the overall convergence of the optimization process since it slows
down the degeneration of the preconditioning through the sample covariance
matrix. We will start the theoretical analysis by considering the ensemble
spread and its collapse behavior. We first note, that unique existence of
solutions of (22) can be verified and the subspace property Proposition 4.3 is
obviously satisfied.

######  Proposition 6.1.

Let $(u_{0}^{(j)})_{j=1,\dots,J}$ be the linearly independent initial ensemble
and define the subspace

$\mathcal{B}=u_{0}^{\perp}+\operatorname{span}\\{e_{0}^{(j)}\mid
j=1,\dots,J\\},$

where $u_{0}^{\perp}=\bar{u}_{0}-\mathcal{P}_{E}\bar{u}_{0}$, where
$\mathcal{P}_{E}=E(E^{\top}E)^{-1}E^{T}$ with
$E=((e_{0}^{(1)})^{\top},\dots,(e_{0}^{(J)})^{\top})\in\mathbb{R}^{n_{x}\times
J}$ is the orthogonal projection on the subspace
$\operatorname{span}\\{e_{0}^{(j)}\mid j=1,\dots,J\\}$. Then for all $T\geq 0$
and $\rho\in[0,T)$ the unique solution $(u_{t}^{(j)},t\in[0,T])$ of the flow
(22) exists and remains in $\mathcal{B}$, i.e. $u_{t}^{(j)}\in\mathcal{B}$ for
all $t\in[0,T]$ and all $j=1,\dots,J$.

The following result states, that the ensemble collapse occurs with reduced
rate through the factor $1-\rho\in(0,1]$.

######  Lemma 6.2.

Let $\rho\in[0,1)$ and $(u_{t}^{(j)},t\geq 0)_{j=1,\dots,J}$ be the solution
of (22) initialized by a linearly independent ensemble
$(u_{0}^{(j)})_{j=1,\dots,J}$. Then the mapping
$V_{e}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ with
$t\mapsto\frac{1}{J}\sum_{j=1}^{J}\|e_{t}\|^{2}$ is bounded by

$V_{e}(t)\leq\frac{1}{\frac{2(1-\rho)\sigma_{\min}}{J}t+V_{e}(0)^{-1}},$

where $\sigma_{\min}$ denotes the smallest eigenvalue of $C_{0}^{-1}$.

######  Proof 6.3.

The evolution of $V_{e}(t)$ is given by

$\displaystyle\frac{{\mathrm{d}}\frac{1}{J}\sum_{j=1}^{J}\|u_{t}^{(j)}-\bar{u}_{t}\|^{2}}{{\mathrm{d}}t}$
$\displaystyle=-(1-\rho)2\|C^{u,G}(u_{t})\Gamma^{-1/2}\|_{\mathrm{HS}}$
$\displaystyle\quad-\frac{(1-\rho)2}{J^{2}}\sum_{j,k=1}^{J}\langle
u_{t}^{(j)}-\bar{u}_{t},u_{t}^{(k)}-\bar{u}_{t}\rangle\langle
u_{t}^{(k)}-\bar{u}_{t},C_{0}^{-1}(u_{t}^{(j)}-\bar{u}_{t})\rangle$
$\displaystyle\leq-\frac{2(1-\rho)\sigma_{\min}}{J}\left(\frac{1}{J}\sum_{j=1}^{J}\|u_{t}^{(j)}-\bar{u}_{t}\|^{2}\right)^{2},$

and the assertion follows again by Lyapunov-type argument.

Similarly, we can derive the lower bound of the ensemble collapse.

######  Lemma 6.4.

Let $\rho\in[0,1)$ and $(u_{t}^{(j)},t\geq 0)_{j=1,\dots,J}$ be the solution
of (22) initialized by a linearly independent ensemble
$(u_{0}^{(j)})_{j=1,\dots,J}$ such that

$\zeta_{0}=\min_{z\in{\mathcal{B}},\ \|z\|=1}\ \langle z,C(u_{0})z\rangle>0.$

For each $z\in{\mathcal{B}}$ with $\|z\|=1$ it holds true that

$\langle z,C(u_{t})z\rangle\geq\frac{1}{(1-\rho)mt+\zeta_{0}^{-1}},$

where $m=2(c_{\mathrm{lip}}^{2}\lambda_{\max}V_{e}(0)+\sigma_{\max})>0$
depends on the smallest and largest eigenvalue $\sigma_{\max}$ and
$\sigma_{\min}$ of $C_{0}^{-1}$, the largest eigenvalue $\lambda_{\max}$ of
$\Gamma^{-1}$ and the Lipschitz constant $c_{\mathrm{lip}}$ of $G$.

######  Proof 6.5.

The time evolution of the sample covariance is given by

$\displaystyle\frac{{\mathrm{d}}C(u_{t})}{{\mathrm{d}}t}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\frac{{\mathrm{d}}e_{t}^{(j)}}{{\mathrm{d}}t}\left(e_{t}^{(j)}\right)^{\top}+\frac{1}{J}\sum_{j=1}^{J}e_{t}^{(j)}\left(\frac{{\mathrm{d}}e_{t}^{(j)}}{{\mathrm{d}}t}\right)^{\top}$
$\displaystyle=-2(1-\rho)C^{u,G}(u_{t})\Gamma^{-1}C^{G,u}(u_{t})-2(1-\rho)C(u_{t})C_{0}^{-1}C(u_{t})$

and similarly, for $z\in{\mathcal{B}}$ with $\|z\|=1$ we obtain

$\displaystyle\langle z,\frac{{\mathrm{d}}C(u_{t})}{{\mathrm{d}}t}z\rangle$
$\displaystyle=-2(1-\rho)\|C^{G,u}(u_{t})z\|_{\Gamma}^{2}-2(1-\rho)\|C(u_{t})z\|_{C_{0}}^{2}$

and therefore

$\frac{{\mathrm{d}}\langle
z,C(u_{t})z\rangle}{{\mathrm{d}}t}\geq-2(1-\rho)(c_{\mathrm{lip}}^{2}\lambda_{\max}V_{e}(0)+\sigma_{\max})\|C(u_{t})z\|^{2}.$

We consider again the smallest eigenvalue $\zeta_{t}\geq 0$ of $C(u_{t})$
restricted to $\mathcal{B}$ with unit-norm eigenvector
$\varphi(t)\in{\mathcal{B}}$. For the time-evolution of $\zeta_{t}$ we obtain

$\frac{{\mathrm{d}}}{{\mathrm{d}}t}\zeta_{t}\geq-2(1-\rho)(c_{\mathrm{lip}}^{2}\lambda_{\max}V_{e}(0)+\sigma_{\max})\zeta_{t}^{2},$

such that the assertion follows.

Our convergence result will be based on the mean of the particle system.
Therefore we will need the following extension of Lemma 4.8.

######  Lemma 6.6.

Suppose that Assumption 3 is satisfied. Given a particle system
$(u^{(j)})_{j=1,\dots,J}$ it holds true that

$\|C^{u,G}(u)\Gamma^{-1}(\bar{G}-y)-C(u)\nabla\Phi(\bar{u})\|\leq
b_{1}J\sqrt{\Phi(\bar{u})}\left(\frac{1}{J}\sum_{k=1}^{J}\|u^{(k)}-\bar{u}\|^{2}\right)^{3/2},$

for some $b_{1}$ independent of $J$.

######  Proof 6.7.

We apply triangular inequality

$\displaystyle\|C^{u,G}(u)\Gamma^{-1}(\bar{G}-y)-C(u)\nabla\Phi(\bar{u})\|$
$\displaystyle\leq\|C^{u,G}(u)\Gamma^{-1}(\bar{G}-G(\bar{u}))\|$
$\displaystyle\quad+\|C^{u,G}(u)\Gamma^{-1}(G(\bar{u})-y)-C(u)\nabla\Phi(\bar{u})\|.$

Using Lipschitz continuity and the inequality
$\|C^{u,G}\|_{\mathrm{HS}}\leq\sqrt{\|C^{u,u}\|_{\mathrm{HS}}\|C^{G,G}\|_{\mathrm{HS}}}$
we can bound the first term

$\|C^{u,G}(u)\Gamma^{-1}(\bar{G}-G(\bar{u}))\|\leq\|\Gamma^{-1}\|c_{\mathrm{lip}}^{2}\left(\frac{1}{J}\sum_{k=1}^{J}\|u^{(k)}-\bar{u}\|^{2}\right)^{3/2}.$

The second term can be bounded similarly to Lemma 4.8 and the assertion
follows.

######  Theorem 6.8.

Suppose Assumption 3 and Assumption 3 are satisfied, let $(u_{t}^{(j)},t\geq
0)_{j=1,\dots,J}$ solve (22) with $\rho\in[0,1)$ and linearly independent
initial ensemble $(u_{0}^{(j)})_{j=1,\dots,J}$ and let
$u_{\ast}\in\mathcal{B}$ be the unique global minimizer of (14). Then there
exist $c_{1},c_{2}>0$ such that

(25)
$\Phi_{R}(\bar{u}_{t})-\Phi_{R}(u_{\ast})\leq\left(\frac{c_{1}}{t+c_{2}}\right)^{\frac{1}{\alpha}},$

where
$0<\alpha<(1-\rho)\frac{L}{\mu}(\sigma_{\max}+c_{\mathrm{lip}}\lambda_{\max}\|C(u_{0})\|_{\mathrm{HS}})$.

######  Proof 6.9.

With the above presented results in Lemma 6.2, Lemma 6.4 and Lemma 6.6, the
proof is similar to the proof of Theorem 3.2.

Through the incorporation of covariance inflation we have reduced the ensemble
collapse by decreasing the constant of the degeneration rate through the
constant inflation parameter $(1-\rho)\in(0,1]$. Resulting from this slow down
of the ensemble collapse the convergence rate of TEKI as optimization method
can be accelerated by the factor $(1-\rho)$.

##  7 Numerical experiment

In our numerical experiment we consider a nonlinear inverse problem based on
the one-dimensional elliptic boundary-value problem in the form of

(26) $\begin{array}[]{r@{}l}-\frac{{\mathrm{d}}}{{\rm
d}s}\left(\exp(u(s))\frac{{\mathrm{d}}}{{\mathrm{d}}s}p(s)\right)&=1,\qquad
s\in[0,1],\\\ p(s)&=0,\qquad s\in\\{0,1\\}.\end{array}$

Given discrete noisy observations $y\in\mathbb{R}^{K}$ of the solution $p\in
H_{0}^{1}([0,1];\mathbb{R})$ we consider the task of recovering the unknown
parameter $u\in L^{\infty}([0,1])$. To be more precise, we consider
$y=(\mathcal{O}\circ S)(u(s))+\eta$, where $\eta\sim\mathcal{N}(0,\Gamma)$ is
additive Gaussian noise, $S:\ L^{\infty}([0,1])\to
H_{0}^{1}([0,1];\mathbb{R})$ denotes the solution operator of (26) and
$\mathcal{O}:\ H_{0}^{1}([0,1];\mathbb{R})\to\mathbb{R}^{K}$ denotes the
observation operator, which provides function values at $n_{y}$ equidistant
observation points in $[0,1]$, i.e. $z(\cdot)\in
H_{0}^{1}([0,1];\mathbb{R})\mapsto\mathcal{O}(z(\cdot))=(z(s_{1}),\dots,z(s_{K}))^{\top}$,
$s_{i}=\frac{i}{K},\ i=1,\dots,K$.

For our numerical results we replace $S$ by an numerical solution operator for
(26) on the grid $\mathcal{D}\subset[0,1]$ with mesh size $h=2^{-8}$ and
restrict $u(\cdot)$ to the computational grid $s_{l}=l\,h$,
$l=1,\ldots,2^{8}-1$, i.e. we consider the finite dimensional parameter space
$\mathcal{X}=\mathbb{R}^{n_{x}}$ with $n_{x}=2^{8}$. Our observations are
taken at $K=2^{5}-1$ observation points and the noise covariance is assumed to
be $\Gamma=0.01\cdot{\mathrm{Id}}\in\mathbb{R}^{K\times K}$.

We set a Gaussian prior $\mathcal{N}(0,C_{0})$ on the unknown parameter with
covariance operator $C_{0}=\beta(-\Delta)^{-1}$, $\beta=10$. The underlying
ground truth is constructed as realization
$u^{\dagger}\sim\mathcal{N}(0,C_{0})$ and shown in Figure 2 with corresponding
taken observations.

?figurename? 2: Underlying ground truth of the unknown parameter $u$ (left)
and the corresponding (noisy) observations (right).

###  7.1 Fixed regularization

In our first experiment, we consider a fixed regularization and two types of
initialization for the implementation of the TEKI flow:

  * •

Basis initialization: we choose the initial ensemble of particles based on
(ordered) eigenvalues and eigenfunctions
$\\{\lambda_{j},z_{j}\\}_{j\in\mathbb{N}}$ of the covariance operator $C_{0}$.

  * •

Random initialization: we choose the initial ensemble of particles based on
i.i.d. realizations of the Gaussian prior $\mathcal{N}(0,C_{0})$.

The TEKI flow with covariance inflation and regularization parameter
$\kappa=0.01^{2}$

(27)
$\begin{split}\frac{{\mathrm{d}}u_{t}^{(j)}}{{\mathrm{d}}t}&=C^{u,G}(u_{t})\Gamma^{-1}(y-G(u_{t}^{(j)}))-C(u_{t})\kappa
C_{0}^{-1}u_{t}^{(j)}\\\ &\quad+\rho
C^{u,G}(u_{t})\Gamma^{-1}(G(u_{t}^{(j)})-\bar{G})+\rho C(u_{t})\kappa
C_{0}^{-1}(u_{t}^{(j)}-\bar{u})\end{split}$

is tested for different choices of inflation scaling
$\rho\in\\{0,0.25,0.5,0.8\\}$ and solved through MATLAB via the function ode45
until final time $T=10^{7}$ with the aim of minimizing

$\Phi_{R}(u)=\frac{1}{2}\|G(u)-y\|_{\Gamma}^{2}+\frac{\kappa}{2}\|u\|_{C_{0}}^{2},\quad
u\in\mathbb{R}^{n_{x}}.$

The reference solution $u_{\ast}^{J}$ of the constrained optimization problem
(19) depending of the ensemble size $J\in\\{5,20,50\\}$ is approximated with
the MATLAB function fmincon. We note that $u_{\ast}^{J}$ depends on the
initial ensemble and differs for basis and random initialization.

In Figure 3 we show that the ensemble collapses with rate $1/t$ and illustrate
the effect of the covariance inflation slowing down the collapse by a
constant. This observation verifies our expectation on the behavior from the
results in Lemma 6.2 and Lemma 6.4. Here, we observe no significant difference
between basis initialization and random initialization. We demonstrate our
main convergence result presented in Theorem 3.2 and Theorem 6.8 through
Figure 4-6, where we observe that the lossfunctions evaluated at the particle
mean converge towards the expected minimum restricted to the subspace spanned
by the initial ensemble spread. We firstly observe, that increasing the
ensemble size, the resulting minimum improves since the initial subspace
increases. Secondly, we observe that the resulting estimates for basis
function initialization lead to lower values in the lossfunction compared to
random initialization. This behavior is expected since the basis functions are
weighted through the ordered eigenvalues. Furthermore, in Figure 7 we
illustrate the improved convergence through increasing variance inflation
parameter $\rho$. Finally, Figure 8 shows the resulting estimators of the
unknown parameter $u$ where we observe smooth estimates for basis function
initialization compared to the estimates resulting from random initialization.

?figurename? 3: Particle spread $\frac{1}{J}\sum_{j=1}^{J}\|e_{t}^{(j)}\|^{2}$
of the TEKI flow for basis function initialization (left) and random
initialization (right) with ensemble size $J=50$.

?figurename? 4: Lossfunction $\Phi_{R}$ evaluated at the particle mean
$\bar{u}_{t}$ of the TEKI flow for basis function initialization (left) and
random initialization (right) with ensemble size $J=5$.

?figurename? 5: Lossfunction $\Phi_{R}$ evaluated at the particle mean
$\bar{u}_{t}$ of the TEKI flow for basis function initialization (left) and
random initialization (right) with ensemble size $J=20$.

?figurename? 6: Lossfunction $\Phi_{R}$ evaluated at the particle mean
$\bar{u}_{t}$ of the TEKI flow for basis function initialization (left) and
random initialization (right) with ensemble size $J=50$.

?figurename? 7: Difference of the lossfunction $\Phi_{R}$ evaluated at the
particle mean $\bar{u}_{t}$ the TEKI flow for basis function initialization
(left) and random initialization (right) with ensemble size $J=50$ and the
optimal solution restricted to the initial subspace $u_{\ast}^{J}$.

?figurename? 8: Resulting estimates of the unknown parameter $u$ for different
ensemble sizes $J\in\\{5,20,50\\}$ of the TEKI flow initialized by basis
functions (left) and random (right).

###  7.2 Adaptive regularization

While in the previous experiment we have tested TEKI with fixed regularization
parameter in order to verify our theoretical findings, we are providing a
second experiment motivated from a more practical point of view. Here, we move
away from the optimistic scenario, where the reference $u^{\dagger}$ is chosen
from the correct prior, instead we assume
$u^{\dagger}\sim\mathcal{N}(0,{C}_{0}^{\dagger})$ with
${C}_{0}^{\dagger}:=\tilde{\beta}(-\Delta)^{-\alpha}$, $\tilde{\beta}=1$,
$\alpha=2$. Hence, the incorporated smoothness information through the
regularization matrix $C_{0}$ may be misleading. We therefore propose to apply
an adaptive regularization schemes introduced [62] for the stochastic
formulation of TEKI. The idea of choosing a ”good” type of Tikhonov
regularization is motivated from the hierarchical approach for Bayesian
inverse problems. In the classical setting we are interested in the posterior
distribution

$u\mid y\sim\mu({\mathrm{d}}x)\propto\exp(-\Phi(x;y))\,\mu_{0}({\mathrm{d}}x)$

whereas in the hierarchical setting we assume that the prior is depending on
some unknown hyperparameter $\Theta\in\mathbb{R}^{n_{\theta}}$ such that we
are interested in the posterior

$(u,\Theta)\mid
y\sim\mu({\mathrm{d}}(x,\theta))\propto\exp(-\Phi(x;y))\mu_{0}({\mathrm{d}}x,\theta)q_{0}({\mathrm{d}}\theta)$

assuming a prior $\Theta\sim q_{0}({\mathrm{d}}\theta)$. We consider the case
where $\mu_{0}$ is given as Gaussian prior with pdf

$\mu_{0}(x,\theta)=\frac{1}{\sqrt{2\pi\det(C_{0}(\theta))}}\exp(-\|x-m_{0}(\theta)\|_{C_{0}(\theta)}^{2},$

where the hyperparameter $\theta$ parametrizes the prior mean $m_{0}(\theta)$
and prior covariance $C_{0}(\theta)$, and we assume an uniform prior on
$\Theta\sim\mathcal{U}((0,M)^{n_{\theta}})$. When computing the joint MAP of
$(u,\Theta)$ the normalizing constant plays an important role such that the
optimization problem reads as

$\min_{x\in\mathbb{R}^{n_{x}},\theta\in(0,M)^{n_{\theta}}}\
\frac{1}{2}\|G(x)-y\|_{\Gamma}^{2}+\frac{1}{2}\|u-m_{0}(\theta)\|_{C_{0}(\theta)}+\frac{1}{2}\log(\det(C_{0}(\theta))).$

In [62] the authors propose to apply a two-level update scheme, where on the
first level one applies TEKI to optimize w.r.t. $x$ given $\theta$, while in
the second level one applies Gradient descent w.r.t. $\theta$ given $x$. In
continuous time the resulting adaptive TEKI flow then reads as

$\displaystyle\frac{{\mathrm{d}}u_{t}^{(j)}}{{\mathrm{d}}t}$
$\displaystyle=C^{u,G}(u_{t})\Gamma^{-1}(y-G(u_{t}^{(j)}))-C(u_{t})C_{0}(\theta_{t}^{(j)})^{-1}u_{t}^{(j)},$
$\displaystyle\frac{{\mathrm{d}}\theta_{t}^{(j)}}{{\mathrm{d}}\theta}$
$\displaystyle=-\nabla_{\theta}\left(\frac{1}{2}\|u_{t}^{(j)}-m_{0}(\theta)\|_{C_{0}(\theta^{(j)})}+\frac{1}{2}\log(\det(C_{0}(\theta^{(j)})))\right).$

Alternatively, we consider the second level as Gradient flow w.r.t.
$\theta_{t}$ given the particle mean $\bar{u}_{t}$, i.e.

(28)
$\begin{split}\frac{{\mathrm{d}}u_{t}^{(j)}}{{\mathrm{d}}t}&=C^{u,G}(u_{t})\Gamma^{-1}(y-G(u_{t}^{(j)}))-C(u_{t})C_{0}(\theta_{t}^{(j)})^{-1}u_{t}^{(j)},\\\
\frac{{\mathrm{d}}\theta_{t}}{{\mathrm{d}}\theta}&=-\nabla_{\theta}\left(\frac{1}{2}\|\bar{u}_{t}-m_{0}(\theta)\|_{C_{0}(\theta)}+\frac{1}{2}\log(\det(C_{0}(\theta)))\right).\end{split}$

We note that this type of adaptive choice of regularization remains
derivative-free w.r.t. the forward model $G$. In the following, we consider
$m_{0}(\theta)=0$ and a parametrization of $C_{0}(\theta)$ of the form

$C_{0}(\theta)=V\begin{pmatrix}1/\theta_{1}&&\\\ &\ddots&\\\
&&1/\theta_{n_{x}}\end{pmatrix}V^{\top},$

where $V\in\mathbb{R}^{n_{x}\times n_{x}}$ is an orthonormal matrix built
through the eigenfunctions of $C_{{\mathrm{fix}}}=(-\Delta)^{-1}$. We compare
the adaptive TEKI flow (28) with the fixed TEKI flow (27) (with $\rho=0$) with
$\kappa\in\\{0.001,1\\}$ and fixed regularization matrix $C_{{\mathrm{fix}}}$.
As mentioned above the underlying ground truth has been drawn from
$u^{\dagger}\sim\mathcal{N}(0,{C}_{0}^{\dagger})$. We initialize all of the
three methods with the same initial ensemble
$u_{0}^{(j)}\sim\mathcal{N}(0,C_{{\mathrm{fix}}}),\ J=100$ and choose
$\theta_{0}$ such that $C_{0}(\theta_{0})=1\cdot C_{{\mathrm{fix}}}$.

?figurename? 9: Resulting estimates of the unknown parameter $u^{\dagger}$ for
different types of regularization of the TEKI flow (left) and corresponding
residuals over time (right).

In Figure 9 we plot the corresponding parameter estimation for all of the
three schemes as well as the residuals $\|\bar{u}_{t}-u^{\dagger}\|^{2}$ for
parameter reconstruction, where $\bar{u}_{t}$ denotes the particle mean of the
applied (fixed/adaptive) TEKI flow. As expected, for too large choices of
$\kappa=1$ we observe over-smoothing in the parameter reconstruction, whereas
for too low choice of $\kappa=0.001$ we observe over-fitting issues. The
adaptive regularization scheme helps to significantly improve the performance
of the resulting reconstruction of the underlying ground truth $u^{\dagger}$.
Finally, we refer to [62] for more details on adaptive regularization schemes
for TEKI.

##  8 Conclusion

In this work, we have quantified the EKI as derivative-free optimization
method restricted to the finite dimensional subspace spanned through the
initial ensemble spread. Our results are based on the deterministic
continuous-time formulation represented by an ODE and the incorporation of
Tikhonov regularization. The ensemble collapse has been quantified from above
and below such that we were able to bound the approximation of gradients
uniform in time and avoiding too fast degeneration of the preconditioner. We
have further improved the convergence result through covariance inflation
without breaking the subspace property. The considered covariance inflation
can be viewed as generalization of the ESRF and gives an intuitive connection
between the EKI and ESRF. In our numerical experiment we have illustrated the
application of EKI as optimization method restricted to the initial subspace
and specified the role of choosing the initial ensemble.

The next step for future work includes the extension of the presented results
to the stochastic formulation of EKI represented as stochastic differential
equation. The key challenge in this scenario will be the quantification of the
ensemble collapse, since the resulting lower and upper bounds on the
eigenvalues will be path-depending. In addition, a detailed complexity
analysis of EKI needs to be done in order to quantify the advantage of
applying EKI as derivative-free scheme compared to other commonly used
optimization methods including derivatives such as gradient descent or
(Quasi-)Newton method.

Acknowledgement. The author is very grateful for helpful discussions with
Claudia Schillings and Jakob Zech which had an significant impact on the
contribution of this work. Moreover, the author thanks Neil Chada and Claudia
Schillings for carefully proofreading this manuscript.

##  ?refname?

  * [1] J. L. Anderson, An adaptive covariance inflation error correction algorithm for ensemble filters, Tellus A: Dynamic Meteorology and Oceanography, 59 (2007), pp. 210–224.
  * [2] , Spatially and temporally varying adaptive covariance inflation for ensemble filters, Tellus A: Dynamic Meteorology and Oceanography, 61 (2009), pp. 72–83.
  * [3] M. Benning and M. Burger, Modern regularization methods for inverse problems, Acta Numerica, 27 (2018), p. 1–111.
  * [4] K. Bergemann and S. Reich, A localization technique for ensemble Kalman filters, Quarterly Journal of the Royal Meteorological Society, 136 (2010), pp. 701–707.
  * [5] , A mollified ensemble Kalman filter, Quarterly Journal of the Royal Meteorological Society, 136 (2010), pp. 1636–1643.
  * [6] D. P. Bertsekas, Nonlinear programming, Athena Scientific, 2nd ed., Sept. 2008.
  * [7] D. Blömker, C. Schillings, P. Wacker, and S. Weissmann, Continuous time limit of the stochastic ensemble Kalman inversion: Strong convergence analysis, ArXiv e-prints, Preprint arXiv:2107.14508 (2021).
  * [8] D. Blömker, C. Schillings, and P. Wacker, A strongly convergent numerical scheme from ensemble Kalman inversion, SIAM Journal on Numerical Analysis, 56 (2018), pp. 2537–2562.
  * [9] D. Blömker, C. Schillings, P. Wacker, and S. Weissmann, Well posedness and convergence analysis of the ensemble Kalman inversion, Inverse Problems, 35 (2019), p. 085007.
  * [10] L. Bungert and P. Wacker, Long-time behaviour and spectral decomposition of the linear ensemble Kalman inversion in parameter space, ArXiv e-prints, Preprint arXiv:2104.13281 (2021).
  * [11] N. K. Chada, A. Jasra, and F. Yu, Multilevel ensemble Kalman-bucy filters, ArXiv e-prints, Preprint arXiv:2011.04342 (2021).
  * [12] N. K. Chada, C. Schillings, and S. Weissmann, On the incorporation of box-constraints for ensemble Kalman inversion, Foundations of Data Science, 1 (2019), pp. 433–456.
  * [13] N. K. Chada, A. M. Stuart, and X. T. Tong, Tikhonov regularization within ensemble Kalman inversion, SIAM Journal on Numerical Analysis, 58 (2020), pp. 1263–1294.
  * [14] N. K. Chada and X. T. Tong, Convergence acceleration of ensemble Kalman inversion in nonlinear settings, Math. Comp., 91 (2022), pp. 1247–1280.
  * [15] A. Chambolle, V. Caselles, D. Cremers, M. Novaga, and T. Pock, An Introduction to Total Variation for Image Analysis, De Gruyter, Berlin, Boston, 16 Jul. 2010.
  * [16] Y. Chen and D. S. Oliver, Ensemble randomized maximum likelihood method as an iterative ensemble smoother, Mathematical Geosciences, 44 (2012), pp. 1–26.
  * [17] A. Chernov, H. Hoel, K. J. H. Law, F. Nobile, and R. Tempone, Multilevel ensemble Kalman filtering for spatio-temporal processes, Numerische Mathematik, 147 (2021), pp. 71–125.
  * [18] J. de Wiljes, S. Reich, and W. Stannat, Long-time stability and accuracy of the ensemble Kalman–bucy filter for fully observed processes and small measurement noise, SIAM Journal on Applied Dynamical Systems, 17 (2018), pp. 1152–1181.
  * [19] P. Del Moral and J. Tugaut, On the stability and the uniform propagation of chaos properties of ensemble Kalman bucy filters, The Annals of Applied Probability, 28 (2018), pp. 790–850.
  * [20] Z. Ding and Q. Li, Ensemble Kalman inversion: mean-field limit and convergence analysis, Statistics and Computing, 31 (2021), p. 9.
  * [21] , Ensemble Kalman sampler: Mean-field limit and convergence analysis, SIAM Journal on Mathematical Analysis, 53 (2021), pp. 1546–1578.
  * [22] Z. Ding, Q. Li, and J. Lu, Ensemble Kalman inversion for nonlinear problems: Weights, consistency, and variance bounds, Foundations of Data Science, 0 (2020), pp. –.
  * [23] A. A. Emerick and A. C. Reynolds, Ensemble smoother with multiple data assimilation, Computers & Geosciences, 55 (2013), pp. 3–15.
  * [24] H. Engl, M. Hanke, and G. Neubauer, Regularization of Inverse Problems, Mathematics and Its Applications, Springer Netherlands, 1996.
  * [25] H. W. Engl, K. Kunisch, and A. Neubauer, Convergence rates for Tikhonov regularisation of non-linear ill-posed problems, Inverse Problems, 5 (1989), pp. 523–540.
  * [26] O. G. Ernst, B. Sprungk, and H.-J. Starkloff, Analysis of the ensemble and polynomial chaos Kalman filters in Bayesian inverse problems, SIAM/ASA Journal on Uncertainty Quantification, 3 (2015), pp. 823–851.
  * [27] G. Evensen, The Ensemble Kalman filter: theoretical formulation and practical implementation, Ocean Dynamics, 53 (2003), pp. 343–367.
  * [28] G. Evensen, Data Assimilation: The Ensemble Kalman Filter, Springer Berlin Heidelberg, 2009.
  * [29] A. Garbuno-Inigo, F. Hoffmann, W. Li, and A. M. Stuart, Interacting Langevin diffusions: Gradient structure and ensemble Kalman sampler, SIAM Journal on Applied Dynamical Systems, 19 (2020), pp. 412–441.
  * [30] A. Garbuno-Inigo, N. Nüsken, and S. Reich, Affine invariant interacting Langevin dynamics for Bayesian inference, SIAM Journal on Applied Dynamical Systems, 19 (2020), pp. 1633–1658.
  * [31] P. A. Guth, C. Schillings, and S. Weissmann, 14 ensemble kalman filter for neural network-based one-shot inversion, In: Optimization and Control for Partial Differential Equations: Uncertainty quantification, open and closed-loop control, and shape optimization, edited by R. Herzog and M. Heinkenschloss and D. Kalise and G. Stadler and E. Trélat, Berlin, Boston: De Gruyter (2022), pp. 393–418.
  * [32] M. Herty and G. Visconti, Kinetic methods for inverse problems, Kinetic & Related Models, 12 (2019), pp. 1109–1130.
  * [33] H. Hoel, K. Law, and R. Tempone, Multilevel ensemble Kalman filtering, SIAM Journal on Numerical Analysis, 54 (2016), pp. 1813–1839.
  * [34] H. Hoel, G. Shaimerdenova, and R. Tempone, Multilevel ensemble Kalman filtering based on a sample average of independent EnKF estimators, Foundations of Data Science, 2 (2020), pp. 351–390.
  * [35] M. Iglesias and Y. Yang, Adaptive regularisation for ensemble Kalman inversion, Inverse Problems, 37 (2021), p. 025008.
  * [36] M. A. Iglesias, Iterative regularization for ensemble data assimilation in reservoir models, Computational Geosciences, 19 (2015), pp. 177–212.
  * [37] , A regularizing iterative ensemble Kalman method for PDE-constrained inverse problems, Inverse Problems, 32 (2016), p. 025002.
  * [38] M. A. Iglesias, K. J. H. Law, and A. M. Stuart, Ensemble Kalman methods for inverse problems, Inverse Problems, 29 (2013), p. 045001.
  * [39] D. Kelly, K. Law, and A. M. Stuart, Well-posedness and accuracy of the ensemble Kalman filter in discrete and continuous time, Nonlinearity, 27 (2014), p. 2579.
  * [40] D. Kelly, A. J. Majda, and X. T. Tong, Nonlinear stability and ergodicity of ensemble based Kalman filters, Nonlinearity, 29 (2016), p. 657.
  * [41] N. B. Kovachki and A. M. Stuart, Ensemble Kalman inversion: a derivative-free technique for machine learning tasks, Inverse Problems, 35 (2019), p. 095005.
  * [42] E. Kwiatkowski and J. Mandel, Convergence of the square root ensemble Kalman filter in the large ensemble limit, SIAM/ASA Journal on Uncertainty Quantification, 3 (2015), pp. 1–17.
  * [43] T. Lange, Derivation of ensemble Kalman–Bucy filters with unbounded nonlinear coefficients, Nonlinearity, 35 (2021), p. 1061.
  * [44] T. Lange and W. Stannat, Mean field limit of ensemble square root filters - discrete and continuous time, Foundations of Data Science, 3 (2021), pp. 563–588.
  * [45] , On the continuous time limit of ensemble square root filters, Communications in Mathematical Sciences, 19 (2021), pp. 1855–1880.
  * [46] , On the continuous time limit of the ensemble Kalman filter, Math. Comput., 90 (2021), pp. 233–265.
  * [47] K. Law, H. Tembine, and R. Tempone, Deterministic mean-field ensemble Kalman filtering, SIAM Journal on Scientific Computing, 38 (2016), pp. A1251–A1279.
  * [48] F. Le Gland, V. Monbet, and V.-D. Tran, Large sample asymptotics for the ensemble Kalman filter, Research Report RR-7014, INRIA, 2009\.
  * [49] A. J. Majda and X. T. Tong, Performance of ensemble Kalman filters in large dimensions, Communications on Pure and Applied Mathematics, 71 (2018), pp. 892–937.
  * [50] F. Parzer and O. Scherzer, On convergence rates of adaptive ensemble Kalman inversion for linear ill-posed problems, ArXiv e-prints, Preprint arXiv:2104.10895 (2021).
  * [51] S. Reich, A dynamical systems framework for intermittent data assimilation, BIT Numerical Mathematics, 51 (2011), pp. 235–249.
  * [52] S. Reich and C. J. Cotter, Probabilistic Forecasting and Bayesian Data Assimilation, 2014.
  * [53] S. Reich and S. Weissmann, Fokker–Planck particle systems for Bayesian inference: Computational approaches, SIAM/ASA Journal on Uncertainty Quantification, 9 (2021), pp. 446–482.
  * [54] L. I. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena, 60 (1992), pp. 259 – 268.
  * [55] O. Scherzer, Convergence criteria of iterative methods based on Landweber iteration for solving nonlinear problems, Journal of Mathematical Analysis and Applications, 194 (1995), pp. 911–933.
  * [56] C. Schillings and A. M. Stuart, Analysis of the ensemble Kalman filter for inverse problems, SIAM Journal on Numerical Analysis, 55 (2017), pp. 1264–1290.
  * [57] , Convergence analysis of ensemble Kalman inversion: the linear, noisy case, Applicable Analysis, 97 (2018), pp. 107–123.
  * [58] A. M. Stuart, Inverse problems: A Bayesian perspective, Acta Numerica, 19 (2010), p. 451–559.
  * [59] X. Tong, A. Majda, and D. Kelly, Nonlinear stability of the ensemble Kalman filter with adaptive covariance inflation, Communications in Mathematical Sciences, 14 (2016), pp. 1283–1313.
  * [60] X. T. Tong, Performance analysis of local ensemble Kalman filter, Journal of Nonlinear Science, 28 (2018), pp. 1397–1442.
  * [61] X. T. Tong and M. Morzfeld, Localization in ensemble Kalman inversion, ArXiv e-prints, Preprint arXiv:2201.10821 (2022).
  * [62] S. Weissmann, N. K. Chada, C. Schillings, and X. T. Tong, Adaptive Tikhonov strategies for stochastic ensemble Kalman inversion, Inverse Problems, 38 (2022), p. 045009.