DEKI for high dimensional inverse problemsShuigen Liu, Sebastian Reich and Xin
T. Tong

# Dropout Ensemble Kalman inversion for high dimensional inverse
problems††thanks: Submitted to the editors DATE.

Shuigen Liu National University of Singapore<EMAIL_ADDRESS><EMAIL_ADDRESS>Sebastian Reich Institut für Mathematik, Universität
Potsdam<EMAIL_ADDRESS>Xin T. Tong22footnotemark: 2

###### Abstract

Ensemble Kalman inversion (EKI) is an ensemble-based method to solve inverse
problems. Its gradient-free formulation makes it an attractive tool for
problems with involved formulation. However, EKI suffers from the “subspace
property”, i.e., the EKI solutions are confined in the subspace spanned by the
initial ensemble. It implies that the ensemble size should be larger than the
problem dimension to ensure EKI’s convergence to the correct solution. Such
scaling of ensemble size is impractical and prevents the use of EKI in high
dimensional problems. To address this issue, we propose a novel approach using
dropout regularization to mitigate the subspace problem. We prove that
dropout-EKI converges in the small ensemble settings, and the computational
cost of the algorithm scales linearly with dimension. We also show that
dropout-EKI reaches the optimal query complexity, up to a constant factor.
Numerical examples demonstrate the effectiveness of our approach.

###### keywords:

Ensemble Kalman inversion, dropout, high dimension, small ensemble

##  1 Introduction

Inverse problems arise in various scientific and engineering problems [29, 17,
2, 9], where one is interested in recovering unknown parameters from indirect
and often noisy observations. Mathematically, the relationship between the
unknown parameter, denoted by $u\in\mathbb{R}^{d_{u}}$, and the observed data,
denoted by $y\in\mathbb{R}^{d_{y}}$, can often be described by the following
model:

(1) $y=\mathcal{G}(u)+\eta.$

Here $\mathcal{G}:\mathbb{R}^{d_{u}}\rightarrow\mathbb{R}^{d_{y}}$ is the
forward map. It describes the physical law behind the observation procedure,
and it is often corrupted by some observation noise, which is modeled by
$\eta\in\mathbb{R}^{d_{y}}$ in Eq. 1.

Directly solving Eq. 1 is often problematic, since the inverse of
$\mathcal{G}$ might not be well-defined, and we need to mitigate the
corruption from noise. A standard way to handle these issues is regularization
[2, 9]. In this paper, we will focus on the well-known Tikhonov regularization
[7] approach, which is also referred to as $l_{2}$ or ridge regularization.
Mathematically, this involves finding the minimizer $u^{*}$ of the loss
function

(2)
$l(u)=\frac{1}{2}\|y-\mathcal{G}(u)\|^{2}+\frac{1}{2}\|\mathcal{C}_{0}^{-1/2}u\|^{2}.$

The first part of Eq. 2 is the data misfit, which measures the fitness of a
candidate solution $u$ when compared with data. The second part of Eq. 2 is
the regularization term, it measures how well $u$ fits some prior information,
e.g. sparsity, intensity and smoothness. Here we choose the Tikhonov
regularization, so the prior information is encoded through a positive-
definite structural matrix $\mathcal{C}_{0}$. It is worth pointing out that
Eq. 2 can be seen as negative log-density of the Bayesian posterior of Eq. 1,
if we assume $u$ has prior $\mathcal{N}(0,\mathcal{C}_{0})$ and
$\eta\sim\mathcal{N}(0,I)$. In this context, the minimizer $u^{*}$ will become
the maximum-a-priori estimator of the posterior density. More details of such
connection can be found in [29, 15, 26].

Minimizing Eq. 2 can be computational challenging in practice. In many
applications, $\mathcal{G}$ is formulated by a PDE, of which the derivative
(or adjoint model) can be too complicate to derive or too large to process.
Ensemble Kalman inversion (EKI) is a derivative-free computational method
well-suited for these applications. In short, it employs an ensemble with a
Kalman-like formula to avoid the computation of derivatives. By iterating the
formula, the EKI ensemble will collapse to a single point which minimizes Eq.
2. EKI has shown impressive inversion skills in various problems [15, 18]. A
series of numerical analyses of EKI have been made recently to justify EKI’s
performance in linear, nonlinear and stochastic settings [26, 27, 3, 8, 33,
24].

One fundamental problem of EKI is its implementation on high dimensional
problems. This is characterized by the “subspace property” of EKI [15, 27],
which indicates that all ensemble members remains in the linear subspace
$\mathcal{V}_{0}$ spanned by the initial ensemble. This implies that if we
want to obtain $u^{*}$, $\mathcal{V}_{0}$ needs to be rich enough and the
ensemble size needs to surpass $d_{u}$ in general. But such scaling in
ensemble size is impractical for large-scale problems. In practice, this
problem can be mitigated by exploiting certain specific inverse problem
structures. For example, [6] discusses how to select the best
$\mathcal{V}_{0}$ so $u^{*}$ is likely to be in $\mathcal{V}_{0}$. In another
recent work, [30] considers localization techniques, which can leverage the
spatial dependence within an inverse problems. On the other hand, these
methods all require certain prior knowledge of the solution structures, so
they cannot be universal.

In this paper, we propose to tackle the subspace issue using dropout technique
originated from the deep learning literature [28, 32, 16, 10]. In the training
process of deep neural networks, the dropout technique randomly dropout
neurons to avoid overfitting. In theory, it is usually understood as a
regularization method [28, 32, 31] or a Bayesian neural network approximation
[10]. We will apply this idea and dropout random ensemble deviation components
in the EKI iterations. This leads to the dropout EKI (DEKI) scheme, which can
be implemented on high dimensional problems with a fixed ensemble size. We
prove that by properly adapting the step sizes, DEKI scheme converges
exponentially fast for linear problems, or strongly convex nonlinear problems
in the final stage. The overall computational cost of DEKI scales linearly
with dimension $d_{u}$, which is optimal for high dimensional problems. We
also investigate the query complexity of the zeroth order algorithms for
inverse problems, and show that the DEKI scheme achieves optimal query
complexity up to a constant in the problem class we consider.

The paper is organized as follows. In Section 2, we derive the DEKI scheme and
compare it with related methods. In Section 3, we provide the convergence
analysis of the DEKI, showing the exponential convergence rate. The query
complexity is also analyzed. In Section 4, we show the numerical performance
of DEKI on two examples, confirming our theoretic results.

###  1.1 Notations

For two matrices $A,B\in\mathbb{R}^{m\times n}$, denote the Hadamard product
$C=A\circ B\in\mathbb{R}^{m\times n}$ where $C_{ij}=A_{ij}B_{ij}$. For two
positive semi-definite (PSD) matrices $A,B\in\mathbb{R}^{m\times m}$, denote
$A\succeq B$ iff $A-B$ is PSD. PSD matrices admit standard eigenvalue
decomposition $A=Q\Lambda Q^{\rm T}$, where $\Lambda$ is a diagonal matrix and
$Q$ is an orthogonal matrix. Define the function of PSD matrix $f(A)$ by

$f(A)=Q~{}\text{diag}\\{f(\Lambda_{11}),\dots,f(\Lambda_{mm})\\}~{}Q^{\rm T}.$

For sequential random variables $X_{1},X_{2},\dots$, denote $\mathbb{E}_{n}$
as the conditional expectation conditioned on the information up until $n$,
i.e. $\mathbb{E}_{n}[Y]=\mathbb{E}[Y|\mathcal{F}_{n}]$ where $\mathcal{F}_{n}$
is the $\sigma$-algebra generated by $X_{1},\dots,X_{n}$.

##  2 Dropout EKI scheme

In this section, we start with a brief review of EKI formulation designed for
Eq. 1. We then implement the dropout technique and derive DEKI. Finally we
provide some comparison of DEKI with related EKI schemes.

###  2.1 Review: EKI and subspace property

When minimizing the regularized loss function Eq. 2, it is instructive to
formulate it as

(3) $l(u)=\frac{1}{2}\left\|z-\mathcal{H}(u)\right\|^{2},\quad\text{where
}\mathcal{H}(u)=\begin{bmatrix}\mathcal{G}(u)\\\
\mathcal{C}_{0}^{-1/2}u\end{bmatrix},~{}z=\begin{bmatrix}y\\\ 0\end{bmatrix}.$

[7] investigated applying vanilla EKI from [15] on Eq. 3. In particular, if we
denote the ensemble as $\\{u_{n}^{(j)}\\}_{j=1}^{J}$, the (Tikhonov) EKI
scheme reads

(4)
$u_{n+1}^{(j)}=u_{n}^{(j)}+h_{n}C_{n}^{up}(h_{n}C_{n}^{pp}+I)^{-1}(z-\mathcal{H}(u_{n}^{(j)})),$

where $C_{n}^{up},C_{n}^{pp}$ are the empirical covariance matrices of the
deviations. Denote

(5)
$\begin{split}&\overline{u}_{n}=\frac{1}{J}\sum_{j=1}^{J}u_{n}^{(j)},\quad\tau_{n}^{(j)}=u^{(j)}_{n}-\overline{u}_{n},\quad\overline{\mathcal{H}(u_{n})}=\frac{1}{J}\sum_{j=1}^{J}\mathcal{H}(u_{n}^{(j)}).\\\
&C_{n}^{uu}=\frac{1}{J-1}\sum_{j=1}^{J}\tau_{n}^{(j)}\otimes\tau_{n}^{(j)},\quad
C_{n}^{up}=\frac{1}{J-1}\sum_{j=1}^{J}\tau_{n}^{(j)}\otimes\left(\mathcal{H}(u_{n}^{(j)})-\overline{\mathcal{H}(u_{n})}\right)\\\
&C_{n}^{pp}=\frac{1}{J-1}\sum_{j=1}^{J}\left(\mathcal{H}(u_{n}^{(j)})-\overline{\mathcal{H}(u_{n})}\right)\otimes\left(\mathcal{H}(u_{n}^{(j)})-\overline{\mathcal{H}(u_{n})}\right).\end{split}$

The $h_{n}$ in Eq. 4 is the “stepsize” of the algorithm. The continuous limit
formulation can be derived if we fix $h_{n}\equiv h$ and let $h\to 0$. [8] has
shown that if $h_{n}$ increases with $n$, the convergence speed can be
improved.

The implementation of EKI involves iterating Eq. 4. In general, the EKI
ensemble will collapse onto a single point, and we can use the ensemble mean
$\bar{u}_{n}$ as a candidate minimizer of Eq. 1. When we have sufficient
ensemble, $J>d_{u}$, the update Eq. 4 can be interpreted as a Gauss-Newton
type iteration [8].

###  2.2 Naive DEKI

Using mathematical induction, it is not difficult to show that
$u_{n}^{(j)}\in\mathcal{V}_{0}=\text{span}\\{u_{0}^{(j)}\\}_{j=1}^{J}$, which
is the “subspace property” discussed in the introduction. At a high level, it
is not hard to break the “subspace property”, since we can modify Eq. 4 so the
new ensembles leave $\mathcal{V}_{0}$. But we also need some guidelines on how
to find the correct modification. We derive these guidelines from a
variational perspective. In specific, the fixed point of the iterations should
be be the minimizer (or a stationary point) of Eq. 2.

For the simplicity of illustration, we will assume the problem is linear, so
$\mathcal{H}(u)=Hu$ for some matrix $H$. In this case,
$C^{up}_{n}=C^{uu}_{n}H^{\rm T}$. Then if $u^{(j)}_{n}$ is the stationary
point of Eq. 4, we find

(6) $0=C^{uu}_{n}H^{\rm T}(h_{n}HC^{uu}_{n}H^{\rm T}+I)^{-1}(z-Hu_{n}^{(j)}).$

Meanwhile, $u^{(j)}_{n}$ is a minimizer of Eq. 2, iff $0=H^{\rm
T}(z-Hu_{n}^{(j)}).$ This can be guaranteed by Eq. 5 if $C^{uu}_{n}$ has full
rank $d_{u}$. But with vanilla EKI, the rank of $C^{uu}_{n}$ is at most $J-1$,
which is less than $d_{u}$ for $J\leq d_{u}$.

Next, we investigate how dropout technique can resolve this issue. Since rank
deficiency of $C^{uu}_{n}$ is the root of the problem, we consider applying
dropout on the ensemble deviations
$\tau_{n}^{(j)}=u_{n}^{(j)}-\overline{u}_{n}$, which leads to a dropout
ensemble $\\{\tilde{u}_{n}^{(j)}\\}_{j=1}^{J}$,

(7)
$\tilde{u}_{n}^{(j)}=\overline{u}_{n}+\tilde{\tau}_{n}^{(j)},\quad\tilde{\tau}_{n}^{(j)}=\rho_{n}\circ\tau_{n}^{(j)},\quad\rho_{n}(s)\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\text{Bernoulli}(\lambda).$

Here $1-\lambda\in(0,1)$ is the dropout rate, and $\circ$ is the Hadamard
product $(a\circ b)_{i}=a_{i}b_{i}$. Note that we will add tilde on symbols
related to the dropout ensemble. In particular, we have the dropout ensemble
covariance:

(8)
$\begin{split}&\tilde{C}_{n}^{uu}=\frac{1}{J-1}\sum_{j=1}^{J}\tilde{\tau}_{n}^{(j)}\otimes\tilde{\tau}_{n}^{(j)},\quad\tilde{C}_{n}^{up}=\frac{1}{J-1}\sum_{j=1}^{J}\tilde{\tau}_{n}^{(j)}\otimes\left(\mathcal{H}(\tilde{u}_{n}^{(j)})-\overline{\mathcal{H}(\tilde{u}_{n})}\right),\\\
&\tilde{C}_{n}^{pp}=\frac{1}{J-1}\sum_{j=1}^{J}\left(\mathcal{H}(\tilde{u}_{n}^{(j)})-\overline{\mathcal{H}(\tilde{u}_{n})}\right)\otimes\left(\mathcal{H}(\tilde{u}_{n}^{(j)})-\overline{\mathcal{H}(\tilde{u}_{n})}\right).\end{split}$

Replacing the empirical covariance $C_{n}^{up}$ and $C_{n}^{pp}$ in Eq. 4 with
their dropout versions will lead to the naive DEKI. From numerical
experiments, we find it can already resolve the “subspace property”, namely it
will converge to accurate numerical inversion result even if $J<d_{u}$,
although the convergence speed might be slow. We will add further numerical
modifications to accelerate the convergence.

But before that, it is instructive to provide the intuition why dropout can
help with the rank deficiency of sample covariance. Similar to other
stochastic algorithms such as stochastic gradient descent, the average effect
of using $\tilde{C}^{uu}_{n}$ can be characterized by the its conditional
average $\mathbb{E}_{n}\tilde{C}_{n}^{uu}$, which is taken conditioned on the
information of the current ensemble $\\{u_{n}^{(j)}\\}_{j=1}^{J}$. With some
elementary probability, it can be derived that

$\mathbb{E}_{n}\tilde{C}_{n}^{uu}=\lambda(1-\lambda)\text{diag}(C_{n}^{uu})+\lambda^{2}C_{n}^{uu}\succeq\lambda(1-\lambda)\text{diag}(C_{n}^{uu}),$

where $\text{diag}(C_{n}^{uu})$ is the diagonal part of $C_{n}^{uu}$, i.e.
$\text{diag}(C_{n}^{uu})(s,t)=\delta_{st}C_{n}^{uu}(s,t)$, where $\delta_{st}$
is the standard Kronecker delta. It is easy to see that
$\text{diag}(C_{n}^{uu})$ is full rank if the diagonal elements of
$C_{n}^{uu}$ are strictly positive.

###  2.3 DEKI with mean–deviation separation

The convergence of naive DEKI is slow because there is a tradeoff between
approximation and optimization: in order to approximate the nonlinear map
$\mathcal{H}$ more accurately, we would need the ensemble covariance collapse
quickly; but if the ensemble collapse too quickly, the ensemble will stop
moving based on Eq. 4. Moreover, analyzing the convergence of naive DEKI is
also very challenging, since $\tilde{C}^{uu}_{n}$ is very difficult to track.

An easy way to resolve this issue is to use separate dynamics for the mean and
deviation update. To illustrate, note that Eq. 4 can also be formulated as

(9a) $\displaystyle\overline{u}_{n+1}=~{}$
$\displaystyle\overline{u}_{n}+h_{n}C_{n}^{up}(I+h_{n}C_{n}^{pp})^{-1}(z-\overline{\mathcal{H}(u_{n})}),$
(9b) $\displaystyle\tau_{n+1}^{(j)}=~{}$
$\displaystyle\tau_{n}^{(j)}+h_{n}C_{n}^{up}(I+h_{n}C_{n}^{pp})^{-1}(\overline{\mathcal{H}(u_{n})}-\mathcal{H}(u_{n}^{(j)})).$

The separation of the two dynamics means we use different $h_{n},C_{n}^{up}$
and $C_{n}^{pp}$ in the DEKI.

Since the ensemble mean will be used as the candidate solution, we hope the
stationary point of its dynamics will also be a stationary point of $l(u)$.
Therefore we will use the dropout covariances in Eq. 9a. We will also use a
step size $\tilde{h}_{n}=\tilde{h}\|C^{uu}_{n}\|^{-1}$, an adaptive choice
adopted from [18]. Note the factor $\|C^{uu}_{n}\|^{-1}$ is added to counter
balance the fact that the ensemble will collapse, i.e. $C^{uu}_{n}\to 0$. This
leads to

(10)
$\overline{u}_{n+1}=\overline{u}_{n}+\tilde{h}_{n}\tilde{C}_{n}^{up}(I+\tilde{h}_{n}\tilde{C}_{n}^{pp})^{-1}(z-\mathcal{H}(\overline{u}_{n})),$

where $\tilde{C}_{n}^{up},\tilde{C}_{n}^{pp}$ are generated by Eq. 8.

As for the dynamics of ensemble deviation Eq. 9b, we seek a modification where
analysis of the ensemble covariance $C^{uu}_{n}$ is available. When
$\mathcal{H}$ is a linear map, $C^{uu}_{n}$ is known to have explicit Kalman
type update formula [15, 26]. For this reason, we seek a linearization of Eq.
9b, in the form of

(11) $\tau_{n+1}^{(j)}=\tau_{n}^{(j)}-h_{n}C_{n}^{uu}H_{n}^{\rm
T}(I+h_{n}H_{n}C_{n}^{uu}H_{n}^{\rm T})^{-1}H_{n}\tau_{n}^{(j)}.$

We set the step size $h_{n}=h\|C^{uu}_{n}\|^{-1}$, where $h$ is a constant
that can be smaller than $\tilde{h}$, so the ensemble will not collapse too
fast.

In order for Eq. 11 to approximate Eq. 9b, a natural choice will be
$H_{n}=C_{n}^{pu}(C_{n}^{uu})^{\dagger}$ where $(C_{n}^{uu})^{\dagger}$ is the
pseudo-inverse of $C_{n}^{uu}$. Notably, Eq. 11 will then be identical to Eq.
9b if $\mathcal{H}$ is linear. In our numerical test, this choice already
works. But for theoretical analysis below, we would require $H_{n}$ to be
bounded in $l_{2}$ norm, which may not hold for all nonlinear maps. On the
other hand, we note that a truncated version is always available as:

(12) $H_{n}:=\mathop{\arg\min}_{\|H\|_{2}\leq
M}\|C_{n}^{pu}(C_{n}^{uu})^{\dagger}-H\|_{\rm F}^{2}.$

$M$ here is a big prefixed constant. The optimization problem of Eq. 12 can be
solved simply by truncation on the singular value, see details in Appendix C.

In summary, the dropout EKI (DEKI) scheme is given by the mean update Eq. 10
using dropout ensemble covariance, and the deviation update Eq. 11 using a
truncated linearization Eq. 12. The algorithm is provided as Algorithm 1. We
note that while its formulation is more involved than EKI Eq. 4, DEKI’s
computational cost per iteration is only a constant multiple of the EKI. Our
numerical analysis below will show that DEKI will obtain a solution of
$\epsilon$ accuracy with $O(d_{u}\log\epsilon^{-1}/J)$ iterations when the
problem is strongly convex.

Algorithm 1 Dropout EKI

1: Input dropout rate $1-\lambda$, reference step size $h,\tilde{h}$.
Initialize $\\{u_{0}^{(j)}\\}_{j=1}^{J},n=0$.

2: while not converge do

3: Update $n\leftarrow n+1$.

4: Generate
$\rho_{n}(s)\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\text{Bernoulli}(\lambda)$
and compute dropout ensemble
$\tilde{u}_{n}^{(j)}=\overline{u}_{n}+\tilde{\tau}_{n}^{(j)},\quad\tilde{\tau}_{n}^{(j)}=\rho_{n}\circ\tau_{n}^{(j)}.$

5: Compute the dropout covariance $\tilde{C}_{n}^{up},\tilde{C}_{n}^{pp}$ by
Eq. 8.

6: Find the linearized map $H_{n}$ by solving Eq. 12.

7: Determine the step sizes
$h_{n}=h\|C^{uu}_{n}\|^{-1},\tilde{h}_{n}=\tilde{h}\|C^{uu}_{n}\|^{-1}$.

8: Update the mean and deviations separately by Eq. 10 and Eq. 11.

9: end while

###  2.4 Related literature

####  2.4.1 EKI methodology

EKI has been introduced in [15] as an extension of the popular ensemble Kalman
filter (EnKF) methodology to the treatment of inverse problems in the context
of optimization. Both the EnKF and EKI have recently been surveyed in [4]
primarily from a mean field perspective. However, practical algorithm use
finite ensemble sizes and establishing convergence and computational
efficiency for nonlinear problems have remained an open problem.

Here we focus first on the finite ensemble size effect and later also address
the issue of computational efficiency. As discussed before, the subspace
property will prevent the ensemble to converge to the optimal solution.
However, some asymptotic behavior of EKI can be analyzed in linear [26, 27]
and nonlinear [7, 33] cases. Roughly speaking, the ensemble converges to the
optimal solution projected onto the subspace spanned by the initial ensemble.
Localization [30] is one method to overcome this issue, which will be
discussed later in detail. Interested readers are referred to [8, 33] for
further recent results on EKI.

EKI can also be understood as a zeroth order optimization method [15, 7, 8].
In particular, it can be viewed as a gradient flow using the empirical
covariance as a preconditioner [26, 33]. It approximates the gradients by the
finite differences, making the method derivative-free. Among zeroth order
optimization methods, it is known that there are methods with complexity
$\mathcal{O}(d_{u}\log\varepsilon^{-1})$ to achieve $\varepsilon$ accuracy for
strongly convex functions, such as the Nesterov random search method [22]. The
complexity of DEKI matches this result, and in fact we show that this rate is
optimal even for linear inverse problems.

####  2.4.2 Dropout in deep learning

Dropout technique is first introduced in neural network (NN) by Hinton et al.
[28] as an implicit regularization technique to avoid overfitting. The idea is
to randomly omit parts of the neurons in a NN during training. This method has
proven to be very effective in improving the performance of NN in a wide range
of tasks [28] and is accepted as a common practice now. There is a long
literature exploring the practical implementation of different dropout methods
in different types of NNs. Popular methods include the original dropout [28],
dropconnect [32], variational dropout [16], Monte Carlo dropout [10] etc.

Despite its rich application in deep learning, theoretic understanding of
dropout is still limited. Most practitioners regard dropout as an implicit
regularization method [28, 32, 31]. By randomly dropping out neurons, dropout
decreases the dependence on individual neurons to avoid overfitting. As a
regularization method, dropout can also be viewed as a noisy perturbation of
the output of the neurons, which encourages the model to learn more robust
features. Another explanation is to interpret dropout as a Bayesian
approximation [10], where a probability distribution is imposed over the
parameters, and dropout is used to approximate a stochastic realization of the
parameters.

The dropout method we use in DEKI is similar to dropout NN only at a high
level. Dropout in DEKI decreases the dependence on different parameters to
avoid spurious long-distance correlation, which is a well-known issue
especially in the high dimensional and small ensemble settings. On the other
hand, the exact reason and benefits of dropout for EKI should be analyzed
independently.

####  2.4.3 Localized EKI

Localization is a crucial technique used in EnKF to deal with the sampling
errors due to a limited ensemble size [13, 14, 11, 23]. The idea is to
artificially reduces the the covariance between distant locations, typically
by applying a distance-based weighting function. [30] investigated the
application of localization in EKI, and showed that the localized EKI (LEKI)
scheme can be applied in the small ensemble settings and overcome the subspace
issue of EKI effectively.

Our scheme bears a resemblance to the LEKI, but avoids the potentially-
involved design of the localization methods, and is computationally cheaper
since DEKI can be implemented using only matrix-vector multiplication while
LEKI cannot. At a high level, DEKI can also be viewed as a stochastic
implementation of the LEKI, see Eq. 23 for details. We also emphasize that
there are other possible ways to apply dropout in EKI, and here we just choose
one version to illustrate the general idea and show theoretic properties.

##  3 Analysis of DEKI

The analysis is divided into two parts. First we show the ensemble collapses
in a controllable way by studying the dynamic of the spectrum of the ensemble
covariance, from which we obtain lower bounds of the diagonal entries. Second
we prove the convergence of DEKI by approximating it locally by a Gauss–Newton
type scheme and derive contraction using the lower bounds obtained before.
Here an important observation is that the dropout emphasizes the diagonal
elements of the ensemble covariance in average, see Eq. 23 for details.

###  3.1 Assumptions

We make the following assumptions on the regularized forward operator
$\mathcal{H}$ Eq. 3, the linearized map $H_{n}$ Eq. 12 and the $l^{2}$-loss
$l(u)$ Eq. 2. The following conditions hold:

  1. 1.

Bounded Hessian. $\mathcal{H}\in C^{2}$ and $\exists H>0$ s.t.

(13)
$\sup_{u}\left(\sum_{s=1}^{d_{z}}\left\|\nabla^{2}\mathcal{H}_{s}(u)\right\|_{2}^{2}\right)^{1/2}\leq
H,$

where $\mathcal{H}_{s}$ denotes the $s$-th component of $\mathcal{H}$.

  2. 2.

Boundedness of $H_{n}$. $\exists 0<\gamma\leq M$ s.t. $\forall
n\in\mathbb{N}$,

(14) $\gamma^{2}I\preceq H_{n}^{\rm T}H_{n}\preceq M^{2}I.$

  3. 3.

$L$-smoothness of $l(u)$. $\exists L>0$, s.t. $\forall
u,v\in\mathbb{R}^{d_{u}}$,

(15) $\|\nabla l(u)-\nabla l(v)\|\leq L\left\|u-v\right\|.$

  4. 4.

Polyak-Łojasiewicz (PL) condition for $l(u)$. $\exists c>0$ s.t. $\forall
u\in\mathbb{R}^{d_{u}}$,

(16) $\left\|\nabla l(u)\right\|^{2}\geq c(l(u)-l_{\min}).$

(1) The assumptions on $\mathcal{H}$ and $H_{n}$ can be obtained from the
according properties of the original map $\mathcal{G}$. For instance, we can
take $H_{n}$ of the form $H_{n}=[G_{n},\mathcal{C}_{0}^{-1/2}]^{\rm T}$, where
$G_{n}$ is a linearization of $\mathcal{G}$ and $\left\|G_{n}\right\|\leq
M_{0}$. If the regularization operator satisfies
$\gamma^{2}I\preceq\mathcal{C}_{0}^{-1}\preceq\Gamma^{2}I$, then Eq. 14 holds
with $\gamma$ and $M=\sqrt{M_{0}^{2}+\Gamma^{2}}$.

(2) The $L$-smoothness and PL condition for $l(u)$ are standard assumptions in
optimization. Note PL condition is a weaker condition than strong convexity,
and is widely adopted in the machine learning literature when lacking
convexity [19].

###  3.2 Controllable ensemble collapse

The linearization Eq. 12 applied to the covariance evolution produces a
behavior similar to that of a linear forward map. It allows us to have good
controls on the ensemble collapse as follows:

######  Proposition 3.1.

Consider the DEKI scheme for the regularized problem Eq. 3. Assume that the
linearized map $H_{n}$ satisfies $\gamma^{2}I\preceq H_{n}^{\rm T}H_{n}\preceq
M^{2}I$. Denote $C_{n}^{uu}$ as the empirical covariance. Choose adaptive step
sizes $h_{n}=\theta\left\|C_{n}^{uu}\right\|^{-1}$, where $\theta\leq M^{-2}$.
Then the modified condition number of $C_{n}^{uu}$ is uniformly bounded:

(17) $\forall
n\in\mathbb{N},\quad\kappa(C_{n}^{uu}):=\frac{\lambda_{1}(C_{n}^{uu})}{\lambda_{r}(C_{n}^{uu})}\leq\bar{\kappa}=\max\left\\{\kappa(C_{0}^{uu}),\frac{3M^{2}}{2\gamma^{2}}\right\\}.$

where $\lambda_{k}(C_{n}^{uu})$ is the $k$-th largest eigenvalues of
$C_{n}^{uu}$ and $r$ is the rank of $C_{n}^{uu}$, which is invariant under
evolution: $\text{rank}(C_{n}^{uu})\equiv r:=\text{rank}(C_{0}^{uu})$.
Moreover, the ensemble collapses exponentially,

(18) $\|C_{n}^{uu}\|\leq\|C_{0}^{uu}\|(1+\gamma^{2}\theta)^{-2n}.$

######  Proof 3.2.

Rewrite the deviation update Eq. 11 as

$\begin{split}\tau_{n+1}^{(j)}=~{}&\left[I-h_{n}C_{n}^{uu}H_{n}^{\rm
T}(I+h_{n}H_{n}C_{n}^{uu}H_{n}^{\rm T})^{-1}H_{n}\right]\tau_{n}^{(j)}\\\
=~{}&(I+h_{n}C_{n}^{uu}H_{n}^{\rm T}H_{n})^{-1}\tau_{n}^{(j)},\end{split}$

where we use the Sherman-Morrison-Woodbury formula (see Lemma B.1) by taking
$A=I,U=h_{n}C_{n}^{uu}H_{n}^{\rm T}$ and $V=H_{n}^{\rm T}$.

Recall by definition
$C_{n}^{uu}=\frac{1}{J-1}\sum_{j=1}^{J}\tau_{n}^{(j)}(\tau_{n}^{(j)})^{\rm
T}$. Therefore,

$C_{n+1}^{uu}=(I+h_{n}C_{n}^{uu}H_{n}^{\rm
T}H_{n})^{-1}C_{n}^{uu}(I+h_{n}H_{n}^{\rm T}H_{n}C_{n}^{uu})^{-1}.$

By assumption, $\gamma^{2}I\preceq H_{n}^{\rm T}H_{n}\preceq M^{2}I.$ Apply
Lemma B.4 to $C_{n+1}^{uu}$, we obtain

(19) $C_{n}^{uu}(I+M^{2}h_{n}C_{n}^{uu})^{-2}\preceq C_{n+1}^{uu}\preceq
C_{n}^{uu}(I+\gamma^{2}h_{n}C_{n}^{uu})^{-2}.$

Notice when $0\leq x\leq a^{-1}$,

$\frac{x}{1+3ax}\leq\frac{x}{(1+ax)^{2}}\leq\frac{x}{1+2ax}.$

Applying Lemma B.2 (3), when $h_{n}=\theta\|C_{n}^{uu}\|^{-1}$ and $\theta\leq
M^{-2}\leq\gamma^{-2}$, it holds

$C_{n}^{uu}(I+3M^{2}h_{n}C_{n}^{uu})^{-1}\preceq C_{n+1}^{uu}\preceq
C_{n}^{uu}(I+2\gamma^{2}h_{n}C_{n}^{uu})^{-1}.$

Denote $\lambda_{k}(C_{n}^{uu})$ as the $k$-th largest eigenvalues of
$C_{n}^{uu}$, then (see Lemma B.2 (2))

(20)
$\frac{\lambda_{k}(C_{n}^{uu})}{1+3h_{n}M^{2}\lambda_{k}(C_{n}^{uu})}\leq\lambda_{k}(C_{n+1}^{uu})\leq\frac{\lambda_{k}(C_{n}^{uu})}{1+2h_{n}\gamma^{2}\lambda_{k}(C_{n}^{uu})}.$

Hence $\lambda_{k}(C_{n}^{uu})>0\Leftrightarrow\lambda_{k}(C_{n+1}^{uu})>0$,
which implies that the rank of $C_{n}^{uu}$ remains constant
$r=\text{rank}(C_{0}^{uu})$. From Eq. 20, we get for $1\leq k\leq r$,

$2h_{n}\gamma^{2}\leq\lambda_{k}^{-1}(C_{n+1}^{uu})-\lambda_{k}^{-1}(C_{n}^{uu})\leq
3h_{n}M^{2}.$

Denote the modified condition number
$\kappa_{n}=\lambda_{1}(C_{n}^{uu})/\lambda_{r}(C_{n}^{uu})$. Then

$\begin{split}\kappa_{n+1}=~{}&\frac{\lambda_{r}^{-1}(C_{n+1}^{uu})}{\lambda_{1}^{-1}(C_{n+1}^{uu})}\leq\frac{\lambda_{r}^{-1}(C_{n}^{uu})+3h_{n}M^{2}}{\lambda_{1}^{-1}(C_{n}^{uu})+2h_{n}\gamma^{2}}\\\
\leq~{}&\max\left\\{\frac{\lambda_{r}^{-1}(C_{n}^{uu})}{\lambda_{1}^{-1}(C_{n}^{uu})},\frac{3M^{2}}{2\gamma^{2}}\right\\}=\max\\{\kappa_{n},3M^{2}/2\gamma^{2}\\}.\end{split}$

By induction, we deduce that

$\kappa_{n}\leq\bar{\kappa}:=\max\\{\kappa_{0},3M^{2}/2\gamma^{2}\\}.$

Finally, applying Lemma B.2 (1) to Eq. 19, we obtain

$\begin{split}\|C_{n+1}^{uu}\|\leq\lambda_{\max}\Big{(}&C_{n}^{uu}(I+\gamma^{2}h_{n}C_{n}^{uu})^{-2}\Big{)}=\frac{\|C_{n}^{uu}\|}{(1+\gamma^{2}h_{n}\|C_{n}^{uu}\|)^{2}}.\\\
\Rightarrow~{}&\|C_{n}^{uu}\|\leq\|C_{0}^{uu}\|(1+\gamma^{2}\theta)^{-2n}.\end{split}$

Note the equality above holds since $\psi(x)=x(1+\gamma^{2}h_{n}x)^{-2}$ is
monotone increasing on $x\in[0,\|C_{n}^{uu}\|]$, due to the fact that when
$h_{n}=\theta\|C_{n}^{uu}\|^{-1}$ and $\theta\leq M^{-2}\leq\gamma^{-2}$,

$\gamma^{2}h_{n}x\leq\gamma^{2}h_{n}\|C_{n}^{uu}\|\leq
1\Rightarrow\psi^{\prime}(x)=\frac{1-\gamma^{2}h_{n}x}{(1+\gamma^{2}h_{n}x)^{2}}\geq
0.$

The convergence analysis also requires the control of the diagonal elements of
the covariance matrix, whose proof is mainly based on the invariance of the
column space of the covariance matrix. We prove

######  Proposition 3.3.

Under the settings of Proposition 3.1. The lower bound for the diagonal
elements of $C_{n}^{uu}$ holds:

(21)
$\begin{split}\min_{s}C_{n}^{uu}&(s,s)\geq\bar{\kappa}^{-1}\left\|C_{n}^{uu}\right\|\min_{s}P(s,s),\end{split}$

where $\bar{\kappa}$ is the upper bound in Eq. 17, and
$P\in\mathbb{R}^{d_{u}\times d_{u}}$ is the the orthogonal projector onto the
column space $\mathcal{V}_{n}=\text{Im}(C_{n}^{uu})$. The projector is
invariant under evolution since $\mathcal{V}_{n}$ is invariant:
$\mathcal{V}_{n}\equiv\mathcal{V}=\text{Im}(C_{0}^{uu})$.

######  Remark 1.

To obtain a meaningful lower bound, we should assume

(22) $\min_{s}P(s,s)>0.$

It is equivalent to $\min_{s}C_{0}^{uu}(s,s)>0$, which is necessary since the
dropout procedure cannot explore directions with zero covariance. This
condition can be achieved easily by using Gaussian initialization.

######  Proof 3.4.

First we prove that the column space of $C_{n}^{uu}$ is invariant. Note
$\text{Im}(C_{n}^{uu})=\text{span}\\{\tau_{n}^{(j)}:j=1,\dots J\\}$. Observe
that for any $v\in\mathbb{R}^{d_{u}}$ that satisfies $\forall 1\leq j\leq
J,v^{\rm T}\tau_{n}^{(j)}=0$, it holds

$v^{\rm T}\tau_{n+1}^{(k)}=v^{\rm
T}\left(\tau_{n}^{(k)}+\frac{h_{n}}{J-1}\sum_{j=1}^{J}\tau_{n}^{(j)}(H_{n}(u_{n}^{(j)})-\overline{H_{n}(u_{n})})^{\rm
T}w_{n}^{(j)}\right)=0,$

where we denote $w_{n}^{(j)}=(I+h_{n}H_{n}C_{n}^{uu}H_{n}^{\rm
T})^{-1}(\overline{H_{n}(u_{n})}-H_{n}(u_{n}^{(j)}))$. Thus by induction, we
obtain $\forall n,1\leq j\leq J,v^{\rm T}\tau_{n}^{(j)}=0$. This implies that

$\text{Im}(C_{n}^{uu})\subset\text{Im}(C_{0}^{uu}).$

While by Proposition 3.1, $\text{rank}(C_{n}^{uu})=\text{rank}(C_{0}^{uu})$,
we conclude that $\text{Im}(C_{n}^{uu})=\text{Im}(C_{0}^{uu})$.

To control the diagonal element of $C_{n}^{uu}$, consider the singular value
decomposition $C_{n}^{uu}=Q_{n}\Lambda_{n}Q_{n}^{\rm T}$, then

$C_{n}^{uu}(s,s)=\sum_{t=1}^{r}\lambda_{t}(C_{n}^{uu})(Q_{n}(s,t))^{2}.$

Notice $P=Q_{n}Q_{n}^{\rm T}$ is the orthogonal projector onto the invariant
colunm space $\mathcal{V}=\text{Im}(C_{0}^{uu})$, and therefore is invariant.
Then we can control the lower bound of the diagonal element as follows:

$\min_{s}C_{n}^{uu}(s,s)\geq\lambda_{r}(C_{n}^{uu})\min_{s}\sum_{t=1}^{r}(Q_{n}(s,t))^{2}\geq\bar{\kappa}^{-1}\|C_{n}^{uu}\|\min_{s}P(s,s).$

###  3.3 Convergence to optimal solution

The convergence of DEKI is mainly guaranteed by the positive definiteness of
average dropout ensemble covariance, this is similar to the localization
effect used in LEKI. In particular, we claim that

(23)
$\bar{C}_{n}^{uu}=\mathbb{E}_{n}\tilde{C}_{n}^{uu}=\lambda(1-\lambda)\text{diag}(C_{n}^{uu})+\lambda^{2}C_{n}^{uu},$

where $\text{diag}(C_{n}^{uu})$ is the diagonal part of $C_{n}^{uu}$, i.e.
$\text{diag}(C_{n}^{uu})(s,t)=\delta_{st}C_{n}^{uu}(s,t)$, and here
$\delta_{st}$ is the standard kronecker delta. If we denote
$\Psi=\lambda^{2}I+\lambda(1-\lambda)E$, where $E(s,t)\equiv 1$, then
$\bar{C}_{n}^{uu}=\Psi\circ C_{n}^{uu}$.

To prove Eq. 23, just notice when $s\neq t$,

$\begin{split}\bar{C}_{n}^{uu}(s,t)=~{}&\mathbb{E}_{n}\left(\frac{1}{J-1}\sum_{j=1}^{J}\rho(s)\rho(t)\tau_{n}^{(j)}(s)\tau_{n}^{(j)}(t)\right)\\\
=~{}&\frac{1}{J-1}\sum_{j=1}^{J}\lambda^{2}\tau_{n}^{(j)}(s)\tau_{n}^{(j)}(t)=\lambda^{2}C_{n}^{uu}(s,t).\end{split}$

And when $s=t$,

$\begin{split}\bar{C}_{n}^{uu}(s,s)=~{}&\mathbb{E}_{n}\left(\frac{1}{J-1}\sum_{j=1}^{J}\rho^{2}(s)(\tau_{n}^{(j)}(s))^{2}\right)\\\
=~{}&\frac{1}{J-1}\sum_{j=1}^{J}\lambda(\tau_{n}^{(j)}(s))^{2}=\lambda
C_{n}^{uu}(s,s).\end{split}$

The convergence analysis is based on the connection of EKI with Gauss-Newton.
To be specific, notice the mean evolution Eq. 10 in DEKI can be approximated
by the following Gauss-Newton type iteration

(24)
$\overline{u}_{n+1}^{\prime}=\overline{u}_{n}+\tilde{h}_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T}(I+\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T})^{-1}(z-\mathcal{H}(\overline{u}_{n})),$

where $G_{n}=\nabla\mathcal{H}(\overline{u}_{n})$ is the Jacobian at
$\overline{u}_{n}$. Note it is the Gauss-Newton update for the loss function

$l_{n}(u)=\|u-\overline{u}_{n}\|_{C_{n}^{uu}}^{2}+h_{n}\|\mathcal{H}(u)-z\|^{2}.$

We can bound the difference of the two updates by the following lemma. It
takes a similar form as Proposition 3.3 in [8], but the upper bound depends on
$d_{u}$ there. Here we improve the bounds so that there is no dimension
dependence.

######  Proposition 3.5.

Under the settings of Theorem 3.6, it holds

(25) $\|\overline{u}_{n+1}-\overline{u}_{n+1}^{\prime}\|\leq
C\tilde{h}_{n}\|\tilde{C}_{n}^{uu}\|^{3/2}\|z-\mathcal{H}(\overline{u}_{n})\|,$

where $C=(J-1)^{3/2}H$ is a constant independent of $d_{u}$ and $H$ is the
bound of the Hessian of $\mathcal{H}$ defined in Eq. 13.

The proof is postponed to Section A.1. The proposition shows that DEKI can be
approximated well by the Gauss-Newton update as the ensemble collapses. It is
expected to obtain exponential convergence for such scheme under the PL
condition Eq. 16. To be specific, we prove the main theorem as follows.

######  Theorem 3.6.

Consider the DEKI scheme for the regularized problem Eq. 3 in the high
dimensional regime $d_{u}\gg J$. Assume Section 3.1 holds, and initialize the
ensemble so that Eq. 22 holds. Choose adaptive step sizes

(26)
$h_{n}=\theta\|C_{n}^{uu}\|^{-1},\quad\tilde{h}_{n}=\mu\|C_{n}^{uu}\|^{-1},$

where $\theta\leq M^{-2},\mu\leq L^{-1}$, and $L,M$ are defined in Section
3.1. Then there exists $n_{0}\leq\mathcal{O}(\log d_{u})$ s.t. when $n>n_{0}$,
the $l^{2}$-loss Eq. 2 converges exponentially,

(27)
$\mathbb{E}_{n_{0}}l(\overline{u}_{n})-l_{\min}\leq(1-\beta)^{n-n_{0}}(l(\overline{u}_{n_{0}})-l_{\min}+C_{2}),$

where $\beta=\mathcal{O}(J/d_{u})$ and $C_{2}$ are some constants that can be
determined explicitly.

(1) Since $\beta=\mathcal{O}(J/d_{u})$, to reach $\varepsilon$-accuracy, the
needed steps are

(28) $N_{\varepsilon}=\mathcal{O}\left(d_{u}\log\varepsilon^{-1}/J\right).$

Such complexity scaling to the dimension $d_{u}$ (linear) is optimal. For more
details, please see Proposition 3.9.

(2) Without the high dimensional assumption $d_{u}\gg J$, one can still show
that

(29)
$\mathbb{E}_{n_{0}}l(\overline{u}_{n})-l_{\min}\leq(1-\beta)^{n-n_{0}}(l(\overline{u}_{n_{0}})-l_{\min})+C_{1}(n-n_{0})\max\\{1-\beta,\delta\\}^{n-1},$

where $\beta,\delta$ are chosen as

(30)
$\beta=c\lambda(1-\lambda)\bar{\kappa}^{-1}\min_{s}P(s,s)\cdot\frac{\mu(1+2\mu
M^{2})}{4(1+\mu M^{2})^{2}},\quad\delta=(1+\gamma^{2}\theta)^{-1}.$

Here $\bar{\kappa}$ is the upper bound in Eq. 17, and
$P\in\mathbb{R}^{d_{u}\times d_{u}}$ is the orthogonal projector onto
$\mathcal{V}=$ Im$(C_{0}^{uu})$ introduced in Proposition 3.3. Note $\beta$ is
the same as in the theorem.

(3) Here $n_{0}\leq\mathcal{O}(\log d_{u})$ is the warm-up steps. Note
$\mathcal{O}(\log d_{u})$ is usually viewed as $\mathcal{O}(1)$ quantity, and
the term $l(\overline{u}_{n_{0}})$ is harmless to our result. From Lemma A.1,
we can see that the $l^{2}$-loss can only have constant growth:
$l(\overline{u}_{n_{0}})\leq C^{\prime}l(\overline{u}_{0})$.

(4) Computational costs comparision. For vanilla EKI, the ensemble size should
be at least $J=\mathcal{O}(d_{u})$, and it takes
$\mathcal{O}(\log\varepsilon^{-1})$ steps to reach $\varepsilon$-accuracy. In
each step, it takes $J$ evaluations of the forward map, and
$\mathcal{O}(J^{2}(d_{u}+d_{y}))$ additional costs due to the matrix-vector
multiplications. Thus the overall costs would be
$\mathcal{O}(d_{u}\log\varepsilon^{-1})$ forward evaluations and
$\mathcal{O}(d_{u}^{2}(d_{u}+d_{y})\log\varepsilon^{-1}))$ additional costs.
For DEKI, it takes $\mathcal{O}(J^{-1}d_{u}\log\varepsilon^{-1})$ steps when
$\varepsilon$ is small, and thus the overall costs will be
$\mathcal{O}(d_{u}\log\varepsilon^{-1})$ forward evaluations, and
$\mathcal{O}(Jd_{u}(d_{u}+d_{y})\log\varepsilon^{-1})$ additional costs. Since
$J\ll d_{u}$, we can see that DEKI reduces the additional cost while using the
same order of forward evaluations.

######  Proof 3.7.

The proof is based on the one step estimation in Lemma 3.8 stated at the end
of this section:

$\mathbb{E}_{n}l(\overline{u}_{n+1})-l_{\min}\leq(1-2\beta+\Delta_{n})(l(\overline{u}_{n})-l_{\min})+\Delta_{n}l_{\min},$

By Proposition 3.1, the ensemble collapses exponentially,

$\|C_{n}^{uu}\|\leq\|C_{0}^{uu}\|\delta^{2n},\quad\delta=(1+\gamma^{2}\theta)^{-1}<1.$

Then we can bound $\Delta_{n}$ by

$\Delta_{n}\leq C_{1}\delta^{n},$

for
$C_{1}=C\mu\|C_{0}^{uu}\|^{1/2}\left[2+2L\mu^{1/2}+LC\mu\|C_{0}^{uu}\|^{1/2}\right]=\mathcal{O}(d_{u})$.

To obtain convergence, we pick $n_{0}$ so that $\forall n\geq
n_{0},\Delta_{n}\leq\beta$. One can take

$n_{0}=\Big{\lceil}\frac{\log(2\beta^{-1}C_{1})}{\log\delta^{-1}}\Big{\rceil}\leq\mathcal{O}(\log
d_{u}).$

Thus when $n>n_{0}$,

$\mathbb{E}_{n}l(\overline{u}_{n+1})-l_{\min}\leq(1-\beta)(l(\overline{u}_{n})-l_{\min})+C_{1}\delta^{n}l_{\min},$

By induction we obtain

(31) $\begin{split}&\mathbb{E}_{n_{0}}l(\overline{u}_{n})-l_{\min}\\\
\leq~{}&(1-\beta)^{n-n_{0}}(l(\overline{u}_{n_{0}})-l_{\min})+C_{1}l_{\min}\sum_{k=n_{0}}^{n-1}(1-\beta)^{n-1-k}\delta^{k}.\end{split}$

It is elementary to see that

$\sum_{k=n_{0}}^{n-1}(1-\beta)^{n-1-k}\delta^{k}\leq\sum_{k=n_{0}}^{n-1}\max\\{1-\beta,\delta\\}^{n-1}=(n-n_{0})\max\\{1-\beta,\delta\\}^{n-1}.$

From which we obtain (denote $C_{1}^{\prime}=l_{\min}C_{1}$)

$\mathbb{E}_{n_{0}}l(\overline{u}_{n})-l_{\min}\leq(1-\beta)^{n-n_{0}}(l(\overline{u}_{n_{0}})-l_{\min})+C_{1}^{\prime}(n-n_{0})\max\\{1-\beta,\delta\\}^{n-1}.$

Next we consider the high dimensional settings, i.e. $J\ll d_{u}$. First
notice

$c,L=\mathcal{O}(M^{2}),\quad\lambda,\bar{\kappa},\mu,\mu,C=\mathcal{O}(1).$

Assume that we initialize $\\{u_{0}^{(j)}\\}_{j=1}^{J}$ with i.i.d. Gaussian
variables. Then in the regime $J\ll d_{u}$, since
$r=\text{rank}(C_{0}^{uu})\approx J-1$,

$\text{tr}(C_{0}^{uu})=\mathcal{O}(d_{u})\Rightarrow\|C_{0}^{uu}\|=\mathcal{O}(d_{u}/J).$

Similarly, $P$ is an orthogonal projector onto the $r$-dimensional subspace,
so that

$\text{tr}(P)=\mathcal{O}(J)\Rightarrow\min_{s}P(s,s)=\mathcal{O}(J/d_{u}).$

From this, we can see that

$\beta=c\lambda(1-\lambda)\bar{\kappa}^{-1}\min_{s}P(s,s)\cdot\frac{\mu(1+2\mu
M^{2})}{4(1+\mu M^{2})^{2}}=\mathcal{O}(J/d_{u}).$

Since $\delta=\mathcal{O}(1)$, it holds that
$1-\beta-\delta=\mathcal{O}(1)>0$, so that we can bound the last term in Eq.
31 with a tighter bound:

$\sum_{k=n_{0}}^{n-1}(1-\beta)^{n-1-k}\delta^{k}=(1-\beta)^{n-1}\sum_{k=n_{0}}^{n-1}\left(\frac{\delta}{1-\beta}\right)^{k}\leq\frac{(1-\beta)^{n}}{1-\beta-\delta}.$

And the convergence result can be improved to be

$\mathbb{E}_{n_{0}}l(\overline{u}_{n})-l_{\min}\leq(1-\beta)^{n-n_{0}}(l(\overline{u}_{n_{0}})-l_{\min}+C_{2}),$

where $C_{2}=(1-\beta)^{n_{0}}C_{1}^{\prime}/(1-\beta-\delta)$.

######  Lemma 3.8.

Under the settings of Theorem 3.6, it holds

(32)
$\mathbb{E}_{n}l(\overline{u}_{n+1})-l_{\min}\leq(1-2\beta+\Delta_{n})(l(\overline{u}_{n})-l_{\min})+\Delta_{n}l_{\min},$

where $\beta,\Delta_{n}$ are determined by

$\beta=c\lambda(1-\lambda)\bar{\kappa}^{-1}\min_{s}P(s,s)\cdot\frac{\mu(1+2\mu
M^{2})}{4(1+\mu M^{2})^{2}}.$
$\Delta_{n}=2C\mu(M+L\mu^{1/2})\|C_{n}^{uu}\|^{1/2}+LC^{2}\mu^{2}\|C_{n}^{uu}\|.$

See Section A.2 for proof.

###  3.4 Query complexity analysis

In this section, we consider the lower bound of query complexity for zeroth
order algorithms to solve inverse problems, which is closely related to the
optimization problems. In the optimization literature, zeroth order algorithms
are those who can only evaluate the loss function $l$ at a query point, but
not the derivatives of $l$. The query complexity is the number of query points
needed for the algorithm to reach a certain accuracy. We mention that since
the first formal study in [20], a lot of results are obtained for different
types of optimization algorithms and problem classes, for instance [21, 5, 1].

It is natural to extend to the inverse problem setting for zeroth order
algorithms that can only evaluate $\mathcal{G}(u)$ at a query point $u$.
Ensemble Kalman type algorithms clearly are zeroth order algorithms, and their
query complexity will be $nJ$, where $n$ is iteration number needed and $J$ is
the ensemble size. Query complexity is important because the evaluation of
$\mathcal{G}$, which is often a simulation of an expensive black box model, is
the main cost of the operation. Theorem 3.6 indicates the query complexity of
DEKI is $\mathcal{O}(d_{u}\log\varepsilon^{-1})$. Next, we show such linear
dependence on $d_{u}$ is optimal.

Mathematically speaking, a zeroth order algorithm $A$ solving Eq. 1 can be
described as a map of the form

$A:(U,Y,y)\mapsto u.$

Here $U=(u_{1},\dots,u_{n})\in\mathbb{R}^{d_{u}\times n}$ are the query
points, $Y=(y_{1},\dots,y_{n})\in\mathbb{R}^{d_{y}\times n}$ are the model
outputs (i.e. $y_{i}=\mathcal{G}(u_{i})$), $y\in\mathbb{R}^{d_{y}}$ is the
observation data that we try to invert, and $u\in\mathbb{R}^{d_{u}}$ is the
output of the algorithm. We say $A(U,Y,y)$ is a zeroth order algorithm with
$n$ admissible query, if the query points $u_{1},\ldots,u_{n}$ are obtained
from some admissible rule, i.e. $u_{i+1}=B(y,u_{t},\mathcal{G}(u_{t}),t\leq
i)$ for some update rule $B$, and thus the output is a map
$A(U,\mathcal{G}(U),Y)$.

Next we show that even for linear inverse problems, with $n\leq d_{u}/2$ query
points, any algorithms cannot identify the optimal solution to Eq. 2
accurately.

######  Proposition 3.9.

Consider solving Eq. 2 for linear forward maps with regularization operator
$\mathcal{C}_{0}^{-1/2}=I$. Then for any zeroth order algorithm $A$ using any
$n$ admissible query $U\in\mathbb{R}^{d_{u}\times n}$ where $n\leq d_{u}/2$,
it holds that

$\max_{G:\|G\|\leq 2}\|A(U,GU,y)-u^{*}(G)\|>0.1\left\|y\right\|,$

where $u^{*}(G)$ is the optimal solution to Eq. 2.

Proofs can be found in Section A.3. The proposition implies that we need at
least $n>d_{u}/2$ queries of the forward map to ensure the algorithm $A$ finds
the optimal solution within $\mathcal{O}(1)$ error. For DEKI, in each step we
use $2J$ queries, so that at least $\mathcal{O}(d_{u}/J)$ steps are needed.
This proves the optimality of DEKI.

##  4 Numerical examples

We will apply DEKI to two numerical examples with comparison to vanilla EKI
and LEKI. For fair comparison, we use the same step size of these EKI type
algorithms.

###  4.1 Linear transport equation

Consider the inverse problem of determining the initial wave based on
observations on a later wave field. Here we consider the toy model where the
wave speed is constant. The one-dimensional linear transport equation is

$\frac{\partial U}{\partial t}(x,t)+a\frac{\partial U}{\partial
x}(x,t)=0,\quad x\in[0,1].$

For simplicity, we use the periodic boundary condition $U(t,0)=U(t,1)$. Then
the analytical solution is $U(x,t)=U(\\{x-at\\},0)$, where $\\{x\\}=x-\lfloor
x\rfloor$ is the fractional part of $x$. The analytical solution is used in
the numerical tests.

Suppose we can observe the wave at time $T=1$ at $d_{y}$ equi-spaced points
$x^{D}_{i}=i/d_{y}~{}(i=1,\dots d_{y})$, i.e. the data is generated by

(33) $y=(y_{1},\dots,y_{d_{y}}),\quad y_{i}=U(T,x^{D}_{i})+\xi_{i},$

where $\xi_{i}\sim\mathcal{N}(0,\sigma^{2})$ and $\sigma=10^{-2}$ is the
observational noise. Introduce a computational grid with $d_{u}$ equally
spaced grid points $x_{i}=i/d_{u}$ and denote the restriction of spatial
functions $U(x)$ to this grid by

$u=(U(x_{1}),\dots U(x_{d_{u}}))\in\mathbb{R}^{d_{u}}.$

The forward map is then defined by

(34) $\mathcal{G}:u_{0}\in\mathbb{R}^{d_{u}}\mapsto y\in\mathbb{R}^{d_{y}},$

where $u_{0}=(u(0,x_{1}),\dots,u(0,x_{d_{u}}))\in\mathbb{R}^{d_{u}}$ is the
discrete initial wave to be inferred from the observations Eq. 33.

In numerical experiments, $u_{0}$ is randomly generated
$u_{0}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,I_{d_{u}})$, and we choose
the regularization operator $\mathcal{C}_{0}^{-1/2}=\gamma d_{u}^{-1/2}I$
where $\gamma=0.1$. We will run the DEKI Algorithm 1 with dropout rate
$\lambda=0.5$. Fix reference step size $\tilde{h}=2.5$ and set the ratio
$h_{n}/\tilde{h}_{n}=0.1$, i.e.,

(35)
$\tilde{h}_{n}=\frac{\tilde{h}}{\left\|C_{n}^{uu}\right\|+\epsilon_{0}},\quad
h_{n}=0.1\tilde{h}_{n}.$

Note we add a small parameter $\epsilon_{0}=10^{-12}$ for numerical stability.
The ensemble is initialized as Gaussian variables
$u_{0}^{(j)}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\gamma^{2}I_{d_{u}})$.

We compare DEKI with vanilla EKI and LEKI. For each experiment, we run
$N=10^{2}$ steps and plot the relative loss

(36) $e_{n}=\frac{l(\overline{u}_{n})-l_{\min}}{\left\|y\right\|^{2}},$

where $l_{\min}$ is obtained by optimizing the loss function. We repeat the
experiments $N_{rep}=100$ times and record the mean and deviation of the
relative loss. For each experiment, the data $y$ is fixed and different
realizations of the initial ensemble and dropout are used.

The comparison is shown in Fig. 1. We can see that in the high dimension and
small ensemble regime, the vanilla EKI fails due to the subspace issue, while
DEKI and LEKI work well. We can also see the linear convergence of DEKI.

As discussed before, LEKI needs to design the localization matrix, which
requires prior knowledge (here it is the transport speed). Details of the LEKI
scheme and the localization design are provided in Appendix D. We show that
without such prior knowledge, the LEKI still fails if the localization is not
well designed. We also show the controllable decay of the maximum and minimum
of the diagonal elements of the covariance matrix, confirming our analysis
result.

Figure 1: Comparison of DEKI with vanilla EKI and LEKI in linear transport
equation model with $d_{u}=120,J=20$. Left: decay of the relative data misfit
Eq. 36, here LEKI w/o stands for LEKI without prior knowledge of the transport
speed. Right: collapse of the ensemble, the solid line represents
$\max_{s}C_{n}^{uu}(s,s)$ and the dashed line represents
$\min_{s}C_{n}^{uu}(s,s)$.

We compute the convergence rate by

(37) $r=\frac{\log e_{m}-\log e_{n}}{m-n}.$

In numerical experiments, it is obtained by averaging over $100$ independent
experiments under different problem dimensions $d_{u}$ and ensemble sizes $J$.
For the repeated experiments here, the problem setup (the transport speed
$a$), data, initial ensemble and dropout are all randomized. Fig. 2 shows the
linear dependence of the convergence rate on $d_{u}^{-1}$ and $J$.

Figure 2: Linear dependence of the convergence rate Eq. 37 of DEKI on
reciprocal of dimension $d_{u}^{-1}$ and ensemble size $J$, tested in linear
transport equation model.

###  4.2 Darcy’s Law

Consider the 2d elliptic equation on $\Omega=[0,1]^{2}$:

(38) $-\nabla\cdot\left(\exp(\kappa(x))\nabla u(x)\right)=f(x),\quad
u|_{\partial\Omega}=0.$

Consider estimating the log-permeability field $\kappa(x)$ from observations
on the field $u(x)$: pick evaluation points $x_{1},\dots,x_{d_{y}}$, and

$y=(y_{1},\dots,y_{d_{y}}),\quad y_{j}=u(x_{j})+\xi_{j},$

where $u(x)$ is the solution of Eq. 38 given parameter $\kappa(x)$ and
$\xi_{j}\sim\mathcal{N}(0,\sigma_{0}^{2})~{}(\sigma_{0}=10^{-3})$ is the
observational noise. Thus the forward map is

(39) $\mathcal{G}:\kappa(x)\mapsto y.$

The testing parameter $\kappa(x)$ is generated from a Gaussian random field
with mean $m$ and covariance operator

$k(x_{1},y_{1};x_{2},y_{2})=\sigma^{2}\exp\left(-\frac{(x_{1}-x_{2})^{2}}{2l_{x}^{2}}-\frac{(y_{1}-y_{2})^{2}}{2l_{y}^{2}}\right).$

In computation, we represent $\kappa$ in the truncated Karhunen–Loéve basis:

$\kappa(x)\approx
m+\sum_{n=1}^{d_{u}}\kappa_{n}\varphi_{n}(x),\quad\varphi_{n}(x)=\sqrt{\lambda_{n}}\psi_{n}(x).$

where $(\lambda_{n},\psi_{n}(x))$ are the eigenpairs of the covariance
operator $k$ arranged in descending order
$\lambda_{1}\geq\lambda_{2}\geq\cdots$. Denote the coordinate representation
$\hat{\kappa}=(\kappa_{1},\dots,\kappa_{d_{u}})\in\mathbb{R}^{d_{u}}$. The
truncated dimension $d_{u}$ is determined by a threshold $\epsilon$:

$d_{u}(\epsilon)=\inf\left\\{k\in\mathbb{Z}:\sum_{1\leq n\leq
k}\lambda_{n}\geq(1-\epsilon)\text{tr}(k)\right\\}.$

Namely, $d_{u}$ is the smallest number of basis functions that can approximate
$\kappa$ with $\epsilon$-accuracy. The regularization term under this basis
can be simply chosen as

$R(\hat{\kappa})=\gamma^{2}\left\|\hat{\kappa}\right\|_{l^{2}}^{2}=\gamma^{2}\sum_{i=1}^{d_{u}}\kappa_{i}^{2}.$

Or equivalently, the regularization operator $\mathcal{C}_{0}^{-1/2}=\gamma
I$.

We test under different setups listed in Table 1, and fix external force
$\displaystyle f(x,y)=13\pi^{2}\sin(2\pi x)\sin(3\pi y)$. The forward problem
is solved using finite element method on a $32\times 32$ grid. The observation
points are chosen as the center of the $8\times 8$ subblocks of $\Omega$, so
that the data dimension $d_{y}=64$.

Table 1: Parameters for different setups in Darcy’s law model | $m$ | $\sigma$ | $l_{x}$ | $l_{y}$ | $\varepsilon$ | $d_{u}$
---|---|---|---|---|---|---
setup 1 | $0$ | $0.1$ | $0.1$ | $0.1$ | $10^{-3}$ | $136$
setup 2 | $0$ | $0.1$ | $0.2$ | $0.05$ | $10^{-3}$ | $142$
setup 3 | $0$ | $0.1$ | $0.15$ | $0.05$ | $10^{-3}$ | $179$
setup 4 | $0$ | $0.1$ | $0.1$ | $0.05$ | $10^{-3}$ | $254$

We will compare DEKI with EKI. We do not use LEKI here, since it is not clear
how to design localization under the KL basis. The ensemble is initialized as
Gaussian variables
$\hat{\kappa}_{0}^{(j)}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\gamma^{2}I_{d_{u}})$,
and choose the regularization parameter $\gamma=0.1$. Fix reference step size
$\tilde{h}=0.5$ and set the ratio $h_{n}/\tilde{h}_{n}=0.1$, see Eq. 35.

In each experiment, we run $N=10^{4}$ steps to show the asymptotic behavior.
Note in practice, such long steps are not necessary. In Fig. 3, we compare the
decay of the relative data misfit Eq. 36 and the relative solution error of
DEKI and EKI in setup 1. Here the relative solution error is defined by

(40) $\text{err}_{n}=\frac{\left\|\kappa_{n}-\kappa_{\rm
true}\right\|_{L^{2}}}{\left\|\kappa_{\rm true}\right\|_{L^{2}}}.$

We repeat each experiment $50$ times under the same settings and plot the mean
and standard deviation of these error. Here different realizations of initial
ensemble and dropout are used for the repeated experiments.

Again we see that the EKI gets stuck after some steps but DEKI works well.
Here the data misfit of DEKI does not decay exponentially, since the PL
condition is not satisfied for this model. Nevertheless, the DEKI still
performs well both as an optimizer and as an inversion method in the high
dimensional settings.

Figure 3: Comparison of DEKI and EKI in Darcy’s law model in setup 1
($d_{u}=136$) using ensemble size $J=15$. Left: relative data misfit Eq. 36;
right: relative solution error Eq. 40.

Plots of the inversion solutions are displayed in Fig. 4. We can see that as
dimension grows, DEKI is still able to infer the permeability field with good
accuracy using small ensemble, while EKI solutions become incorrect.

Figure 4: Comparison of the inversion solutions of EKI and DEKI in Darcy’s law
model in different setups with ensemble size $J=15$. Row 1: setup 1
$(l_{x}=0.1,l_{y}=0.1,d_{u}=136)$; Row 2: setup 2
$(l_{x}=0.2,l_{y}=0.05,d_{u}=142)$; Row 3: setup 3
$(l_{x}=0.15,l_{y}=0.05,d_{u}=179)$; Row 4: setup 4
$(l_{x}=0.1,l_{y}=0.05,d_{u}=254)$.

##  5 Conclusion

We investigate the application of the dropout technique in EKI method and
propose the DEKI scheme. We show that this method mitigates the subspace issue
of the vanilla EKI, and performs well in high dimensional inverse problems
using small ensemble, both theoretically and numerically. With careful
analysis on the dimension dependence, we prove that DEKI converges
exponentially for strongly convex problems and the computational complexity
scales linearly with dimension. We also show that such scaling is optimal by
conducting the query complexity analysis adapting from the optimization
literature, which is new for inverse problems.

The idea of incorporating dropout technique into EKI follows the trial of
EnKBF with dropout in [25]. The numerical analysis here is a first step to
understand the effectiveness of dropout in these methods. It would be
interesting to investigate if the idea can be extended and work as a simpler
but effective technique for other ensemble methods in high dimensional
problems. It is also interesting to look at the interplay of our DEKI method
with traditional zeroth order optimization methods, providing new perspective
to understand them and generating new algorithms.

## Acknowledgments

The work of SL has been partially funded by Singapore MOE grant
A-8000459-00-00. The work of SR has been partially funded by Deutsche
Forschungsgemeinschaft (DFG) - Project-ID 318763901 - SFB1294. SR would also
like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge,
for support and hospitality during the programme The Mathematical and
Statistical Foundation of Future Data-Driven Engineering where work on this
paper was undertaken. This work was supported by EPSRC grant no EP/R014604/1.
The work of XT has been funded by Singapore MOE grants A-0004263-00-00 and
A-8000459-00-00.

## References

  * [1] A. Agarwal, M. J. Wainwright, P. Bartlett, and P. Ravikumar, Information-theoretic lower bounds on the oracle complexity of convex optimization, in Advances in Neural Information Processing Systems, Y. Bengio, D. Schuurmans, J. Lafferty, C. Williams, and A. Culotta, eds., vol. 22, Curran Associates, Inc., 2009.
  * [2] M. Benning and M. Burger, Modern regularization methods for inverse problems, Acta Numer., 27 (2018), pp. 1–111, https://doi.org/10.1017/s0962492918000016.
  * [3] D. Blömker, C. Schillings, P. Wacker, and S. Weissmann, Well posedness and convergence analysis of the ensemble Kalman inversion, Inverse Problems, 35 (2019), pp. 085007, 32, https://doi.org/10.1088/1361-6420/ab149c.
  * [4] E. Calvello, S. Reich, and A. M. Stuart, Ensemble Kalman methods: A mean field perspective, arXiv preprint arXiv:2209.11371, (2022).
  * [5] Y. Carmon, J. C. Duchi, O. Hinder, and A. Sidford, Lower bounds for finding stationary points I, Math. Program., 184 (2020), pp. 71–120, https://doi.org/10.1007/s10107-019-01406-y.
  * [6] N. K. Chada, M. A. Iglesias, L. Roininen, and A. M. Stuart, Parameterizations for ensemble Kalman inversion, Inverse Problems, 34 (2018), p. 055009.
  * [7] N. K. Chada, A. M. Stuart, and X. T. Tong, Tikhonov regularization within ensemble Kalman inversion, SIAM J. Numer. Anal., 58 (2020), pp. 1263–1294, https://doi.org/10.1137/19M1242331.
  * [8] N. K. Chada and X. T. Tong, Convergence acceleration of ensemble Kalman inversion in nonlinear settings, Math. Comp., 91 (2022), pp. 1247–1280, https://doi.org/10.1090/mcom/3709.
  * [9] H. W. Engl, M. Hanke, and A. Neubauer, Regularization of inverse problems, vol. 375 of Mathematics and its Applications, Kluwer Academic Publishers Group, Dordrecht, 1996.
  * [10] Y. Gal and Z. Ghahramani, Dropout as a Bayesian approximation: Representing model uncertainty in deep learning, in Proceedings of The 33rd International Conference on Machine Learning, vol. 48, New York, USA, 2016, PMLR, pp. 1050–1059.
  * [11] T. M. Hamill, J. S. Whitaker, and C. Snyder, Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter, Monthly Weather Review, 129 (2001), pp. 2776 – 2790, https://doi.org/https://doi.org/10.1175/1520-0493(2001)129<2776:DDFOBE>2.0.CO;2.
  * [12] R. A. Horn and C. R. Johnson, Matrix analysis, Cambridge University Press, Cambridge, second ed., 2013.
  * [13] P. L. Houtekamer and H. L. Mitchell, Data Assimilation Using an Ensemble Kalman Filter Technique, Monthly Weather Review, 126 (1998), p. 796, https://doi.org/10.1175/1520-0493(1998)126<0796:DAUAEK>2.0.CO;2.
  * [14] P. L. Houtekamer and H. L. Mitchell, A Sequential Ensemble Kalman Filter for Atmospheric Data Assimilation, Monthly Weather Review, 129 (2001), p. 123, https://doi.org/10.1175/1520-0493(2001)129<0123:ASEKFF>2.0.CO;2.
  * [15] M. A. Iglesias, K. J. H. Law, and A. M. Stuart, Ensemble Kalman methods for inverse problems, Inverse Problems, 29 (2013), pp. 045001, 20, https://doi.org/10.1088/0266-5611/29/4/045001.
  * [16] D. P. Kingma, T. Salimans, and M. Welling, Variational dropout and the local reparameterization trick, in Advances in Neural Information Processing Systems, vol. 28, Curran Associates, Inc., 2015, https://proceedings.neurips.cc/paper_files/paper/2015/file/bc7316929fe1545bf0b98d114ee3ecb8-Paper.pdf.
  * [17] A. Kirsch, An introduction to the mathematical theory of inverse problems, vol. 120 of Applied Mathematical Sciences, Springer, Cham, [2021] ©2021, https://doi.org/10.1007/978-3-030-63343-1.  Third edition [of 1479408].
  * [18] N. B. Kovachki and A. M. Stuart, Ensemble Kalman inversion: A derivative-free technique for machine learning tasks, Inverse Problems, 35 (2019), pp. 095005, 35, https://doi.org/10.1088/1361-6420/ab1c3a.
  * [19] C. Liu, L. Zhu, and M. Belkin, Loss landscapes and optimization in over-parameterized non-linear systems and neural networks, Appl. Comput. Harmon. Anal., 59 (2022), pp. 85–116, https://doi.org/10.1016/j.acha.2021.12.009.
  * [20] A. S. Nemirovsky and D. B. a. Yudin, Problem complexity and method efficiency in optimization, Wiley-Interscience Series in Discrete Mathematics, John Wiley & Sons, Inc., New York, 1983.  Translated from the Russian and with a preface by E. R. Dawson.
  * [21] Y. Nesterov, Introductory lectures on convex optimization, vol. 87 of Applied Optimization, Kluwer Academic Publishers, Boston, MA, 2004, https://doi.org/10.1007/978-1-4419-8853-9.  A basic course.
  * [22] Y. Nesterov and V. Spokoiny, Random gradient-free minimization of convex functions, Found. Comput. Math., 17 (2017), pp. 527–566, https://doi.org/10.1007/s10208-015-9296-2.
  * [23] E. Ott, B. R. Hunt, I. Szunyogh, A. V. Zimin, E. J. Kostelich, M. Corazza, E. Kalnay, D. Patil, and J. A. Yorke, A local ensemble Kalman filter for atmospheric data assimilation, Tellus A: Dynamic Meteorology and Oceanography, 56 (2004), pp. 415–428, https://doi.org/10.3402/tellusa.v56i5.14462.
  * [24] F. Parzer and O. Scherzer, On convergence rates of adaptive ensemble Kalman inversion for linear ill-posed problems, Numer. Math., 152 (2022), pp. 371–409, https://doi.org/10.1007/s00211-022-01314-y.
  * [25] J. Pidstrigach and S. Reich, Affine-invariant ensemble transform methods for logistic regression, Found. Comput. Math., 23 (2023), pp. 675–708, https://doi.org/10.1007/s10208-022-09550-2.
  * [26] C. Schillings and A. M. Stuart, Analysis of the ensemble Kalman filter for inverse problems, SIAM J. Numer. Anal., 55 (2017), pp. 1264–1290, https://doi.org/10.1137/16M105959X.
  * [27] C. Schillings and A. M. Stuart, Convergence analysis of ensemble Kalman inversion: The linear, noisy case, Appl. Anal., 97 (2018), pp. 107–123, https://doi.org/10.1080/00036811.2017.1386784.
  * [28] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, Dropout: A simple way to prevent neural networks from overfitting, Journal of Machine Learning Research, 15 (2014), pp. 1929–1958, http://jmlr.org/papers/v15/srivastava14a.html.
  * [29] A. M. Stuart, Inverse problems: A Bayesian perspective, Acta Numer., 19 (2010), pp. 451–559, https://doi.org/10.1017/S0962492910000061.
  * [30] X. T. Tong and M. Morzfeld, Localized ensemble Kalman inversion, Inverse Problems, 39 (2023), pp. Paper No. 064002, 38.
  * [31] S. Wager, S. Wang, and P. Liang, Dropout training as adaptive regularization, in Proceedings of the 26th International Conference on Neural Information Processing Systems - Volume 1, NIPS’13, Red Hook, NY, USA, 2013, Curran Associates Inc., p. 351–359.
  * [32] L. Wan, M. Zeiler, S. Zhang, Y. Le Cun, and R. Fergus, Regularization of neural networks using DropConnect, in Proceedings of the 30th International Conference on Machine Learning, vol. 28, Atlanta, Georgia, USA, 2013, PMLR, pp. 1058–1066.
  * [33] S. Weissmann, Gradient flow structure and convergence analysis of the ensemble Kalman inversion for nonlinear forward models, Inverse Problems, 38 (2022), pp. Paper No. 105011, 30.

##  Appendix A Missing proofs

###  A.1 Proof for Proposition 3.5

Denote for simplicity that

$\begin{split}&T_{n}=\frac{1}{\sqrt{J-1}}[\tilde{\tau}_{n}^{(1)},\dots,\tilde{\tau}_{n}^{(J)}],\quad\hat{Y}_{n}=G_{n}T_{n},\\\
&Y_{n}=\frac{1}{\sqrt{J-1}}[\tilde{\mathfrak{r}}_{n}^{(1)},\dots,\tilde{\mathfrak{r}}_{n}^{(J)}],\quad\tilde{\mathfrak{r}}_{n}^{(j)}=\mathcal{H}(\tilde{u}_{n}^{(j)})-\overline{\mathcal{H}(\tilde{u}_{n})}.\end{split}$

Then it holds $\tilde{C}_{n}^{uu}=T_{n}T_{n}^{\rm
T},\tilde{C}_{n}^{up}=T_{n}Y_{n}^{\rm T},\tilde{C}_{n}^{pp}=Y_{n}Y_{n}^{\rm
T}$, and we can write Eq. 10 as

(41) $\overline{u}_{n+1}=\overline{u}_{n}+\tilde{h}_{n}T_{n}Y_{n}^{\rm
T}(I+\tilde{h}_{n}Y_{n}Y_{n}^{\rm T})^{-1}(z-\mathcal{H}(\overline{u}_{n})),$

and similarly for Eq. 24:

(42)
$\overline{u}_{n+1}^{\prime}=\overline{u}_{n}+\tilde{h}_{n}T_{n}\hat{Y}_{n}^{\rm
T}(I+\tilde{h}_{n}\hat{Y}_{n}^{\rm
T}\hat{Y}_{n})^{-1}(z-\mathcal{H}(\overline{u}_{n})).$

Denote $Z_{t}=(1-t)Y_{n}+t\hat{Y}_{n}$, and define

$S(t)=Z_{t}^{\rm T}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm T})^{-1},$

Take difference of Eq. 41 and Eq. 42,

(43)
$\begin{split}\|\overline{u}_{n+1}-\overline{u}_{n+1}^{\prime}\|=~{}&\|\tilde{h}_{n}T_{n}(S(1)-S(0))(z-\mathcal{H}(\overline{u}_{n})\|\\\
\leq~{}&\tilde{h}_{n}\|T_{n}\|\|S(1)-S(0)\|\|z-\mathcal{H}(\overline{u}_{n})\|.\end{split}$

By direct computation and the Sherman-Morrison-Woodbury formula (Lemma B.1),

$\begin{split}S^{\prime}(t)=~{}&(Z_{t}^{\rm
T})^{\prime}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm T})^{-1}\\\
&-\tilde{h}_{n}Z_{t}^{\rm T}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm
T})^{-1}[Z_{t}^{\prime}Z_{t}^{\rm T}+Z_{t}(Z_{t}^{\rm
T})^{\prime}](I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm T})^{-1}\\\
=~{}&\left[I-\tilde{h}_{n}Z_{t}^{\rm T}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm
T})^{-1}Z_{t}\right](Z_{t}^{\rm T})^{\prime}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm
T})^{-1}\\\ &-\tilde{h}_{n}Z_{t}^{\rm T}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm
T})^{-1}Z_{t}^{\prime}Z_{t}^{\rm T}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm T})^{-1}\\\
=~{}&(I+\tilde{h}_{n}Z_{t}^{\rm T}Z_{t})^{-1}(Z_{t}^{\rm
T})^{\prime}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm T})^{-1}\\\
&-\tilde{h}_{n}Z_{t}^{\rm T}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm
T})^{-1}Z_{t}^{\prime}Z_{t}^{\rm T}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm
T})^{-1}.\end{split}$

Notice since $Z_{t}Z_{t}^{\rm T}$ is PSD, apply Lemma B.2 and we get

$\begin{split}&\|Z_{t}^{\rm T}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm T})^{-1}\|\\\
=~{}&\|(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm T})^{-1}Z_{t}Z_{t}^{\rm
T}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm T})^{-1}\|^{1/2}\\\
\leq~{}&\|(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm
T})^{-1}\|^{1/2}\tilde{h}_{n}^{-1/2}\|\tilde{h}_{n}Z_{t}Z_{t}^{\rm
T}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm
T})^{-1}\|^{1/2}\leq\tilde{h}_{n}^{-1/2}.\end{split}$

Therefore,

$\begin{split}\left\|S^{\prime}(t)\right\|\leq~{}&\|(I+\tilde{h}_{n}Z_{t}^{\rm
T}Z_{t})^{-1}\|\left\|(Z_{t}^{\rm
T})^{\prime}\right\|\|(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm T})^{-1}\|\\\
&+\tilde{h}_{n}\|Z_{t}^{\rm T}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm
T})^{-1}\|\|Z_{t}^{\prime}\|\|Z_{t}^{\rm T}(I+\tilde{h}_{n}Z_{t}Z_{t}^{\rm
T})^{-1}\|\\\
\leq~{}&\|Z_{t}^{\prime}\|+\tilde{h}_{n}(\tilde{h}_{n}^{-1/2})^{2}\left\|Z_{t}^{\prime}\right\|=2\|Y_{n}-\hat{Y}_{n}\|.\end{split}$

By definition,

(44)
$\begin{split}Y_{n}-\hat{Y}_{n}=~{}&\frac{1}{\sqrt{J-1}}[\tilde{\mathfrak{r}}_{n}^{(1)}-G_{n}\tilde{\tau}_{n}^{(1)},\dots,\tilde{\mathfrak{r}}_{n}^{(J)}-G_{n}\tilde{\tau}_{n}^{(J)}]\\\
=~{}&\frac{1}{\sqrt{J-1}}[r_{n}^{(1)},\dots,r_{n}^{(J)}]\left(I-\frac{1}{J}\mathbf{1}\mathbf{1}^{\rm
T}\right),\end{split}$

where $\mathbf{1}=(1,1,\dots,1)^{\rm T}$ and

$r_{n}^{(j)}=\mathcal{H}(\tilde{u}_{n}^{(j)})-\mathcal{H}(\overline{u}_{n})-G_{n}\tilde{\tau}_{n}^{(j)}=\begin{bmatrix}\mathcal{G}(\tilde{u}_{n}^{(j)})-\mathcal{G}(\overline{u}_{n})-\nabla\mathcal{G}(\overline{u}_{n})\tilde{\tau}_{n}^{(j)}\\\
0\end{bmatrix}.$

Using the bounds $\left\|\nabla^{2}\mathcal{G}\right\|\leq H$ in Eq. 13, we
obtain that

$\|r_{n}^{(j)}\|\leq\frac{1}{2}H\|\tilde{\tau}_{n}^{(j)}\|^{2}.$

Notice $I-\frac{1}{J}\mathbf{1}\mathbf{1}^{\rm T}$ is an orthogonal projector
and thus from Eq. 44,

$\begin{split}\|Y_{n}-\hat{Y}_{n}\|^{2}\leq~{}&\frac{1}{J-1}\|[r_{n}^{(1)},\dots,r_{n}^{(J)}]\|^{2}\\\
\leq~{}&\frac{1}{J-1}\sum_{j=1}^{J}\|r_{n}^{(j)}\|^{2}\leq\frac{H^{2}}{4(J-1)}\sum_{j=1}^{J}\|\tilde{\tau}_{n}^{(j)}\|^{4}.\end{split}$

Therefore,

$\begin{split}\|Y_{n}-\hat{Y}_{n}\|\leq~{}&\frac{H}{2\sqrt{J-1}}\left(\sum_{j=1}^{J}\|\tilde{\tau}_{n}^{(j)}\|^{4}\right)^{1/2}\\\
\leq~{}&\frac{H}{2\sqrt{J-1}}\sum_{j=1}^{J}\|\tilde{\tau}_{n}^{(j)}\|^{2}\\\
=~{}&\frac{H}{2}\sqrt{J-1}\text{tr}(\tilde{C}_{n}^{uu})\leq\frac{1}{2}(J-1)^{3/2}H\|\tilde{C}_{n}^{uu}\|.\end{split}$

The last inequality holds since $\text{rank}(\tilde{C}_{n}^{uu})\leq J-1$.
Combining the above results into Eq. 43 and noticing
$\|\tilde{C}_{n}^{uu}\|=\|T_{n}T_{n}^{\rm T}\|=\|T_{n}\|^{2}$, we obtain

$\begin{split}\|\overline{u}_{n+1}-\overline{u}_{n+1}^{\prime}\|\leq~{}&\tilde{h}_{n}\|T_{n}\|\|z-\mathcal{H}(\overline{u}_{n})\|\int_{0}^{1}\|S^{\prime}(t)\|\mathrm{d}t\\\
\leq~{}&\tilde{h}_{n}\|\tilde{C}_{n}^{uu}\|^{1/2}\|z-\mathcal{H}(\overline{u}_{n})\|\cdot(J-1)^{3/2}H\|\tilde{C}_{n}^{uu}\|.\end{split}$

Thus the conclusion holds for $C=(J-1)^{3/2}H$.

###  A.2 Proof for Lemma 3.8

Since $l$ is $L$-smooth,

(45)
$\begin{split}l(\overline{u}_{n+1})\leq~{}&l(\overline{u}_{n})+\left\langle\nabla
l(\overline{u}_{n}),\overline{u}_{n+1}-\overline{u}_{n}\right\rangle+\frac{1}{2}L\left\|\overline{u}_{n+1}-\overline{u}_{n}\right\|^{2}\\\
=~{}&l(\overline{u}_{n})+\left\langle\nabla
l(\overline{u}_{n}),\overline{u}^{\prime}_{n+1}-\overline{u}_{n}\right\rangle+\frac{1}{2}L\left\|\overline{u}_{n+1}^{\prime}-\overline{u}_{n}\right\|^{2}\\\
+&\left\langle\nabla
l(\overline{u}_{n})+L(\overline{u}_{n+1}^{\prime}-\overline{u}_{n}),\overline{u}_{n+1}-\overline{u}^{\prime}_{n+1}\right\rangle+\frac{1}{2}L\|\overline{u}_{n+1}-\overline{u}^{\prime}_{n+1}\|^{2}.\end{split}$

where $\overline{u}_{n+1}^{\prime}$ is defined in Eq. 24. The second term
provides contraction, and the rest terms are the residues to be controlled.
Note though the third term is a higher order term
$(\mathcal{O}(\tilde{h}_{n}^{2}))$, it is not small due to the adaptive step
size. Thus we will combine the second and third terms together to derive
sufficient decay. And the last two terms is the nonlinear residue, which will
be controlled using Proposition 3.5.

By definition, we can write the second term in Eq. 45 as

(46) $\begin{split}&\left\langle\nabla
l(\overline{u}_{n}),\overline{u}_{n+1}^{\prime}-\overline{u}_{n}\right\rangle\\\
=&-\left\langle G_{n}^{\rm
T}(z-\mathcal{H}(\overline{u}_{n})),\tilde{h}_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T}(I+\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T})^{-1}(z-\mathcal{H}(\overline{u}_{n}))\right\rangle\\\ =&-\left\langle
z-\mathcal{H}(\overline{u}_{n}),\psi_{1}(\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T})(z-\mathcal{H}(\overline{u}_{n}))\right\rangle,\end{split}$

where $\psi_{1}(x)=x(1+x)^{-1}$. For the third term in Eq. 45,

(47)
$\begin{split}&\frac{1}{2}L\left\|\overline{u}_{n+1}^{\prime}-\overline{u}_{n}\right\|^{2}\\\
=~{}&\frac{1}{2}L\|\tilde{h}_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T}(I+\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T})^{-1}(z-\mathcal{H}(\overline{u}_{n}))\|^{2}\\\
\leq~{}&\frac{1}{2}L\tilde{h}_{n}\|\tilde{C}_{n}^{uu}\|\|(\tilde{h}_{n}\tilde{C}_{n}^{uu})^{1/2}G_{n}^{\rm
T}(I+\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T})^{-1}(z-\mathcal{H}(\overline{u}_{n}))\|^{2}\\\
=~{}&\frac{1}{2}L\tilde{h}_{n}\|\tilde{C}_{n}^{uu}\|\left\langle
z-\mathcal{H}(\overline{u}_{n}),\psi_{2}(\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T})(z-\mathcal{H}(\overline{u}_{n}))\right\rangle,\end{split}$

where $\psi_{2}(x)=x(1+x)^{-2}$. Notice $\tilde{C}_{n}^{uu}$ is a submatrix of
$C_{n}^{uu}\Rightarrow\|\tilde{C}_{n}^{uu}\|\leq\|C_{n}^{uu}\|$. By assumption
$\mu\leq L^{-1}$, and thus the choice $\tilde{h}_{n}=\mu\|C_{n}^{uu}\|^{-1}$
implies that

$\frac{1}{2}L\tilde{h}_{n}\|\tilde{C}_{n}^{uu}\|\leq\frac{1}{2}L\mu\leq\frac{1}{2}.$

Therefore, combine Eq. 46 and Eq. 47 and we obtain

(48) $\begin{split}&\left\langle\nabla
l(\overline{u}_{n}),\overline{u}_{n+1}^{\prime}-\overline{u}_{n}\right\rangle+\frac{1}{2}L\left\|\overline{u}_{n+1}^{\prime}-\overline{u}_{n}\right\|^{2}\\\
\leq&-\left\langle
z-\mathcal{H}(\overline{u}_{n}),\psi(\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T})(z-\mathcal{H}(\overline{u}_{n}))\right\rangle,\end{split}$

where

$\psi(x)=\psi_{1}(x)-\frac{1}{2}\psi_{2}(x)=\frac{x}{1+x}-\frac{x}{2(1+x)^{2}}=\frac{x(1+2x)}{2(1+x)^{2}},$

Notice the upper bound holds

$\|\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T}\|\leq\tilde{h}_{n}M^{2}\|\tilde{C}_{n}^{uu}\|\leq\mu M^{2}.$

It is elementary to show that $\psi(x)/x$ is monotone decreasing, and thus
when $0<x\leq\mu M^{2}$, it holds $\psi(x)\geq\mu^{-1}M^{-2}\psi(\mu M^{2})x$.
Therefore,

$\displaystyle\mathbb{E}_{n}\psi(\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T})$ $\displaystyle\succeq~{}$ $\displaystyle\mu^{-1}M^{-2}\psi(\mu
M^{2})\cdot\mathbb{E}_{n}\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm T}$
(Lemma B.2) $\displaystyle\succeq~{}$ $\displaystyle\mu^{-1}M^{-2}\psi(\mu
M^{2})\cdot\lambda(1-\lambda)\tilde{h}_{n}\min_{s}C_{n}^{uu}(s,s)G_{n}G_{n}^{\rm
T}$ by Eq. 23 $\displaystyle\succeq~{}$ $\displaystyle\mu^{-1}M^{-2}\psi(\mu
M^{2})\cdot\lambda(1-\lambda)\tilde{h}_{n}\bar{\kappa}^{-1}\|C_{n}^{uu}\|\min_{s}P(s,s)G_{n}G_{n}^{\rm
T}$ (Proposition 3.3) $\displaystyle=~{}$ $\displaystyle M^{-2}\psi(\mu
M^{2})\lambda(1-\lambda)\bar{\kappa}^{-1}\min_{s}P(s,s)G_{n}G_{n}^{\rm T}.$

Plug it in Eq. 48 and use the PL condition $\|\nabla l(u)\|^{2}\geq
c(l(u)-l_{\min})$, we obtain

(49) $\begin{split}&\mathbb{E}_{n}\left[\left\langle\nabla
l(\overline{u}_{n}),\overline{u}_{n+1}^{\prime}-\overline{u}_{n}\right\rangle+\frac{1}{2}L\left\|\overline{u}_{n+1}^{\prime}-\overline{u}_{n}\right\|^{2}\right]\\\
\leq&-\mathbb{E}_{n}\left\langle
z-\mathcal{H}(\overline{u}_{n}),\psi(\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T})(z-\mathcal{H}(\overline{u}_{n}))\right\rangle\\\ \leq&-M^{-2}\psi(\mu
M^{2})\lambda(1-\lambda)\bar{\kappa}^{-1}\min_{s}P(s,s)\|G_{n}^{\rm
T}(z-\mathcal{H}(\overline{u}_{n}))\|^{2}\\\
\leq&-2\beta(l(\bar{u}_{n})-l_{\min}),\end{split}$

where we denote $\beta=\frac{1}{2}cM^{-2}\psi(\mu
M^{2})\lambda(1-\lambda)\bar{\kappa}^{-1}\min_{s}P(s,s)$.

Next we control the nonlinear residues. For the fourth term in Eq. 45, first
notice

(50)
$\begin{split}\|\overline{u}_{n+1}^{\prime}-\overline{u}_{n}\|\leq~{}&\|\tilde{h}_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T}(I+\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T})^{-1}\|\|z-\mathcal{H}(\overline{u}_{n})\|\\\
\leq~{}&\mu^{1/2}\|(\tilde{h}_{n}\tilde{C}_{n}^{uu})^{1/2}G_{n}^{\rm
T}(I+\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T})^{-1}\|\|z-\mathcal{H}(\overline{u}_{n})\|\\\
\leq~{}&\mu^{1/2}\|z-\mathcal{H}(\overline{u}_{n})\|,\end{split}$

where we uses the following inequality: since
$G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm T}\succeq 0$,

$\begin{split}&\|(\tilde{h}_{n}\tilde{C}_{n}^{uu})^{1/2}G_{n}^{\rm
T}(I+\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm T})^{-1}\|^{2}\\\
=~{}&\|(I+\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T})^{-1}\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T}(I+\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm T})^{-1}\|\\\
\leq~{}&\|(I+\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T})^{-1}\|\|\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm
T}(I+\tilde{h}_{n}G_{n}\tilde{C}_{n}^{uu}G_{n}^{\rm T})^{-1}\|\leq
1.\end{split}$

Using Eq. 50, Eq. 25 and $\tilde{h}_{n}\leq\mu\|\tilde{C}_{n}^{uu}\|^{-1}$, we
get

$\begin{split}&\left\langle\nabla
l(\overline{u}_{n})+L(\overline{u}_{n+1}^{\prime}-\overline{u}_{n}),\overline{u}_{n+1}-\overline{u}^{\prime}_{n+1}\right\rangle\\\
\leq~{}&\left(\|G^{\rm
T}_{n}(z-\mathcal{H}(\overline{u}_{n}))\|+L\|\overline{u}_{n+1}^{\prime}-\overline{u}_{n}\|\right)\|\overline{u}_{n+1}-\overline{u}^{\prime}_{n+1}\|\\\
\leq~{}&C(M+L\mu^{1/2})\tilde{h}_{n}\|\tilde{C}_{n}^{uu}\|^{3/2}\|z-\mathcal{H}(\overline{u}_{n})\|^{2}\\\
\leq~{}&C(M+L\mu^{1/2})\mu\|\tilde{C}_{n}^{uu}\|^{1/2}\|z-\mathcal{H}(\overline{u}_{n})\|^{2}.\end{split}$

The last term can be controlled in a similar way:

$\begin{split}\frac{1}{2}L\|\overline{u}_{n+1}-\overline{u}^{\prime}_{n+1}\|^{2}\leq~{}&\frac{1}{2}LC^{2}\tilde{h}_{n}^{2}\|\tilde{C}_{n}^{uu}\|^{3}\|z-\mathcal{H}(\overline{u}_{n})\|^{2}\\\
\leq~{}&\frac{1}{2}LC^{2}\mu^{2}\|\tilde{C}_{n}^{uu}\|\|z-\mathcal{H}(\overline{u}_{n})\|^{2}.\end{split}$

Summing up the above two inequalities, we obtain

(51) $\begin{split}\left\langle\nabla
l(\overline{u}_{n})+L(\overline{u}_{n+1}^{\prime}-\overline{u}_{n}),\overline{u}_{n+1}-\overline{u}^{\prime}_{n+1}\right\rangle+\frac{1}{2}L\|\overline{u}_{n+1}-\overline{u}^{\prime}_{n+1}\|^{2}&\\\
\leq\frac{1}{2}\Delta_{n}\|z-\mathcal{H}(\overline{u}_{n})\|^{2}=\Delta_{n}l(\overline{u}_{n})&,\end{split}$

where
$\Delta_{n}=2C\mu(M+L\mu^{1/2})\|C_{n}^{uu}\|^{1/2}+LC^{2}\mu^{2}\|C_{n}^{uu}\|$.
Note here we use the relation $\|\tilde{C}_{n}^{uu}\|\leq\|C_{n}^{uu}\|$.

Combine Eq. 45, Eq. 49 and Eq. 51, and we prove our claim.

###  A.3 Proof for Proposition 3.9

We will prove the lower bound by construction. For simplicity, normalize
$\left\|y\right\|=1$. We start from some arbitrary
$\left\|G^{\prime}\right\|\leq 1$, and consider the matrix class

$P=\\{\left\|G\right\|\leq 2:GU=G^{\prime}U\\}.$

Note $\forall G\in P,A(U,GU,y)\equiv A(U,G^{\prime}U,y)$. We will construct
$G_{0},G_{B}\in P$ s.t.

(52) $\|u^{*}(G_{0})-u^{*}(G_{B})\|\geq 0.2.$

Then at least one $G_{*}\in\\{G_{0},G_{B}\\}$ would satisfy
$\|A(U,G_{*}U,y)-u^{*}(G_{*})\|\geq 0.1$, since

$\begin{split}&\|A(U,G_{0}U,y)-u^{*}(G_{0})\|+\|A(U,G_{B}U,y)-u^{*}(G_{B})\|\\\
\geq~{}&\|\left(A(U,G_{0}U,y)-u^{*}(G_{0})\right)-\left(A(U,G_{B}U,y)-u^{*}(G_{B})\right)\|\\\
=~{}&\|u^{*}(G_{0})-u^{*}(G_{B})\|\geq 0.2.\end{split}$

And the conclusion follows.

The proof is divided into two parts. First we express $P$ and $u^{*}(G)$ in an
explicit form. Then we construct $G_{0},G_{B}\in P$ that satisfy Eq. 52.

Explicit form of $P$ and $u^{*}(G)$. $U$ admits a singular value decomposition

$U=\Sigma\begin{pmatrix}\Lambda&0\\\
0&0\end{pmatrix}V=\Sigma\begin{pmatrix}I_{r}&0\\\
0&0\end{pmatrix}\begin{pmatrix}\Lambda&0\\\ 0&I\end{pmatrix}V,$

where $r$ is the rank of $U$, and notice $r\leq n\leq d_{u}/2$. For any $G\in
P$, denote

$G\Sigma=(A,B),~{}G^{\prime}\Sigma=(C,D),\quad A,C\in\mathbb{R}^{d_{y}\times
r},~{}B,D\in\mathbb{R}^{d_{y}\times(d_{u}-r)}.$

Denote $R=\begin{pmatrix}\Lambda&0\\\ 0&I\end{pmatrix}V$, and notice

$\begin{split}GU=G^{\prime}U\Leftrightarrow~{}&G\Sigma\Sigma^{\rm
T}UR^{-1}=G^{\prime}\Sigma\Sigma^{\rm T}UR^{-1}\\\
\Leftrightarrow~{}&(A,B)\begin{pmatrix}I_{r}&0\\\
0&0\end{pmatrix}=(C,D)\begin{pmatrix}I_{r}&0\\\
0&0\end{pmatrix}\Leftrightarrow A=C.\end{split}$

So that $G=(C,B)\Sigma^{\rm T}$. Then we can write

$P=\\{(C,B)\Sigma^{\rm
T}:B\in\mathbb{R}^{d_{y}\times(d_{u}-r)},~{}\left\|(C,B)\right\|\leq 2\\},$

By direct computation, for $G=(C,B)\Sigma^{\rm T}\in P$, it holds

(53)
$\begin{split}u^{*}(G)=~{}&\mathop{\arg\min}_{u}\frac{1}{2}\left\|Gu-y\right\|^{2}+\frac{1}{2}\left\|u\right\|^{2}=(G^{\rm
T}G+I)^{-1}G^{\rm T}y\\\ =~{}&\Sigma\left[(C,B)^{\rm
T}(C,B)+I\right]^{-1}(C,B)^{\rm T}y.\end{split}$

Construction of $G_{0}$ and $G_{B}$. Consider the two matrices

$G_{0}=(C,0)\Sigma^{\rm T},\quad G_{B}=(C,B)\Sigma^{\rm T},$

where $B$ is to be determined later. By Eq. 53,

$u^{*}(G_{0})=\Sigma z_{0},\quad z_{0}=\begin{bmatrix}(C^{\rm
T}C+I)^{-1}C^{\rm T}y\\\ 0\end{bmatrix}.$ $u^{*}(G_{B})=\Sigma z_{C},\quad
z_{C}=\begin{pmatrix}C^{\rm T}C+I&C^{\rm T}B\\\ B^{\rm T}C&B^{\rm
T}B+I\end{pmatrix}^{-1}\begin{pmatrix}C^{\rm T}y\\\ B^{\rm T}y\end{pmatrix}.$

Notice

$\|u^{*}(G_{0})-u^{*}(G_{B})\|=\left\|z_{0}-z_{C}\right\|\geq\left\|z_{C,2}\right\|,$

where $z_{C,2}$ is the second block of $z_{C}$, which can be computed directly
by

$z_{C,2}=\left(I+B^{\rm T}\left(I+CC^{\rm T}\right)^{-1}B\right)^{-1}B^{\rm
T}\left(I+CC^{\rm T}\right)^{-1}y.$

Denote $CC^{\rm T}=Q\Lambda Q^{\rm T}$ as the eigenvalue decomposition, and
denote $\tilde{B}=Q^{\rm T}B$, then

$z_{C,2}=\left(I+\tilde{B}^{\rm
T}\left(I+\Lambda^{2}\right)^{-1}\tilde{B}\right)^{-1}\tilde{B}^{\rm
T}\left(I+\Lambda^{2}\right)^{-1}Q^{\rm T}y.$

Denote $r^{\prime}=r\wedge d_{y}$, and note $\Lambda\in\mathbb{R}^{d_{y}\times
d_{y}}$ must be of the form

$\Lambda=\begin{pmatrix}\Lambda_{r^{\prime}}&0\\\
0&0\end{pmatrix},\quad\Lambda_{r^{\prime}}\in\mathbb{R}^{r^{\prime}\times
r^{\prime}}.$

We will choose $\tilde{B}\in\mathbb{R}^{d_{y}\times(d_{u}-r)}$ of the form

$\tilde{B}=\begin{pmatrix}\Lambda^{\prime}&0\\\
E&0\end{pmatrix},\quad\Lambda^{\prime}\in\mathbb{R}^{r^{\prime}\times
r^{\prime}},~{}E\in\mathbb{R}^{(d_{y}-r^{\prime})\times r^{\prime}},$

where $\Lambda^{\prime}$ is diagonal. Note the decomposition is possible since
$d_{u}-r\geq r\geq r^{\prime}$. Also notice when $r\geq d_{y}$,
$\Lambda,\tilde{B}$ only have the upper left block. This does not affect our
deduction, since we will set $E=0$ in this case and the computation is still
correct.

Accordingly decompose $Q^{\rm T}y=(z_{1},z_{2})$ where
$z_{1}\in\mathbb{R}^{r^{\prime}},z_{2}\in\mathbb{R}^{d_{y}-r^{\prime}}$.
Notice $1=\|Q^{\rm
T}y\|^{2}=\left\|z_{1}\right\|^{2}+\left\|z_{2}\right\|^{2}$. By direct
cputation, $z_{C,2}=(z^{\prime},0)$, where

$z^{\prime}=\left(I+\Lambda^{\prime}(I+\Lambda_{r^{\prime}}^{2})^{-1}\Lambda^{\prime}+E^{\rm
T}E\right)^{-1}\left(\Lambda^{\prime}(I+\Lambda_{r^{\prime}}^{2})^{-1}z_{1}+E^{\rm
T}z_{2}\right).$

When $\left\|z_{2}\right\|^{2}\geq 1/2$, we can take
$\Lambda^{\prime}=0,E=\left\|z_{2}\right\|^{-1}z_{2}e_{1}^{\rm T}$ where
$e_{1}=(1,0,\cdots,0)$. (note it is possible since $\left\|z_{2}\right\|>0$
implies $r^{\prime}<d_{y}$) and then

$z^{\prime}=\left(I+E^{\rm T}E\right)^{-1}E^{\rm
T}z_{2}=\frac{1}{2}\left\|z_{2}\right\|e_{1}\Rightarrow\left\|z^{\prime}\right\|=\frac{1}{2}\left\|z_{2}\right\|\geq\frac{1}{2\sqrt{2}}.$

Otherwise, $\left\|z_{1}\right\|^{2}=1-\left\|z_{2}\right\|^{2}\geq 1/2$, we
can take $E=0$, $\Lambda^{\prime}=I$, and then

$z^{\prime}=\left(I+(I+\Lambda_{r^{\prime}}^{2})^{-1}\right)^{-1}(I+\Lambda_{r^{\prime}}^{2})^{-1}z_{1}=(2I+\Lambda_{r^{\prime}}^{2})^{-1}z_{1}.$
$\Rightarrow\left\|z^{\prime}\right\|\geq\frac{\left\|z_{1}\right\|}{2+\left\|\Lambda_{r^{\prime}}\right\|^{2}}\geq\frac{1}{3\sqrt{2}},\quad\text{since
}\left\|\Lambda_{r^{\prime}}\right\|=\|CC^{\rm T}\|\leq\|G^{\prime}\|^{2}\leq
1.$

So that we prove in any case,

$\left\|u^{*}(G_{0})-u^{*}(G_{B})\right\|\geq\left\|z_{C,2}\right\|=\left\|z^{\prime}\right\|\geq\frac{1}{3\sqrt{2}}>\frac{1}{5}.$

The last thing to verify is that $G_{B}$ satisfies the norm constraint
$\left\|G_{B}\right\|\leq 2$. Notice

$\left\|G_{B}\right\|^{2}=\left\|(C,B)\right\|^{2}=\left\|(C,B)(C,B)^{\rm
T}\right\|=\left\|CC^{\rm T}+BB^{\rm
T}\right\|\leq\left\|G^{\prime}\right\|^{2}+\|\tilde{B}\tilde{B}^{\rm T}\|.$

In the above cases,
$\|\tilde{B}\|=\left\|E\right\|=\left\|z_{2}\right\|^{-1}\left\|z_{2}\right\|\left\|e_{1}\right\|=1$
or $\|\tilde{B}\|=\left\|\Lambda^{\prime}\right\|=\left\|I\right\|=1$. So that
$\left\|G_{B}\right\|^{2}\leq 2<4$. This completes our proof.

###  A.4 Stability

######  Lemma A.1.

DEKI with step size $\tilde{h}_{n}=\mu\|C_{n}^{uu}\|^{-1}$ is stable in the
sense that

(54)
$\|z-\mathcal{H}(\overline{u}_{n+1})\|\leq(1+M\sqrt{\mu/2})\|z-\mathcal{H}(\overline{u}_{n})\|.$

######  Proof A.2.

From Eq. 41, we obtain

$\begin{split}\|\overline{u}_{n+1}-\overline{u}_{n}\|\leq~{}&\tilde{h}_{n}\|T_{n}Y_{n}^{\rm
T}(I+\tilde{h}_{n}Y_{n}Y_{n}^{\rm
T})^{-1}(z-\mathcal{H}(\overline{u}_{n}))\|\\\
\leq~{}&\tilde{h}_{n}\|T_{n}Y_{n}^{\rm T}(I+\tilde{h}_{n}Y_{n}Y_{n}^{\rm
T})^{-1}\|\|z-\mathcal{H}(\overline{u}_{n})\|.\end{split}$

Notice $\|T_{n}^{\rm T}T_{n}\|=\left\|T_{n}T_{n}^{\rm
T}\right\|=\|\tilde{C}_{n}^{uu}\|$, and thus

$\begin{split}&\|T_{n}Y_{n}^{\rm T}(I+\tilde{h}_{n}Y_{n}Y_{n}^{\rm
T})^{-1}\|^{2}\\\ =~{}&\|(I+\tilde{h}_{n}Y_{n}Y_{n}^{\rm
T})^{-1}Y_{n}T_{n}^{\rm T}T_{n}Y_{n}^{\rm T}(I+\tilde{h}_{n}Y_{n}Y_{n}^{\rm
T})^{-1}\|\\\ \leq~{}&\|C_{n}^{uu}\|\|(I+\tilde{h}_{n}Y_{n}Y_{n}^{\rm
T})^{-1}Y_{n}Y_{n}^{\rm T}(I+\tilde{h}_{n}Y_{n}Y_{n}^{\rm
T})^{-1}\|.\end{split}$

Using Lemma B.2 (3) with $\psi(x)=x(1+x)^{-2}\leq 1/2$, we obtain

$\tilde{h}_{n}(I+\tilde{h}_{n}Y_{n}Y_{n}^{\rm
T})^{-1}Y_{n}Y_{n}(I+\tilde{h}_{n}Y_{n}Y_{n}^{\rm T})^{-1}\preceq I/2.$

Combining the above estimates, we obtain

$\begin{split}\|\overline{u}_{n+1}-\overline{u}_{n}\|\leq~{}&\sqrt{\frac{\tilde{h}_{n}\|C_{n}^{uu}\|}{2}}\|z-\mathcal{H}(\overline{u}_{n})\|=\sqrt{\frac{\mu}{2}}\|z-\mathcal{H}(\overline{u}_{n})\|.\end{split}$

Since $\mathcal{H}$ is $M$-Lipchitz,

$\begin{split}\|z-\mathcal{H}(\overline{u}_{n+1})\|\leq~{}&\|z-\mathcal{H}(\overline{u}_{n})\|+M\|\overline{u}_{n+1}-\overline{u}_{n}\|\\\
\leq~{}&(1+M\sqrt{\mu/2})\|z-\mathcal{H}(\overline{u}_{n})\|.\end{split}$

##  Appendix B Some matrix analysis results

Here we list some lemmas on the matrix analysis. More details can be found in
textbooks, for instance [12].

######  Lemma B.1.

(Sherman-Morrison-Woodbury) Let $A\in\mathbb{R}^{n\times n}$ be an invertible
matrix. Let $U,V\in\mathbb{R}^{n\times m}$. If $I+V^{\rm T}A^{-1}U$ is
nonsingular, then

(55) $(A+UV^{\rm T})^{-1}=A^{-1}-A^{-1}U(I+V^{\rm T}A^{-1}U)^{-1}V^{\rm
T}A^{-1}.$

######  Lemma B.2.

(Comparison) Let $A,B\in\mathbb{R}^{n\times n}$ be two positive semidefinite
matrices such that $A\preceq B$, then it holds that

  1. (1)

$\lambda_{k}(A)\leq\lambda_{k}(B)$, where $\lambda_{k}(M)$ is the $k$-th
largest eigenvalue of a matrix $M$.

  2. (2)

$\forall C\in\mathbb{R}^{m\times n},CAC^{\rm T}\preceq CBC^{\rm T}$.

  3. (3)

$\forall\psi,\varphi~{}\text{s.t.
}\psi(x)\leq\varphi(x),~{}\psi(A)\preceq\varphi(A)$. As a corollary, when
$\lambda_{k}(A)\in[a,b]$ for all $k$ and $\psi(x)\leq\varphi(x)$ when
$x\in[a,b]$, then $\psi(A)\preceq\varphi(A)$.

  4. (4)

$\forall\psi:\mathbb{R}\rightarrow\mathbb{R}$ that is monotone increasing,
$\psi(A)\preceq\psi(B)$. Similarly, for monotone decreasing function $\psi$,
$\psi(A)\succeq\psi(B)$.

######  Proof B.3.

For (1), from Courant minimax principle,

$\lambda_{k}(A)=\max_{V\subset\mathbb{R}^{n},\dim V=k}\min_{0\neq u\in
V}\frac{u^{\rm T}Au}{u^{\rm T}u}.$

Since $A\preceq B\Rightarrow u^{\rm T}Au\leq u^{\rm T}Bu$, we conclude that
$\lambda_{k}(A)\leq\lambda_{k}(B)$. For (2), notice

$\forall u\in\mathbb{R}^{m},~{}u^{\rm T}CAC^{\rm T}u=(C^{\rm T}u)^{\rm
T}A(C^{\rm T}u)\leq(C^{\rm T}u)^{\rm T}B(C^{\rm T}u)=u^{\rm T}CBC^{\rm T}u.$

For (3), denote the eigen decomposition of $A$ is $A=Q\Lambda Q^{\rm T}$, then
by definition $\psi(A)=Q\psi(\Lambda)Q^{\rm T}$ and
$\varphi(A)=Q\varphi(\Lambda)Q^{\rm T}$. Notice

$\psi(\Lambda)=\text{diag}\\{\psi(\lambda_{1}),\dots,\psi(\lambda_{n})\\}\preceq\text{diag}\\{\varphi(\lambda_{1}),\dots,\varphi(\lambda_{n})\\}=\varphi(\Lambda),$

since $\psi(\lambda_{k})\leq\varphi(\lambda_{k})$. By (2), we conclude that
$\psi(A)\preceq\varphi(A)$. For the corollary, just notice we can modify
$\psi(x)=\varphi(x)$ outside $[a,b]$ and $\psi(A),\varphi(A)$ are unchanged.
For (4), notice from (1), $\lambda_{k}(A)\leq\lambda_{k}(B)$, the conclusion
follows by (2) and the monotonicity of $\psi$ using the same claim.

######  Lemma B.4.

Let $A,B\in\mathbb{R}^{n\times n}$ be two positive semidefinite matrices. If
there exists $0<m\leq M,\text{s.t. }mI\preceq A\preceq MI$, then it holds

(56) $B(I+MB)^{-2}\preceq(I+BA)^{-1}B(I+AB)^{-1}\preceq B(I+mB)^{-2}.$

######  Proof B.5.

We first prove the inequality for nonsingular $B$. Notice

$\begin{split}&(I+BA)^{-1}B(I+AB)^{-1}\\\
=~{}&(B^{1/2}(I+B^{1/2}AB^{1/2})B^{-1/2})^{-1}B(B^{-1/2}(I+B^{1/2}AB^{1/2})B^{1/2})^{-1}\\\
=~{}&B^{1/2}(I+B^{1/2}AB^{1/2})^{-2}B^{1/2}.\end{split}$

Using Lemma B.2 (2), since $mI\preceq A\preceq MI$, we obtain

$mB\preceq B^{1/2}AB^{1/2}\preceq MB,$

As a result, using Lemma B.2 (4) for $\psi(x)=(1+x)^{-2}$,

$\begin{split}B^{1/2}(I+MB)^{-2}B^{1/2}\preceq~{}&B^{1/2}(I+B^{1/2}AB^{1/2})^{-2}B^{1/2}\\\
\preceq~{}&B^{1/2}(I+mB)^{-2}B^{1/2}.\end{split}$

The conclusion follows by observing that the bounds are of functions of $B$
are thus the sequence of multiplication is unimportant.

If $B$ is nonsingular, we can use the same argument to show

$B_{\epsilon}(I+MB_{\epsilon})^{-2}\preceq(I+B_{\epsilon}A)^{-1}B_{\epsilon}(I+AB_{\epsilon})^{-1}\preceq
B_{\epsilon}(I+mB_{\epsilon})^{-2}.$

with $B_{\epsilon}=B+\varepsilon I$ and take limit $\varepsilon\rightarrow 0$.

##  Appendix C Implementation of the linearization

Denote the deviation ensemble of $u_{n}$ and $\mathcal{G}u_{n}$ as

$T_{n}=[\tau_{n}^{(1)},\dots,\tau_{n}^{(J)}]\in\mathbb{R}^{d_{u}\times
J},\quad
Y_{n}=[\mathfrak{r}_{n}^{(1)},\dots,\mathfrak{r}_{n}^{(J)}]\in\mathbb{R}^{d_{y}\times
J},$

where
$\tau_{n}^{(j)}=u_{n}^{(j)}-\overline{u}_{n},\mathfrak{r}_{n}^{(j)}=\mathcal{G}u_{n}^{(j)}-\overline{\mathcal{G}u_{n}}$.
Denote the rank of $T_{n}$ and $Y_{n}$ as $r,r^{\prime}$ respectively. Perform
the reduced QR decomposition

$T_{n}=V_{n}Q_{n},~{}V_{n}\in\mathbb{R}^{d_{u}\times r},\quad
Y_{n}=W_{n}R_{n},~{}W_{n}\in\mathbb{R}^{d_{y}\times r^{\prime}}.$

Denote $\Lambda_{n}=Q_{n}Q_{n}^{\rm T}$, and then

$(J-1)C_{n}^{uu}=T_{n}T_{n}^{\rm T}=V_{n}Q_{n}Q_{n}^{\rm T}V_{n}^{\rm
T}=V_{n}\Lambda_{n}V_{n}^{\rm T}.$ $(J-1)C_{n}^{pu}=Y_{n}T_{n}^{\rm
T}=W_{n}R_{n}Q_{n}^{\rm T}V_{n}^{\rm T}.$

Therefore,

$C_{n}^{pu}(C_{n}^{uu})^{\dagger}=W_{n}R_{n}Q_{n}^{\rm
T}\Lambda_{n}^{-1}V_{n}^{\rm T}.$

Introduce

$\hat{R}_{n}=\mathcal{T}_{M}(R_{n}Q_{n}^{\rm
T}\Lambda_{n}^{-1})\Lambda_{n}Q_{n},\quad\hat{Y}_{n}=W_{n}\hat{R}_{n},$

where $\mathcal{T}_{M}$ is the truncation operator on the singular value:

$\mathcal{T}_{M}A=U\text{diag}\\{\Lambda_{ii}\wedge M\\}V^{\rm
T},\quad\text{if }A=U\Lambda V^{\rm T}.$

Then it suffices to define $\hat{\mathcal{G}}u_{n}^{(j)}$ by

$\hat{\mathcal{G}}u_{n}^{(j)}=\overline{\mathcal{G}u_{n}}+\hat{\mathfrak{r}}_{n}^{(j)},\quad\hat{\mathfrak{r}}_{n}^{(j)}=W_{n}\hat{R}_{n}(j),$

where $\hat{R}_{n}(j)$ is the $j$-th column of $\hat{R}_{n}$.

##  Appendix D LEKI and localization

LEKI method is proposed in [30]. We use the following discrete-time scheme:

$u_{n+1}^{(j)}=u_{n}^{(j)}+h_{n}\tilde{C}_{n}^{uz}\left(\mathcal{H}(u_{n}^{(j)})-z\right)+\sigma^{2}\xi_{n}^{(j)}.$

The localized covariance is defined via a localization matrix
$\Psi\in\mathbb{R}^{d_{z}\times d_{u}}$ and

$\tilde{C}_{n}^{uz}=\Psi\circ C_{n}^{uz}.$

And $\xi_{n}^{(j)}$ is some artificial noise chosen following [30] so that

$\frac{1}{J-1}\sum_{j=1}^{J}\left(\xi_{n}^{(j)}\otimes\tau_{n}^{(j)}+\tau_{n}^{(j)}\otimes\xi_{n}^{(j)}\right)=\Sigma_{n},$

where the diagonal entries of $\Sigma_{n}$ are all $1$. The noise level is
fixed $\sigma=10^{-3}$.

For the linear transport model, we compare two different localization designs.
Both uses the Gaspari Cohn function $\psi$ [14], and

$\Psi(i,j)=\psi(|i-i_{o}(j)|/r_{loc}),$

where $i_{o}(j)$ is the coordinate of the parameter that mostly influences the
$j$-th observation, and $r_{loc}$ is the localization radius, see [30]. The
one with knowledge of the speed uses $r_{loc}=1.5$ and

$i_{o}(j)=\lfloor d_{y}(i/d_{u}+a)\rfloor.$

And the one without the knowledge simply chooses $i_{o}(j)=j$ with the same
localization radius.