# Adjoint-Free 4D-Var Methods Via Line Search Optimization For Non-Linear Data
Assimilation

Elias D. Nino-Ruiz Electronic address<EMAIL_ADDRESS>Corresponding
author Applied Math and Computer Science Laboratory
Department of Computer Science
Universidad del Norte, Colombia  J. Diaz-Rodriguez Electronic address:
<EMAIL_ADDRESS>Department of Mathematics and Statistics
York University, Canada

###### Abstract

This paper proposes two practical implementations of Four-Dimensional
Variational (4D-Var) Ensemble Kalman Filter (4D-EnKF) methods for non-linear
data assimilation. Our formulations’ main idea is to avoid the intrinsic need
for adjoint models in the context of 4D-Var optimization and, even more, to
handle non-linear observation operators during the assimilation of
observations. The proposed methods work as follows: snapshots of an ensemble
of model realizations are taken at observation times, these snapshots are
employed to build control spaces onto which analysis increments can be
estimated. Via the linearization of observation operators at observation
times, a line-search based optimization method is proposed to estimate optimal
analysis increments. The convergence of this method is theoretically proven as
long as the dimension of control-spaces equals model one. In the first
formulation, control spaces are given by full-rank square root approximations
of background error covariance matrices via the Bickel and Levina precision
matrix estimator. In this context, we propose an iterative Woodbury matrix
formula to perform the optimization steps efficiently. The last formulation
can be considered as an extension of the Maximum Likelihood Ensemble Filter to
the 4D-Var context. This employs pseudo-square root approximations of prior
error covariance matrices to build control spaces. Experimental tests are
performed by using the Lorenz 96 model. The results reveal that, in terms of
Root-Mean-Square-Error values, both methods can obtain reasonable estimates of
posterior error modes in the 4D-Var optimization problem. Moreover, the
accuracies of the proposed filter implementations can be improved as the
ensemble sizes are increased.

##  1 Introduction

To be concise, the goal of Data Assimilation (DA) is to correct imperfect
model trajectories $\left\\{{\bf x}^{b}_{k}\right\\}_{k=0}^{M}$ via real-noisy
observations $\left\\{{\bf y}_{k}\right\\}_{k=0}^{M}$ [21, 24, 41], where
${\bf x}^{b}_{k}\in\mathbb{R}^{n\times 1}$ is well-known as the forecast
state, $k$ stands for time index where observations are collected, for $0\leq
k\leq M$, $n$ denotes the state size, the size of the assimilation window
reads $M$, ${\bf y}_{k}\in\mathbb{R}^{m\times 1}$ is a vector of observations,
and $m$ is the number of observations. Forecasts are typically obtained by
numerical models which minimizes the behaviour of the ocean and/or the
atmosphere:

$\displaystyle{\bf x}_{k}=\mathcal{M}_{t_{k-1}\rightarrow t_{k}}\left({\bf
x}_{k-1}\right)\,\,\text{ for $0\leq k-1\leq M$, and ${\bf
x}_{k}\in\mathbb{R}^{n\times 1}$}\,,$ (1)

where $\mathcal{M}:\mathbb{R}^{n\times 1}\rightarrow\mathbb{R}^{n\times 1}$ is
the non-linear model operator (typically, given in the form of partial and
ordinary differential equations). In the context of Four-Dimensional
Variational (4D-Var) data assimilation methods [20],cost functions of the form
[48, 49]:

$\displaystyle\displaystyle\mathcal{J}({\bf
x}_{0})=\frac{1}{2}\cdot\left\|{\bf x}_{0}-{\bf x}_{0}^{b}\right\|^{2}_{{\bf
B}^{-1}_{0}}+\frac{1}{2}\cdot\sum_{k=0}^{M-1}\left\|{\bf
y}_{k}-\mathcal{H}_{k}\left({\bf x}_{k}\right)\right\|^{2}_{{\bf
R}_{k}^{-1}}\,,$ (2)

are considered to estimate the initial state ${\bf
x}_{0}\in\mathbb{R}^{n\times 1}$ which best fit a given assimilation window
$\left\\{{\bf y}_{k}\right\\}_{k=0}^{M}$. The optimization problem to solve
reads:

$\displaystyle\displaystyle{\bf x}^{a}_{0}=\arg\underset{{\bf
x}_{0}}{\min}\,\mathcal{J}\left({\bf x}_{0}\right)\,,\text{subject to ${\bf
x}_{k}=\mathcal{M}_{t_{k-1}\rightarrow t_{k}}\left({\bf x}_{k-1}\right)$}\,.$
(3)

A potential drawback of 4D-Var based methods is the use of adjoint models. In
practice, these models can be labor-intensive to develop and computationally
expensive to run [43]. Thus, ensemble-based methods can be employed under
4D-Var operational scenarios to avoid the use of adjoint models [17, 32]. The
main idea behind these approaches is to build ensemble sub-spaces onto which
analysis increments can be estimated [53]. For non-linear observation
operators, the underlying error distribution of background errors can be non-
Gaussian and as a direct consequence, no closed-form expression are available
for computing posterior modes in (3) [1, 48]. Hence, iterative methods such as
Trust Region [33] and Line Search [34] can be employed to estimate posterior
modes of the implicit error distribution in (2). We think that there is an
opportunity to derive a line-search-based optimization method for the solution
of (3) via control-space approximations (i.e., based on ensembles). In this
manner, an ensemble of model realizations can be employed to capture the error
dynamics across all observation times and, even more, to build spaces onto
which analysis increments can be estimated. Simple ensemble anomalies can span
Control-spaces or robust square-root approximations of background-error-
covariance matrices via a modified Cholesky decomposition [3, 42]. In the last
estimator, its special structure can allow matrix-free 4D-EnKF
implementations.

This paper’s outline is as follows: in Section 2, theoretical and practical
aspects of variational and ensemble data assimilation methods are briefly
discussed in the linear and the non-linear cases. In Section 4, two efficient
implementations of 4D-EnKF formulations are proposed via line-search
optimization, the first one makes use of a modified Cholesky decomposition to
build control spaces while the last one uses an ensemble of anomalies to
estimate posterior increments. Their computational costs are detailed as well
as the convergence of the line-search based optimization method. In Section 5,
experimental tests are performed to assess the accuracy of the proposed filter
implementations, the Lorenz 96 model is employed as our surrogate model for
the experiments. Conclusions of this research are stated in Section 6.

##  2 Preliminaries

This Section discusses some related topics for the development of our
variational filter.

##  3 The Ensemble Kalman Filter

The Ensemble Kalman Filter (EnKF) [29, 9] is a sequential ensemble based
method which employes an ensemble of model realizations to estimate moments of
prior error distributions [19, 47]:

$\displaystyle{\bf X}_{k}^{b}=\left[{\bf x}^{b[1]}_{k},\,{\bf
x}^{b[2]}_{k},\,\ldots,\,{\bf x}^{b[N]}_{k}\right]\in\mathbb{R}^{n\times
N}\,,$ (4)

where ${\bf x}^{b[e]}_{k}\in\mathbb{R}^{n\times 1}$ denotes the $e$-th
ensemble member, for $1\leq e\leq N$, at time $k$, for $0\leq k\leq M$. The
empirical moments of (4) are employed to estimate the forecast state ${\bf
x}^{b}_{k}$:

$\displaystyle\displaystyle\overline{\bf
x}^{b}_{k}=\frac{1}{N}\cdot\sum_{e=1}^{N}{\bf
x}^{b[e]}_{k}\in\mathbb{R}^{n\times 1}\,,$ (5)

and the backgrund error covariance matrix ${\bf B}$:

$\displaystyle\displaystyle{\bf
P}^{b}_{k}=\frac{1}{N-1}\cdot\boldsymbol{\Delta}{\bf
X}^{b}_{k}\cdot\left[\boldsymbol{\Delta}{\bf
X}^{b}_{k}\right]^{T}\in\mathbb{R}^{n\times n}\,,$ (6)

where the matrix of member deviations reads:

$\displaystyle\displaystyle\boldsymbol{\Delta}{\bf X}^{b}_{k}={\bf
X}^{b}_{k}-\overline{\bf x}^{b}_{k}\cdot{\bf 1}^{T}\in\mathbb{R}^{n\times
N}\,.$ (7)

The posterior ensemble can then be built by using synthetic observations:

$\displaystyle{\bf X}^{a}_{k}={\bf X}^{b}_{k}+\boldsymbol{\Delta}{\bf
X}^{a}_{k}\,,$ (8)

where the analysis updates can be obtained by solving the following linear
system:

$\displaystyle\left[\left[{\bf P}^{b}_{k}\right]^{-1}+{\bf H}_{k}^{T}\cdot{\bf
R}_{k}^{-1}\cdot{\bf H}_{k}\right]\cdot\boldsymbol{\Delta}{\bf X}^{a}_{k}={\bf
H}_{k}^{T}\cdot{\bf R}_{k}^{-1}\cdot{\bf D}^{s}_{k}\in\mathbb{R}^{n\times
N}\,,$ (9)

the innovation matrix reads ${\bf D}_{k}^{s}\in\mathbb{R}^{m\times N}$ whose
$e$-th column ${\bf y}_{k}-{\bf H}_{k}\cdot{\bf
x}^{b[e]}_{k}+\boldsymbol{\varepsilon}_{k}^{[e]}\in\mathbb{R}^{m\times 1}$ is
a synthetic observation with
$\boldsymbol{\varepsilon}_{k}^{[e]}\sim\mathcal{N}\left({\bf 0}_{m},\,{\bf
R}_{k}\right)$. In operational DA scenarios, the ensemble size can be lesser
than model dimensions by order of magnitudes and as a direct consequence,
sampling errors impact the quality of analysis increments. To counteract the
effects of sampling noise, localizations methods are commonly employed [13,
5]. In the context of covariance tappering, the use of the spatial-
predecessors concept can be employed to obtain sparse estimators of precision
matrices [25]. The predecessors of model component $i$, from now on
$P(i,\,r)$, for $1\leq i\leq n$ and a radius of influence
$r\in\mathbb{Z}^{+}$, are given by the set of components whose labels are
lesser than that of the $i$-th one.

In the EnKF based on a modified Cholesky decomposition (EnKF-MC) [37, 38] the
following estimator is employed to approximate the precision covariance matrix
of the background error distribution [4]:

$\displaystyle\widehat{\bf B}^{-1}_{k}=\widehat{\bf
L}_{k}^{T}\cdot\widehat{\bf D}^{-1}_{k}\cdot\widehat{\bf
L}_{k}\in\mathbb{R}^{n\times n}\,,$ (10)

where the Cholesky factor ${\bf L}_{k}\in\mathbb{R}^{n\times n}$ is a lower
triangular matrix,

$\displaystyle\left\\{\widehat{\bf
L}_{k}\right\\}_{i,v}=\begin{cases}-\beta_{i,v,k}&\,,\,v\in P(i,r)\\\
1&\,,\,i=v\\\ 0&\,,\,otherwise\end{cases}\,,$ (11)

whose non-zero sub-diagonal elements $\beta_{i,v,k}$ are obtained by fitting
models of the form,

$\displaystyle\displaystyle{{\bf x}_{[i]}^{T}}_{k}=\sum_{v\in
P(i,\,r)}\beta_{i,v,k}\cdot{{\bf
x}_{[v]}^{T}}_{k}+{\boldsymbol{\gamma}_{i}}_{k}\in\mathbb{R}^{N\times
1}\,,\text{ $1\leq i\leq n$}\,,$ (12)

where ${{\bf x}_{[i]}^{T}}_{k}\in\mathbb{R}^{N\times 1}$ denotes the $i$-th
row (model component) of the ensemble (4), components of vector
${\boldsymbol{\gamma}_{i}}_{k}\in\mathbb{R}^{N\times 1}$ are samples from a
zero-mean Normal distribution with unknown variance $\sigma^{2}_{k}$, and
${\bf D}_{k}\in\mathbb{R}^{n\times n}$ is a diagonal matrix whose diagonal
elements read,

$\displaystyle\left\\{{\bf D}_{k}\right\\}_{i,i}$ $\displaystyle=$
$\displaystyle\widehat{\bf var}\left({{\bf x}_{[i]}^{T}}_{k}-\sum_{v\in
P(i,\,r)}\beta_{i,v,k}\cdot{{\bf x}_{[j]}^{T}}_{k}\right)^{-1}$ (13)
$\displaystyle\approx$ $\displaystyle{\bf
var}\left({\boldsymbol{\gamma}_{i}}_{k}\right)^{-1}=\frac{1}{\sigma_{k}^{2}}>0\,,\text{
with $\left\\{{\bf D}_{k}\right\\}_{1,1}=\widehat{\bf var}\left({{\bf
x}_{[1]}^{T}}_{k}\right)^{-1}$}\,,$ (14)

where ${\bf var}(\bullet)$ and $\widehat{\bf var}(\bullet)$ denote the actual
and the empirical variances, respectively. The analysis equations can then be
written as discussed in [38].

###  3.1 Maximum Likelihood Ensemble Filter

The filter derivation in the EnKF context is based on the assumptions that
prior and observational errors are Gaussian. However, for non-linear
observation operators (i.e., information brought by satellite radiances),
observational errors are not Gaussian anymore and as a direct consequence,
EnKF based methods can fail to obtain reasonable estimates of posterior error
distributions. The Maximum Likelihood Ensemble Filter (MLEF) [56] is a
numerical method that minimizes the Three-Dimensional Variational (3D-Var)
cost function:

$\displaystyle\displaystyle\mathcal{J}({\bf
x}_{k})=\frac{1}{2}\cdot\left\|{\bf x}_{k}-\overline{\bf
x}_{k}^{b}\right\|^{2}_{{\bf B}^{-1}_{k}}+\frac{1}{2}\cdot\left\|{\bf
y}_{k}-\mathcal{H}_{k}\left({\bf x}_{k}\right)\right\|^{2}_{{\bf
R}_{k}^{-1}}\,,$ (15)

onto the ensemble space. For this purpose, analysis increments are constrained
to the space spanned by the ensemble of anomalies (7), this is:

$\displaystyle\displaystyle{\bf x}_{k}=\overline{\bf
x}^{b}_{k}+\boldsymbol{\Delta}{\bf X}_{k}\cdot{\bf w}\,,$ (16)

which results in a cost function of the form:

$\displaystyle\displaystyle\mathcal{J}(\overline{\bf
x}^{b}_{k}+\boldsymbol{\Delta}{\bf X}_{k}\cdot{\bf w})$ $\displaystyle=$
$\displaystyle\frac{N-1}{2}\cdot\left\|{\bf w}\right\|^{2}$ (17)
$\displaystyle+$ $\displaystyle\frac{1}{2}\cdot\left\|{\bf
y}_{k}-\mathcal{H}_{k}\left(\overline{\bf x}^{b}_{k}+\boldsymbol{\Delta}{\bf
X}_{k}\cdot{\bf w}\right)\right\|^{2}_{{\bf R}_{k}^{-1}}\,,$

and corresponding optimization problem:

$\displaystyle\displaystyle{\bf w}^{*}=\arg\underset{{\bf
w}}{\min}\,\mathcal{J}(\overline{\bf x}^{b}_{k}+\boldsymbol{\Delta}{\bf
X}_{k}\cdot{\bf w})\,.$ (18)

Note that, since all computations are performed onto the ensemble space, the
computational effort of the MLEF becomes:

$\displaystyle\displaystyle\mathcal{O}\left(N^{3}+n\cdot m\cdot N\right)\,,$
(19)

once all analysis increments are mapped back onto the model space (and the
posterior ensemble is built about the estimated posterior mode), which is
linear with regard to the model dimension.

###  3.2 Four Dimensional Variational Ensemble Kalman Filter

In Four Dimensional Variational Ensemble Kalman Filter (4D-EnKF) methods [27,
30], the model trajectory in (3) is constrained to the space spanned by the
background ensemble members, this is:

$\displaystyle\displaystyle{\bf x}_{k}=\overline{\bf
x}^{b}_{k}+\boldsymbol{\Delta}{\bf X}_{k}\cdot{\bf w}\,,$ (20)

where ${\bf w}\in\mathbb{R}^{N\times 1}$ is a vector in redundant coordinates
to be determined later. Making the assumption that all sources of errors are
Gaussian (i.e., the observation operator is linear), and by replacing (20)
into the equation (2) one obtains:

$\displaystyle\mathcal{J}({\bf x}_{0})$ $\displaystyle=$
$\displaystyle\mathcal{J}\left(\overline{\bf
x}^{b}_{0}+\boldsymbol{\Delta}{\bf X}_{0}\cdot{\bf w}\right)$ (21)
$\displaystyle=$ $\displaystyle\widehat{\mathcal{J}}({\bf
w})=\frac{(N-1)}{2}\cdot\left\|{\bf
w}\right\|^{2}+\frac{1}{2}\cdot\sum_{k=0}^{M}\left\|{\bf d}_{k}-{\bf
Q}_{k}\cdot{\bf w}\right\|^{2}_{{\bf R}_{k}^{-1}}\,,$

where ${\bf d}_{k}={\bf y}_{k}-{\bf H}_{k}\cdot\overline{\bf
x}^{b}_{k}\in\mathbb{R}^{m\times 1}$ is the innovation vector and ${\bf
Q}_{k}={\bf H}_{k}\cdot\boldsymbol{\Delta}{\bf
X}_{k}^{b}\in\mathbb{R}^{m\times N}$. The optimal values of the control
variable ${\bf w}$ is then sought in order to estimate the initial analysis
state:

$\displaystyle\displaystyle{\bf w}^{*}=\arg\,\underset{{\bf
w}}{\min}\,\widehat{\mathcal{J}}({\bf w})\,.$ (22)

The gradient of (21) equals:

$\displaystyle\displaystyle\nabla_{{\bf w}}\widehat{\mathcal{J}}({\bf w})$
$\displaystyle=$ $\displaystyle(N-1)\cdot{\bf w}-\sum_{k=0}^{M}{\bf
Q}_{k}^{T}\cdot{\bf R}_{k}^{-1}\cdot\left[{\bf d}_{k}-{\bf Q}_{k}\cdot{\bf
w}\right]$ (23) $\displaystyle=$ $\displaystyle\left[(N-1)\cdot{\bf
I}+\sum_{k=0}^{M}{\bf Q}_{k}^{T}\cdot{\bf R}_{k}^{-1}\cdot{\bf
Q}_{k}\right]\cdot{\bf w}$ $\displaystyle-$ $\displaystyle\sum_{k=0}^{M}{\bf
Q}_{k}^{T}\cdot{\bf R}_{k}^{-1}\cdot{\bf d}_{k}\in\mathbb{R}^{N\times 1}\,,$

and by setting this gradient to zero, the optimal weights (22) reads:

$\displaystyle\displaystyle{\bf w}^{*}=\left[(N-1)\cdot{\bf
I}+\sum_{k=0}^{M}{\bf Q}_{k}^{T}\cdot{\bf R}_{k}^{-1}\cdot{\bf
Q}_{k}\right]^{-1}\cdot\sum_{k=0}^{M}{\bf Q}_{k}^{T}\cdot{\bf
R}_{k}^{-1}\cdot{\bf d}_{k}\in\mathbb{R}^{N\times 1}\,,$ (24)

from which the initial analysis state can be estimated:

$\displaystyle\displaystyle\overline{\bf x}_{0}^{a}=\overline{\bf
x}_{0}^{b}+\boldsymbol{\Delta}{\bf X}_{0}^{b}\cdot{\bf w}^{*}\,.$ (25)

Similar to the MLEF method, the computational effort of 4D-EnKF based methods
are linearly bounded by the model dimension; it can be easily seen that all
computations are performed onto the ensemble space. Note that, in this
formulation, posterior estimates can be highly biased regarding the actual
state when the observation operator is non-linear.

###  3.3 Numerical Optimization

In numerical optimization [52, 54], commonly, problems of the form:

$\displaystyle{\bf x}=\arg\,\underset{{\bf x}}{\min}\,f({\bf x})\,,$

where $f({\bf x}):\mathbb{R}^{n\times 1}\rightarrow\mathbb{R}$ is a smooth
function, are solved by iterative schemes such as

$\displaystyle{\bf x}^{(u+1)}={\bf x}^{(u)}+\boldsymbol{\Delta}{\bf
x}^{(u)}\,,$ (26)

wherein $\boldsymbol{\Delta}{\bf x}^{(u)}$ is a search direction, for
instance, the steepest descent direction [44, 16, 11, 26] the Newton’s step
[2, 15, 40], or quasi-Newton based method [45, 39, 28]. A concise survey of
Newton based methods can be consulted in [23]. Since step sizes in steepest
descent directions can be too large, line search [14, 50, 18] or trust region
[6, 33, 7] methods can be successfully used in order to optimally choose the
step size and therefore, to ensure global convergence of iterates to
stationary points. This holds as long as some assumptions over functions,
gradients, and (potentially) Hessians are preserved [46]. In the context of
line search, the following assumptions are commonly done:

  1. 1.

The function $f({\bf x})$ has lower bound on $\Omega_{0}=\\{{\bf
x}\in\mathbb{R}^{n\times 1},\,f({\bf x})\leq f({\bf x}_{0})\\}$, where ${\bf
x}_{0}\in\mathbb{R}^{n\times 1}$ is available.

  2. 2.

The gradient $\nabla f({\bf x})$ is Lipschitz continuous in an open convex set
$B$ which contains $\Omega_{0}$,

$\displaystyle\left\|\nabla f({\bf x})-\nabla f({\bf z})\right\|\leq
L\cdot\left\|{\bf x}-{\bf z}\right\|,\,\text{ for ${\bf x},\,{\bf z}\in B$,
and $L>0$}.$

These together with iterates of the form,

$\displaystyle\displaystyle{\bf x}^{(u+1)}={\bf
x}^{(u)}+\alpha\cdot\boldsymbol{\Delta}{\bf x}^{(u)}\,,$ (27)

ensure global convergence [55] as long as $\alpha$ is chosen as an
(approximated) minimizer of

$\displaystyle\displaystyle\alpha^{*}=\arg\,\underset{\alpha\geq
0}{\min}\,f\left({\bf x}^{(u)}+\alpha\cdot\boldsymbol{\Delta}{\bf
x}^{(u)}\right)\,.$ (28)

##  4 Proposed Method

In this section, we develop an efficient and practical implementation of a
4D-EnKF method via a modified Cholesky decomposition for non-linear
observation operator.

###  4.1 Filter Derivation

Our main goal is to find an initial analysis ensemble ${\bf
X}^{a}_{0}\in\mathbb{R}^{n\times N}$ whose model trajectories account for the
information brought by observations within a given assimilation window. To
accomplish this, we proceed as follows: given an assimilation window of $M$
observations (with no loss of generality, evenly distributed in time), by
using the numerical model (1), the initial background ensemble ${\bf
X}_{0}^{b}$ is propagated in time and snapshots of its time evolution are
taken at observation times. Denote by $k$, the $k$-th step where an
observation is available, for $0\leq k\leq M$, at each observation time, the
background ensemble (4) can be employed to estimate a full-rank square-root
approximation of the precision covariance matrix ${\bf B}_{k}^{-1}$ via a
modified Cholesky decomposition (10):

$\displaystyle\displaystyle\widehat{\bf B}^{-1/2}_{k}=\widehat{\bf
L}^{T}_{k}\cdot\widehat{\bf D}^{-1/2}_{k}\in\mathbb{R}^{n\times n}\,.$ (29)

This square-root estimation can then be employed to build a control space onto
which analysis increments can be estimated:

$\displaystyle{\bf x}_{k}-\overline{\bf x}^{b}_{k}\in{\bf
range}\left\\{\widehat{\bf B}^{1/2}_{k}\right\\}\,,$

or equivalently:

$\displaystyle\displaystyle{\bf x}_{k}=\overline{\bf x}^{b}_{k}+\widehat{\bf
B}^{1/2}_{k}\cdot\boldsymbol{\alpha}\in\mathbb{R}^{n\times 1}\,,$ (30)

where $\boldsymbol{\alpha}\in\mathbb{R}^{n\times 1}$ is a vector in redundant
coordinates to be determined later. Since the background error covariance
matrix onto the control space (30) is nothing but the identity matrix, the
following error statistics hold for the prior weights
$\boldsymbol{\alpha}^{b}$:

$\displaystyle\displaystyle\boldsymbol{\alpha}^{b}\sim\mathcal{N}\left({\bf
0},\,{\bf I}\right)\,.$

Due to this, the 4D-Var cost function (2) onto the space (30) can be written
as follows:

$\displaystyle\mathcal{J}\left({\bf x}_{0}\right)$ $\displaystyle=$
$\displaystyle\mathcal{J}\left(\overline{\bf x}^{b}_{0}+\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}\right)=\frac{1}{2}\cdot\left\|\boldsymbol{\alpha}\right\|^{2}$
(31) $\displaystyle+$ $\displaystyle\frac{1}{2}\cdot\sum_{k=0}^{M}\left\|{\bf
y}_{k}-\mathcal{H}_{k}\left(\overline{\bf x}^{b}_{k}+\widehat{\bf
B}^{1/2}_{k}\cdot\boldsymbol{\alpha}\right)\right\|_{{\bf R}_{k}^{-1}}^{2}\,,$

The adjoint-free optimization problem to solve reads:

$\displaystyle\displaystyle\boldsymbol{\alpha}^{a}=\arg\,\underset{\boldsymbol{\alpha}}{\min}\,\mathcal{J}\left(\overline{\bf
x}^{b}_{0}+\widehat{\bf B}^{1/2}_{0}\cdot\boldsymbol{\alpha}\right)\,.$ (32)

The next section develops a Line-Search method for solving the optimization
problem (32).

###  4.2 An Adjoint-Free 4D-Var Method Via Modified Cholesky Decomposition
(4D-Var-MC)

To approximate a solution of (32), we propose iterates of the form:

$\displaystyle\displaystyle{\bf x}^{(u+1)}_{k}={\bf x}^{(u)}_{k}+\widehat{\bf
B}^{1/2}_{k}\cdot\boldsymbol{\alpha}_{u}\,,\text{ for $0\leq u\leq U$}\,,$
(33)

where $u$ denotes iteration step, $U$ stands for the number of iterations, and
$\boldsymbol{\alpha}_{u}\in\mathbb{R}^{n\times 1}$ is a vector in redundant
coordinates to be determinated later. The background trajectory is set as the
initial solution ${\bf x}^{(0)}_{k}=\overline{\bf x}^{b}_{k}$, for $0\leq
k\leq M-1$, and therefore $\boldsymbol{\alpha}_{0}={\bf 0}$. Note that, at any
iteration $u$,

$\displaystyle\displaystyle{\bf x}_{k}^{(u+1)}$ $\displaystyle=$
$\displaystyle{\bf x}^{(0)}_{k}+\sum_{v=0}^{u}\widehat{\bf
B}^{1/2}_{k}\cdot\boldsymbol{\alpha}_{v}=\overline{\bf
x}^{b}_{k}+\sum_{v=0}^{u-1}\widehat{\bf
B}^{1/2}_{k}\cdot\boldsymbol{\alpha}_{v}+\widehat{\bf
B}^{1/2}_{k}\cdot\boldsymbol{\alpha}_{u}\,$ (34) $\displaystyle=$
$\displaystyle\overline{\bf x}^{b}_{k}+\widehat{\bf
B}^{1/2}_{k}\cdot\left[\sum_{v=0}^{u-1}\boldsymbol{\alpha}_{v}+\boldsymbol{\alpha}_{u}\right]=\overline{\bf
x}^{b}_{k}+\widehat{\bf
B}^{1/2}_{k}\cdot\left[\boldsymbol{\beta}_{u}+\boldsymbol{\alpha}_{u}\right]\,$

where we have let
$\boldsymbol{\beta}_{u}=\sum_{v=0}^{u-1}\boldsymbol{\alpha}_{v}$. By replacing
(33) and (34) in (2) one obtains:

$\displaystyle\displaystyle\widehat{\mathcal{J}}\left(\boldsymbol{\alpha}_{u}\right)=\mathcal{J}\left({\bf
x}^{(u)}_{k}+\widehat{\bf B}^{1/2}_{k}\cdot\boldsymbol{\alpha}_{u}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\cdot\left\|\boldsymbol{\beta}_{u}+\boldsymbol{\alpha}_{u}\right\|^{2}$
$\displaystyle+$ $\displaystyle\frac{1}{2}\cdot\sum_{k=0}^{M-1}\left\|{\bf
y}_{k}-\mathcal{H}_{k}\left({\bf x}^{(u)}_{k}+\widehat{\bf
B}^{1/2}_{k}\cdot\boldsymbol{\alpha}_{u}\right)\right\|^{2}_{{\bf
R}_{k}^{-1}}\,.$

At iteration $u$, we choose $\boldsymbol{\alpha}_{u}$ in (33) as the value
that solves the following optimization problem:

$\displaystyle\displaystyle\boldsymbol{\alpha}_{u}^{*}=\arg\,\underset{\boldsymbol{\alpha}_{u}}{\min}\,\mathcal{J}\left({\bf
x}^{(u)}_{k}+\widehat{\bf B}^{1/2}_{k}\cdot\boldsymbol{\alpha}_{u}\right)\,.$
(35)

A solution of this problem can be approximated as follows: consider the
linearization of $\mathcal{H}_{k}$ about ${\bf x}^{(u)}_{k}$, at iteration $u$

$\displaystyle\displaystyle\mathcal{H}_{k}\left({\bf
x}_{k}\right)\approx\widehat{\mathcal{H}}^{(u)}_{k}\left({\bf x}_{k}\right)$
$\displaystyle=$ $\displaystyle\mathcal{H}_{k}\left({\bf
x}^{(u)}_{k}\right)+{\bf H}^{(u)}_{k}\cdot\left[{\bf x}_{k}-{\bf
x}^{(u)}_{k}\right]$ $\displaystyle=$ $\displaystyle\mathcal{H}_{k}\left({\bf
x}^{(u)}_{k}\right)+{\bf H}^{(u)}_{k}\cdot\widehat{\bf
B}^{1/2}_{k}\cdot\boldsymbol{\alpha}_{u}\,,$

where ${\bf H}^{(u)}_{k}\in\mathbb{R}^{m\times n}$ is the Jacobian of
$\mathcal{H}_{k}({\bf x}_{k})$ at ${\bf x}^{(u)}_{k}$, then, at iteration $u$,
we make use of the next quadratic approximation:

$\displaystyle\displaystyle\widehat{\mathcal{J}}\left(\boldsymbol{\alpha}_{u}\right)$
$\displaystyle\approx$
$\displaystyle\widetilde{\mathcal{J}}(\boldsymbol{\alpha}_{u})$
$\displaystyle=$
$\displaystyle\frac{1}{2}\cdot\left\|\boldsymbol{\beta}_{u}+\boldsymbol{\alpha}_{u}\right\|^{2}+\frac{1}{2}\cdot\sum_{k=0}^{M-1}\left\|{\bf
y}_{k}-\widehat{\mathcal{H}}^{(u)}_{k}\left({\bf x}^{(u)}_{k}+\widehat{\bf
B}^{1/2}_{k}\cdot\boldsymbol{\alpha}_{u}\right)\right\|^{2}_{{\bf
R}_{k}^{-1}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\cdot\left\|\boldsymbol{\beta}_{u}+\boldsymbol{\alpha}_{u}\right\|^{2}+\frac{1}{2}\cdot\sum_{k=0}^{M-1}\left\|{\bf
d}^{(u)}_{k}-{\bf Q}_{k}^{(u)}\cdot\boldsymbol{\alpha}_{u}\right\|^{2}_{{\bf
R}_{k}^{-1}}$

where ${\bf d}^{(u)}_{k}={\bf y}_{k}-\mathcal{H}_{k}\left({\bf
x}_{k}^{(u)}\right)\in\mathbb{R}^{m\times 1}$, and ${\bf Q}_{k}^{(u)}={\bf
H}^{(u)}_{k}\cdot\widehat{\bf B}^{1/2}_{k}\in\mathbb{R}^{m\times n}$. Thus,

$\displaystyle\displaystyle\nabla\widetilde{\mathcal{J}}_{\boldsymbol{\alpha}_{u}}(\boldsymbol{\alpha}_{u})$
$\displaystyle=$
$\displaystyle\boldsymbol{\beta}_{u}+\boldsymbol{\alpha}_{u}-\sum_{k=0}^{M}\left[{\bf
Q}_{k}^{(u)}\right]^{T}\cdot{\bf R}^{-1}_{k}\cdot\left[{\bf d}^{(u)}_{k}-{\bf
Q}_{k}^{(u)}\cdot\boldsymbol{\alpha}_{u}\right]$ $\displaystyle=$
$\displaystyle\left[{\bf I}+\sum_{k=0}^{M}\left[{\bf
Q}_{k}^{(u)}\right]^{T}\cdot{\bf R}^{-1}_{k}\cdot{\bf
Q}_{k}^{(u)}\right]\cdot\boldsymbol{\alpha}_{u}$ $\displaystyle-$
$\displaystyle\left[-\boldsymbol{\beta}_{u}+\sum_{k=0}^{M}\left[{\bf
Q}_{k}^{(u)}\right]^{T}\cdot{\bf R}^{-1}_{k}\cdot{\bf
d}^{(u)}_{k}\right]\in\mathbb{R}^{N\times 1}\,,$

and

$\displaystyle\displaystyle\nabla\widetilde{\mathcal{J}}_{\boldsymbol{\alpha}_{u},\boldsymbol{\alpha}_{u}}(\boldsymbol{\alpha}_{u})$
$\displaystyle=$ $\displaystyle{\bf I}+\sum_{k=0}^{M}\left[{\bf
Q}_{k}^{(u)}\right]^{T}\cdot{\bf R}^{-1}_{k}\cdot{\bf Q}_{k}^{(u)}\,.$

We consider the search direction as the solution of
$\nabla\widetilde{\mathcal{J}}_{\boldsymbol{\alpha}_{u}}(\hat{\boldsymbol{\alpha}}_{u})={\boldsymbol{0}}$,
for computing partial analysis increments onto the control space, this is

$\displaystyle\hat{\boldsymbol{\alpha}}_{u}=\left[\nabla\widetilde{\mathcal{J}}_{\boldsymbol{\alpha}_{u},\boldsymbol{\alpha}_{u}}(\hat{\boldsymbol{\alpha}}_{u})\right]^{-1}\cdot\left[-\boldsymbol{\beta}_{u}+\sum_{k=0}^{M}\left[{\bf
Q}_{k}^{(u)}\right]^{T}\cdot{\bf R}^{-1}_{k}\cdot{\bf d}^{(u)}_{k}\right]\,,$

but, since we employ a linear approximation of the observation operator to
approximate a descent direction of (3) , we perform a line search to estimate
the partial analysis increments:

$\displaystyle\displaystyle\rho^{*}_{u}=\arg\,\underset{\rho_{u}\in[0,\,1]}{\min}\,\mathcal{J}\left({\bf
x}^{(u)}_{k}+\widehat{\bf
B}^{1/2}_{k}\cdot\left[\rho_{u}\cdot\hat{\boldsymbol{\alpha}}_{u}\right]\right)\,,$
(36)

and therefore, the next iterates reads:

$\displaystyle{\bf x}^{(u+1)}_{k}={\bf x}^{(u)}_{k}+{\bf
B}^{1/2}_{k}\cdot\left[\rho^{*}_{u}\cdot\hat{\boldsymbol{\alpha}}_{u}\right]$
(37)

where
$\boldsymbol{\alpha}_{u}^{*}\approx\rho^{*}_{u}\cdot\hat{\boldsymbol{\alpha}}_{u}$,
this yields to
$\boldsymbol{\beta}_{u}^{*}\approx\sum_{v=1}^{u-1}\,\rho^{*}_{v}\cdot\hat{\boldsymbol{\alpha}}_{v}$.
This process is repeated until some stopping criteria is satisfied (i.e., a
maximum number of iterations $U$ is reached). The resulting analysis state
reads:

$\displaystyle\displaystyle\overline{\bf x}_{0}^{a}=\overline{\bf
x}_{0}^{b}+{\bf
B}^{1/2}_{0}\cdot\sum_{u=1}^{U}\,\rho^{*}_{u}\cdot\hat{\boldsymbol{\alpha}}_{u}\,,$

which is nothing but a linear combination of the approximated search
directions among iterations. Thus, we can approximate
$\boldsymbol{\alpha}^{a}$ in (32) as follows:

$\displaystyle\displaystyle\boldsymbol{\alpha}^{a}\approx\sum_{u=1}^{U}\,\rho^{*}_{u}\cdot\hat{\boldsymbol{\alpha}}_{u}\,.$

We can employ any line-search rule to solve the optimization problem (36).
Putting it all together, the optimization steps are detailed in the Algorithm
1. The posterior ensemble onto the control space can then be built by using a
square root approximation of the information matrix, its can be easily shown
that the posterior error statistics read:

$\displaystyle\displaystyle\boldsymbol{\alpha}^{a[e]}\sim\mathcal{N}\left(\boldsymbol{\alpha}^{a},\,\left[{\bf
I}+\sum_{k=0}^{M}\left[{\bf Q}_{k}^{(U)}\right]^{T}\cdot{\bf
R}_{k}^{-1}\cdot{\bf Q}^{(U)}_{k}\right]^{-1}\right)\,,\text{ for $1\leq e\leq
N$}\,,$ (38)

with corresponding analysis members in the model space:

$\displaystyle{\bf x}^{a[e]}_{0}=\overline{\bf x}^{a}_{0}+\widehat{\bf
B}_{0}^{1/2}\cdot\boldsymbol{\alpha}^{a[e]}\,.$

The initial ensemble members are then propagated in time from which the
optimal (ensemble) trajectory is recovered.

1:function find_initial_analysis_ensemble($\left\\{{\bf y}_{k},\overline{\bf
x}_{k}^{b},\widehat{\bf B}^{1/2}_{k},{\bf R}^{-1}_{k}\right\\}^{\bf M}_{k=0}$)

2: $\boldsymbol{\alpha}^{a}\leftarrow{\boldsymbol{0}}$

3: for $k\leftarrow 1\to M$ do

4: ${\bf x}^{(0)}_{k}\leftarrow\overline{\bf x}_{k}^{b}$

5: end for

6: for $u\leftarrow 0\to U$ do

7: $\boldsymbol{\beta}^{*}_{u}\leftarrow\boldsymbol{\alpha}^{a}$

8: for $k\leftarrow 1\to M$ do

9: ${\bf d}^{(u)}_{k}\leftarrow{\bf y}_{k}-\mathcal{H}_{k}\left({\bf
x}_{k}^{(u)}\right)$

10: ${\bf Q}^{(u)}_{k}\leftarrow{\bf H}^{(u)}_{k}\cdot\widehat{\bf
B}^{1/2}_{k}$ $\triangleright$ ${\bf H}^{(u)}_{k}$ is the Jacobian of
$\mathcal{H}_{k}({\bf x}_{k})$ at ${\bf x}^{(u)}_{k}$

11: end for

12:
$\widetilde{\mathcal{J}}_{\boldsymbol{\alpha}_{u},\boldsymbol{\alpha}_{u}}(\hat{\boldsymbol{\alpha}}_{u})\leftarrow{\bf
I}+\sum_{k=0}^{M}\left[{\bf Q}_{k}^{(u)}\right]^{T}\cdot{\bf
R}^{-1}_{k}\cdot{\bf Q}_{k}^{(u)}$

13:
$\hat{\boldsymbol{\alpha}}_{u}\leftarrow\left[\nabla\widetilde{\mathcal{J}}_{\boldsymbol{\alpha}_{u},\boldsymbol{\alpha}_{u}}(\hat{\boldsymbol{\alpha}}_{u})\right]^{-1}\cdot\left[-\boldsymbol{\beta}_{u}+\sum_{k=0}^{M}\left[{\bf
Q}_{k}^{(u)}\right]^{T}\cdot{\bf R}^{-1}_{k}\cdot{\bf d}^{(u)}_{k}\right]$

14:
$\rho^{*}_{u}\leftarrow\arg\,\underset{\rho_{u}\in[0,\,1]}{\min}\,\mathcal{J}\left({\bf
x}^{(u)}_{k}+\widehat{\bf
B}^{1/2}_{k}\cdot\left[\rho_{u}\cdot\hat{\boldsymbol{\alpha}}_{u}\right]\right)$

15: for $k\leftarrow 1\to M$ do

16: ${\bf x}^{(u+1)}_{k}\leftarrow{\bf x}^{(u)}_{k}+{\bf
B}^{1/2}_{k}\cdot\left[\rho^{*}_{u}\cdot\hat{\boldsymbol{\alpha}}_{u}\right]$

17: end for

18:
$\boldsymbol{\alpha}_{u}^{*}\leftarrow\rho^{*}_{u}\hat{\boldsymbol{\alpha}}_{u}$

19:
$\boldsymbol{\alpha}^{a}\leftarrow\boldsymbol{\alpha}^{a}+\boldsymbol{\alpha}_{u}^{*}$

20: end for

21: $\overline{\bf x}_{0}^{a}\leftarrow{\bf x}^{(U)}_{0}$

22: for $e\leftarrow 1\to N$ do

23:
$\boldsymbol{\alpha}^{a[e]}\sim\mathcal{N}\left(\boldsymbol{\alpha}^{a},\,\left[{\bf
I}+\sum_{k=0}^{M}\left[{\bf Q}_{k}^{(U)}\right]^{T}\cdot{\bf
R}_{k}^{-1}\cdot{\bf Q}^{(U)}_{k}\right]^{-1}\right)$

24: ${\bf x}^{a[e]}_{0}\leftarrow\overline{\bf x}^{a}_{0}+\widehat{\bf
B}_{0}^{1/2}\cdot\boldsymbol{\alpha}^{a[e]}$

25: end for

26: return $\\{{\bf x}^{a[e]}_{0}\\}^{N}_{e=1}$

27:end function

Algorithm 1 An Adjoint-Free 4D-Var Method Via Modified Cholesky Decomposition
(4D-Var-MC)

Matrix inversions in eht 4D-Var-MC can be avoided by considering at each
optimization step the matrix-free 4D-Var formulation [36].

###  4.3 An Adjoint-Free 4D-Var Method Via MLEF (4D-Var-MLEF)

The proposed 4D-Var-MC can be similarly derived onto the ensemble space by
employing the matrix of anomalies (7) as our set of basis vectors at
observation times. We then can propose iterates of the form:

$\displaystyle\displaystyle{\bf x}^{(u+1)}_{k}={\bf
x}^{(u)}_{k}+\boldsymbol{\Delta}{\bf
X}_{k}\cdot\boldsymbol{\gamma}_{u}\,,\text{ for $0\leq u\leq U$}\,,$ (39)

where $\boldsymbol{\gamma}_{u}\in\mathbb{R}^{N\times 1}$. We perform analogous
derivations as in (34)–(4.2), replacing the control space $\widehat{\bf
B}^{1/2}_{k}$ by the ensemble space $\boldsymbol{\Delta}{\bf X}_{k}$. Search
directions are then approximated by considering linearizations of the form:

$\displaystyle\displaystyle\mathcal{H}_{k}\left({\bf
x}_{k}\right)\approx\widehat{\mathcal{H}}^{(u)}_{k}\left({\bf x}_{k}\right)$
$\displaystyle=$ $\displaystyle\mathcal{H}_{k}\left({\bf
x}^{(u)}_{k}\right)+{\bf H}^{(u)}_{k}\cdot\boldsymbol{\Delta}{\bf
X}_{k}\cdot\boldsymbol{\alpha}_{u}\,,$

and similar to (36), a line-search method can be employed to control step
lengths among iterations. This algorithm becomes an adjoint-free method based
on the MLEF method discussed in section 3.1; we name it the 4D-Var-MLEF. Note
that, since all computations are performed onto the ensemble space, the
computational cost of estimating analysis increments together with analysis
updates reads:

$\displaystyle\displaystyle\mathcal{O}\left(N^{3}+N^{2}\cdot m+N\cdot m\cdot
n\right)\,.$

As a drawback, the 4D-Var-MLEF suffers from the MLEF limitations regarding
global convergence.

###  4.4 Global Convergence of the 4D-Var-MC

For proving the convergence of the line-search method during the assimilation
of observations in the 4D-Var-MC, we consider the model (1), (2), and we
assume

$\displaystyle\displaystyle\nabla\mathcal{J}\left({\bf
x}^{(u)}\right)^{T}\cdot\left[\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right]<0,\,\text{ for $1\leq
k\leq K$}\,.$ (40)

The next Theorem states the necessary conditions to ensure global convergence
for the line search method during the assimilation of observations in the
4D-Var-MC method.

######  Theorem 4.1.

If (1), (2), and (40) hold, then the line search in the 4D-Var-MC generates an
infinite sequence $\left\\{{\bf x}^{(u)}\right\\}_{u=0}^{\infty}$, then

$\displaystyle\displaystyle\underset{u\rightarrow\infty}{\lim}\left[\frac{-\nabla\mathcal{J}\left({\bf
x}^{(u)}\right)^{T}\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}}{\left\|\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right\|}\right]^{2}=0$ (41)

holds.

###### Proof.

By Taylor series and the Mean Value Theorem we know that,

$\displaystyle\mathcal{J}\left({\bf x}^{(u)}_{0}+\rho_{u}^{*}\cdot{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right)$ $\displaystyle=$
$\displaystyle\mathcal{J}\left({\bf x}^{(u)}_{0}\right)$ $\displaystyle+$
$\displaystyle\rho^{*}\cdot\int_{0}^{1}\,\nabla\mathcal{J}\left({\bf
x}^{(u)}_{0}+\rho^{*}_{u}\cdot t\cdot{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right)^{T}$ $\displaystyle\cdot$
$\displaystyle{\bf B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\cdot dt\,,$

where $\boldsymbol{\alpha}^{*}_{u}$ is given by (35), and therefore,

$\displaystyle\mathcal{J}\left({\bf x}^{(u)}_{0}\right)-\mathcal{J}\left({\bf
x}^{(u+1)}_{0}\right)$ $\displaystyle\geq$
$\displaystyle-\rho^{*}_{u}\cdot\int_{0}^{1}\,\nabla\mathcal{J}\left({\bf
x}^{(u)}_{0}+\rho^{*}_{u}\cdot t\cdot{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right)^{T}$ $\displaystyle\cdot$
$\displaystyle{\bf Q}^{(u)}\cdot{\bf w}^{*}\cdot dt$

for any ${\bf x}^{(u+1)}$ on the ray ${\bf x}^{(u)}+\rho_{u}\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}$, with $\rho_{u}\in[0,\,1]$, we
have

$\displaystyle\displaystyle\mathcal{J}\left({\bf
x}^{(u)}\right)-\mathcal{J}\left({\bf
x}^{(u+1)}\right)\geq\mathcal{J}\left({\bf
x}^{(u)}\right)-\mathcal{J}\left({\bf x}^{(u)}+\rho^{*}_{u}\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right)\,,$

hence:

$\displaystyle\displaystyle\mathcal{J}\left({\bf x}^{(u)}\right)$
$\displaystyle-$ $\displaystyle\mathcal{J}\left({\bf
x}^{(u+1)}\right)\geq-\rho^{*}_{u}\cdot\nabla\mathcal{J}\left({\bf
x}^{(u)}\right)^{T}\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}_{u}^{*}$ $\displaystyle-$
$\displaystyle\rho^{*}_{u}\cdot\int_{0}^{1}\,\left[\nabla\mathcal{J}\left({\bf
x}^{(u)}+\rho^{*}_{u}\cdot t\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right)-\nabla\mathcal{J}\left({\bf
x}^{(u)}\right)\right]^{T}$ $\displaystyle\cdot$ $\displaystyle\widehat{\bf
B}^{1/2}\cdot\boldsymbol{\alpha}^{*}_{u}\cdot dt\,,$

by the Cauchy Schwarz inequality, we have

$\displaystyle\displaystyle\mathcal{J}\left({\bf x}^{(u)}\right)$
$\displaystyle-$ $\displaystyle\mathcal{J}\left({\bf
x}^{(u+1)}\right)\geq-\rho^{*}_{u}\cdot\nabla\mathcal{J}\left({\bf
x}^{(u)}\right)^{T}\cdot\widehat{\bf B}^{1/2}\cdot\boldsymbol{\alpha}^{*}_{u}$
$\displaystyle-$
$\displaystyle\rho^{*}_{u}\cdot\int_{0}^{1}\,\left\|\nabla\mathcal{J}\left({\bf
x}^{(u)}+\rho^{*}_{u}\cdot t\cdot\widehat{\bf
B}^{1/2}\cdot\boldsymbol{\alpha}^{*}_{u}\right)-\nabla\mathcal{J}\left({\bf
x}^{(u)}\right)\right\|$ $\displaystyle\cdot$
$\displaystyle\left\|\widehat{\bf
B}^{1/2}\cdot\boldsymbol{\alpha}^{*}_{u}\right\|\cdot dt$ $\displaystyle\geq$
$\displaystyle-\rho^{*}_{u}\cdot\nabla\mathcal{J}\left({\bf
x}^{(u)}\right)^{T}\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}$ $\displaystyle-$
$\displaystyle\rho^{*}_{u}\cdot\int_{0}^{1}\,L\cdot\left\|\rho^{*}_{u}\cdot
t\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right\|\cdot\left\|\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}_{u}^{*}\right\|\cdot dt$ $\displaystyle=$
$\displaystyle-\rho^{*}_{u}\cdot\nabla\mathcal{J}\left({\bf
x}^{(u)}\right)^{T}\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}$ $\displaystyle-$
$\displaystyle\rho^{*}_{u}\cdot L\cdot\left\|\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right\|\cdot\int_{0}^{1}\left\|t\cdot\rho^{*}_{u}\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right\|\cdot dt$ $\displaystyle=$
$\displaystyle-\rho^{*}_{u}\cdot\nabla\mathcal{J}\left({\bf
x}^{(u)}\right)^{T}\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}-\frac{1}{2}\cdot{\rho^{*}_{u}}^{2}\cdot
L\cdot\left\|\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right\|^{2}\,,$

since the gradient is Lipschitz continuous, we can choose

$\displaystyle\displaystyle\rho^{*}_{u}=-\frac{\nabla\mathcal{J}\left({\bf
x}^{(u)}\right)^{T}\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}}{L\cdot\left\|\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right\|^{2}}\,,$

therefore,

$\displaystyle\displaystyle\mathcal{J}\left({\bf x}^{(u)}\right)$
$\displaystyle-$ $\displaystyle\mathcal{J}\left({\bf
x}^{(u+1)}\right)\geq\frac{\left[\nabla\mathcal{J}\left({\bf
x}^{(u)}\right)^{T}\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right]^{2}}{L\cdot\left\|\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right\|^{2}}$ $\displaystyle-$
$\displaystyle\frac{1}{2}\cdot\frac{\left[-\nabla\mathcal{J}\left({\bf
x}^{(u)}\right)^{T}\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right]^{2}}{L\cdot\left\|\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right\|^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{2\cdot L}\cdot\left[-\frac{\nabla\mathcal{J}\left({\bf
x}^{(u)}\right)^{T}\cdot\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}}{\left\|\widehat{\bf
B}^{1/2}_{0}\cdot\boldsymbol{\alpha}^{*}_{u}\right\|}\right]^{2}\,.$

By (1), and (40), it follows that $\left\\{\mathcal{J}\left({\bf
x}^{(u)}\right)\right\\}_{u=0}^{\infty}$ is a monotone decreasing number
sequence and it has a bound below, therefore $\left\\{\mathcal{J}\left({\bf
x}^{(u)}\right)\right\\}_{u=0}^{\infty}$ has a limit, and consequently (41)
holds. ∎

##  5 Experimental Settings

We make use of the Lorenz-96 model [31] as our surrogate model in order to
tests the compared EnKF implementations. The Lorenz-96 model is described by
the following set of ordinary differential equations [10]:

$\displaystyle\displaystyle\frac{dx_{j}}{dt}=\begin{cases}\left(x_{2}-x_{n-1}\right)\cdot
x_{n}-x_{1}+F&\text{ for $j=1$},\\\ \left(x_{j+1}-x_{j-2}\right)\cdot
x_{j-1}-x_{j}+F&\text{ for $2\leq j\leq n-1$},\\\
\left(x_{1}-x_{n-2}\right)\cdot x_{n-1}-x_{n}+F&\text{ for $j=n$},\end{cases}$
(42)

where $F$ is an external force and $n=40$ is the number of model components.
Periodic boundary conditions are assumed. When $F=8$ units, the model exhibits
chaotic behavior, which makes it a relevant surrogate problem for atmospheric
dynamics [22, 12]. A time unit in the Lorenz-96 represents 120 hours (5 days)
in the atmosphere.

The experimental settings are described below, they are similar to those
presented in [35]:

  * •

An initial random solution is integrated over a long time period in order to
obtain an initial condition ${\bf x}^{*}_{-2}\in\mathbb{R}^{n\times 1}$
dynamically consistent with the model (42).

  * •

A perturbed background solution ${\bf x}^{b}_{-2}$ is obtained at time
$t_{-2}$ by drawing a sample from the Normal distribution,

$\displaystyle\displaystyle{\bf x}^{b}_{-2}\sim\mathcal{N}\left({\bf
x}^{*}_{-2},\,0.05^{2}\cdot{\bf I}\right)\,,$

this solution is then integrated for 10 time units (equivalent to 50 days in
the atmosphere) in order to obtain a background solution ${\bf x}^{b}_{-1}$
consistent with the numerical model.

  * •

An initial perturbed ensemble is built about the background state by taking
samples from the distribution,

$\displaystyle\displaystyle{\bf x}^{b[i]}_{-1}\sim\mathcal{N}\left({\bf
x}^{b}_{-1},\,0.05^{2}\cdot{\bf I}\right)\,,\text{ for $1\leq i\leq N$}\,.$

These samples are propagated over a time period of 10 units in order to make
them consistent with the model dynamics (42), after which the initial ensemble
members ${\bf x}^{b[i]}_{0}\in\mathbb{R}^{n\times 1}$ are obtained. We build
an initial pool of members. Two-dimensional projections of such collection are
shown in figure 1.

(a) 2D projection of the inital members.

(b) 2D initial background error distribution.

Figure 1: Initial pool of background members for the experiments. Two
dimensional projections based on the two leading directions are shown.

  * •

The number of assimilation steps reads $M=500$. Observations are collected
every 12 hours and their error statistics are described by the Gaussian
distribution. We use the Dormand and Prince method [8] with absolute error
tolerance of $10^{-7}$ in order to propagate the numerical model every 0.5
time units,

$\displaystyle\displaystyle{\bf y}_{k}\sim\mathcal{N}\left({\bf
H}_{k}\cdot{\bf x}^{*}_{k}\,,0.01^{2}\cdot{\bf I}\right)\,,\text{ for $1\leq
k\leq M$}\,.$ (43)

We try different number of observations per assimilation step, these are given
by

$\displaystyle\displaystyle m=p\cdot n\,,$

where $p=\\{70\%\,,100\%\\}$, which corresponds with $m\in\\{28,\,40\\}$
observations. Observations from the model state are randomly chosen at the
different assimilation steps. We employ the non-linear observation operator
[51]:

$\displaystyle\displaystyle\left\\{\mathcal{H}\left({\bf
x}\right)\right\\}_{j}=\frac{\left\\{{\bf
x}\right\\}_{j}}{2}\cdot\left[\left(\frac{\left|\left\\{{\bf
x}\right\\}_{j}\right|}{2}\right)^{\gamma-1}+1\right]\,,$ (44)

where $j$ denotes the $j$-th observed component from the model state.
Likewise, $1\leq\gamma\leq 7$. Note that, the experiments are performed under
perfect model assumptions.

  * •

We evenly range the inflation factor $\rho\in[1.1,\,1.9]$ with increments of
$0.1$.

  * •

Different ensemble sizes are tried during the experiments,
$N\in\\{20,\,40\\}$.

  * •

The L–2 norm of errors is utilized as a measure of accuracy, at the
assimilation step $k$, it is defined as follows:

$\displaystyle\lambda_{k}=\left\|{\bf x}^{*}_{k}-\overline{\bf
x}^{a}_{k}\right\|_{2}=\sqrt{\left[{\bf x}^{*}_{k}-\overline{\bf
x}^{a}_{k}\right]^{T}\cdot\left[{\bf x}^{*}_{k}-\overline{\bf
x}^{a}_{k}\right]}\,,$ (45)

where ${\bf x}^{*}_{k}$ and ${\bf x}^{a}_{k}$ are the reference and the
analysis solutions, respectively.

  * •

The Root-Mean-Square-Error (RMSE) measures the performance of a filter, in
average, for a given number of assimilation steps $M$:

$\displaystyle\displaystyle\lambda=\sqrt{\frac{1}{M}\cdot\sum_{k=1}^{M}\lambda_{k}^{2}}\,.$

  * •

The number of optimization steps $U$ equals 10. We choose this value for two
main reasons: model propagations are involved during optimization steps and
even more, for either of the compared filter implementations, the 4D-EnKF-MC
or the 4D-MLEF, after 10 iterations no reduction in cost-function values can
be evidenced.

For each configuration of the parameters, we perform 30 experiments to get
insights about the error behavior and his associate statistics. We employ the
first two statistical moments (mean and variance) of underlying error
distributions to analyze the numerical results. The analysis is divided into
two groups, (nearly) linear observation operators where $1\leq\gamma\leq 3$
and highly non-linear observation operators $4\leq\gamma\leq 7$.

###  5.1 Results with (non-) linear observation operators $1\leq\gamma\leq 3$

For different values of $N$ and $1\leq\gamma\leq 3$, the results are shown in
figures 2 and 3 for $p=70\%$ and $p=100\%$ (full observational networks),
respectively. We set $r=2$ to build the set of basis vectors in the 4D-EnKF-MC
context. As expected, the accuracy of both filter implementations, the
4D-EnKF-MC, and the 4D-MLEF, are better than that of the NO Data Assimilation
(NODA) forecast, which is reported as well. Note that increments in the
ensemble size have a clear impact on the 4D-EnKF-MC performance; the larger
the ensemble size, the better (in terms of accuracy) the posterior estimate,
for instance, the L–2 error norms are decreased by some order of magnitudes.
This behavior is not evident in the 4D-MLEF method wherein slight improvements
can be evidenced. Note that both filter implementations benefit from having
more observations to digest during assimilation steps. This is, errors can be
decreased as the number of observations is increased if the assimilation
window keeps constant. Again, this patron is more evident in the 4D-EnKF-MC
method than in the 4D-MLEF filter. This can be explained as follows: during
assimilation steps, the 4D-EnKF-MC method builds control spaces whose
dimensions equal the model one. Thus, all information brought by the
observations can be captured onto that space (whose dimension is larger than
the observation one), not the 4D-MLEF whose space dimension relies on the
ensemble size.

###  5.2 Results with highly non-linear observation operators $4\leq\gamma\leq
7$

For different values of $N$ and $4\leq\gamma\leq 7$, the results are shown in
figures 4 and 5 for $p=70\%$ and $p=100\%$ (full observational networks),
respectively. Again, we set $r=2$ to build the set of basis vectors in the
4D-EnKF-MC context. As expected, all filter implementations provide better
results than those obtained by the NODA forecast. In all cases, the results
proposed by the 4D-EnKF-MC are better than those of the 4D-MLEF in terms of
L–2 error norms. For $p=70\%$, the results of the compared filter
implementations get closer as the degree $\gamma$ increases. Nevertheless, as
the number of observations increases, the 4D-EnKF-MC provides better results
than those of the 4D-MLEF. Note that, in the 4D-EnKF-MC, the quality of
analysis corrections improves as the number of observations increases.
Increments in the ensemble size have a clear impact in the quality of both
filters, the larger the ensemble size, the better the quality of analysis
corrections.

###  5.3 Further discussions

For all parameter configurations, the results are shown in the Tables 1 and 2
for $p=70\%$ and $p=100\%$, respectively. In the 4D-EnKF-MC method, besides
$r=2$, the values of 6 and 18 are tried as the radius of influence for
building the set of basis vectors during assimilation steps, in this manner,
we can have an idea about how the radius size impacts the quality of the
initial analysis increments. The first thing to note is that, as the radius of
influence increases, the initial analysis corrections’ quality degrades. This
can be expected given the special form of the error dynamics in the numerical
model (42) wherein fast decay of error correlations can be possible. For
instance, the model dynamics strongly correlate errors among neighboring
components. Thus, having large radii of influences can impact the quality of
posterior initial members. Note that the behavior of 4D-MLEF is similar for
all $\gamma$ values, while for the 4D-EnKF-MC RMSE, values degrade as the
$\gamma$ increases. However, for all scenarios, the results of the 4D-EnKF-MC
are better than those of the 4D-MLEF. In both filter implementations, their
performance can be improved by increasing the number of ensemble members. For
instance, in the 4D-MLEF, the larger the ensemble size, the more degrees of
freedom are available during optimization steps. In the 4D-EnKF-MC method,
under such circumstances, more information about the error dynamics is
available to fit the models (12).

| $N=20$ | $N=60$
---|---|---
$\gamma=1$ |  |
$\gamma=2$ |  |
$\gamma=3$ |  |
Figure 2: Averages (dashed lines) and standard deviations (shaded regions) of error norms for $p=70\%$. We let $N\in\\{20,\,60\\}$. Likewise, $\gamma\in\\{1,\,2,\,3\\}$. | $N=20$ | $N=60$
---|---|---
$\gamma=1$ |  |
$\gamma=2$ |  |
$\gamma=3$ |  |
Figure 3: Averages (dashed lines) and standard deviations (shaded regions) of error norms for 30 realizations of experiments and the values of parameters $p=100\%,N\in\\{20,\,60\\}$, and $1\leq\gamma\leq 3$. | $N=20\%$ | $N=60\%$
---|---|---
$\gamma=4$ |  |
$\gamma=5$ |  |
$\gamma=6$ |  |
$\gamma=7$ |  |
Figure 4: Averages (dashed lines) and standard deviations (shaded regions) of error norms for 30 realizations of experiments and the values of parameters $p=70\%,N\in\\{20,\,60\\}$, and $4\leq\gamma\leq 7$. The results are shown for the compared filter implementations as well as the NO Data Assimilation (NODA) forecast. | $N=20\%$ | $N=60\%$
---|---|---
$\gamma=4$ |  |
$\gamma=5$ |  |
$\gamma=6$ |  |
$\gamma=7$ |  |
Figure 5: Averages (dashed lines) and standard deviations (shaded regions) of error norms for 30 realizations of experiments and the values of parameters $p=100\%,N\in\\{20,\,60\\}$, and $4\leq\gamma\leq 7$. $\gamma$ | $N$ | $r$ | 4D-EnKF-MC | 4D-MLEF | NODA
---|---|---|---|---|---
1 | 20 | 2 | 0.158 | 22.397 | 31.328
6 | 0.549 | 22.967
18 | 18.233 | 19.743
60 | 2 | 0.144 | 22.206
6 | 0.173 | 22.167
18 | 0.629 | 22.697
2 | 20 | 2 | 0.276 | 22.944 | 31.387
6 | 2.172 | 26.318
18 | 18.725 | 19.968
60 | 2 | 0.148 | 22.609
6 | 0.234 | 23.261
18 | 1.665 | 25.636
3 | 20 | 2 | 11.230 | 23.642 | 31.387
6 | 14.081 | 22.601
18 | 20.093 | 20.985
60 | 2 | 7.730 | 22.077
6 | 10.534 | 21.661
18 | 12.979 | 21.499
5 | 20 | 2 | 20.736 | 24.138 | 31.387
6 | 20.065 | 24.411
18 | 20.970 | 22.574
60 | 2 | 19.225 | 21.930
6 | 18.550 | 21.596
18 | 18.736 | 21.772
6 | 20 | 2 | 22.629 | 24.236 | 31.387
6 | 21.977 | 23.924
18 | 21.386 | 23.040
60 | 2 | 20.540 | 22.365
6 | 20.449 | 22.200
18 | 20.927 | 21.826
7 | 20 | 2 | 23.209 | 23.837 | 31.387
6 | 22.368 | 23.527
18 | 21.650 | 23.476
60 | 2 | 20.947 | 22.432
6 | 21.088 | 22.275
18 | 22.512 | 22.551
Table 1: Averages of RMSE values for 30 realizations of experiments and the parameter values: $p=70\%,N\in\\{20,\,60\\}$, and $1\leq\gamma\leq 7$. $\gamma$ | $N$ | $r$ | 4D-EnKF-MC | 4D-MLEF | NODA
---|---|---|---|---|---
1 | 20 | 2 | 0.143 | 22.396 | 31.463
6 | 0.520 | 22.948
18 | 17.963 | 19.680
60 | 2 | 0.110 | 22.159
6 | 0.145 | 22.160
18 | 0.536 | 22.522
2 | 20 | 2 | 0.280 | 22.838 | 31.404
6 | 1.890 | 24.636
18 | 18.591 | 19.748
60 | 2 | 0.131 | 22.349
6 | 0.202 | 22.691
18 | 1.476 | 23.618
3 | 20 | 2 | 11.261 | 23.013 | 31.404
6 | 14.209 | 22.643
18 | 20.046 | 20.998
60 | 2 | 8.117 | 22.882
6 | 10.734 | 21.926
18 | 12.720 | 21.607
5 | 20 | 2 | 20.463 | 23.668 | 31.404
6 | 20.310 | 24.062
18 | 21.019 | 22.643
60 | 2 | 19.090 | 22.293
6 | 18.432 | 21.718
18 | 18.842 | 21.701
6 | 20 | 2 | 22.188 | 23.586 | 31.404
6 | 22.128 | 23.080
18 | 21.343 | 23.024
60 | 2 | 20.197 | 22.196
6 | 20.216 | 22.026
18 | 20.955 | 22.054
7 | 20 | 2 | 23.676 | 24.389 | 31.404
6 | 23.089 | 23.566
18 | 21.626 | 23.394
60 | 2 | 21.098 | 22.340
6 | 21.197 | 22.437
18 | 21.785 | 22.226
Table 2: Averages of RMSE values for 30 realizations of experiments and the
parameter values: $p=100\%,N\in\\{20,\,60\\}$, and $1\leq\gamma\leq 7$.

##  6 Conclusions

We propose two adjoint free 4D-Var methods for non-linear data assimilation.
In both methods, control spaces are built to avoid the intrinsic need for
adjoint models in the 4D-Var context. The techniques work as follows: based on
an ensemble of model realizations, control spaces are built at observation
times, then model states are constrained to the space spanned by those control
spaces. The analysis increments are computed onto the proposed spaces by
linearizing observation operators together with employing line-search
optimization methods. The global convergence of the method can be proven as
long as the control-space dimension equals that of the model one. We also
propose a matrix-free method to compute search directions whose computational
effort is linearly bounded by the model dimension; this feature is attractive
under data assimilation operational settings. Experimental tests are performed
by using the Lorenz-96 model and by employing a non-linear observation
operator. The results reveal that both filters can compute initial posterior
modes of the error distribution with reasonable accuracy in terms of RMSE
values. Even more, the 4D-EnKF-MC can be employed to mitigate the impact of
sampling noise as well as to guarantee the convergence of line-search methods
during optimization steps.

## Acknowledgments

This work was supported in part by award UN 2018-38, and by the Applied Math
and Computer Science Lab at Universidad del Norte, Colombia .

## Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * [1] Arulampalam, M.S., Maskell, S., J. Neil, G., Clapp, T.: A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking.  Signal Processing, IEEE Transactions on 50, 174 – 188 (2002)
  * [2] Battiti, R.: First-and second-order methods for learning: between steepest descent and newton’s method.  Neural computation 4(2), 141–166 (1992)
  * [3] Bickel, P.J., Levina, E., et al.: Covariance regularization by thresholding.  The Annals of Statistics 36(6), 2577–2604 (2008)
  * [4] Bickel, P.J., Levina, E., et al.: Regularized estimation of large covariance matrices.  The Annals of Statistics 36(1), 199–227 (2008)
  * [5] Chen, Y., Oliver, D.S.: Cross-covariances and localization for enkf in multiphase flow data assimilation.  Computational Geosciences 14(4), 579–601 (2010)
  * [6] Conn, A.R., Gould, N.I., Toint, P.L.: Trust region methods, vol. 1.  Siam (2000)
  * [7] Curtis, F.E., Robinson, D.P., Samadi, M.: A trust region algorithm with a worst-case iteration complexity of for nonconvex optimization.  Mathematical Programming 162(1-2), 1–32 (2017)
  * [8] Dormand, J.R., Prince, P.J.: A family of embedded runge-kutta formulae.  Journal of computational and applied mathematics 6(1), 19–26 (1980)
  * [9] Evensen, G.: The ensemble kalman filter: Theoretical formulation and practical implementation.  Ocean dynamics 53(4), 343–367 (2003)
  * [10] Fertig, E.J., Harlim, J., Hunt, B.R.: A comparative study of 4d-var and a 4d ensemble kalman filter: Perfect model simulations with lorenz-96.  Tellus A 59(1), 96–100 (2007)
  * [11] Fletcher, R., Reeves, C.M.: Function minimization by conjugate gradients.  The computer journal 7(2), 149–154 (1964)
  * [12] Gottwald, G.A., Melbourne, I.: Testing for chaos in deterministic systems with noise.  Physica D: Nonlinear Phenomena 212(1), 100–110 (2005)
  * [13] Greybush, S.J., Kalnay, E., Miyoshi, T., Ide, K., Hunt, B.R.: Balance and ensemble kalman filter localization techniques.  Monthly Weather Review 139(2), 511–522 (2011)
  * [14] Grippo, L., Lampariello, F., Lucidi, S.: A nonmonotone line search technique for newton’s method.  SIAM Journal on Numerical Analysis 23(4), 707–716 (1986)
  * [15] Grippo, L., Lampariello, F., Lucidi, S.: A truncated newton method with nonmonotone line search for unconstrained optimization.  Journal of Optimization Theory and Applications 60(3), 401–419 (1989)
  * [16] Hager, W.W., Zhang, H.: A survey of nonlinear conjugate gradient methods.  Pacific journal of Optimization 2(1), 35–58 (2006)
  * [17] Harlim, J., Hunt, B.R.: Four-dimensional local ensemble transform kalman filter: numerical experiments with a global circulation model.  Tellus A: Dynamic Meteorology and Oceanography 59(5), 731–748 (2007)
  * [18] Hosseini, S., Huang, W., Yousefpour, R.: Line search algorithms for locally lipschitz functions on riemannian manifolds.  SIAM Journal on Optimization 28(1), 596–619 (2018)
  * [19] Houtekamer, P.L., Mitchell, H.L.: Data assimilation using an ensemble kalman filter technique.  Monthly Weather Review 126(3), 796–811 (1998)
  * [20] Huang, X.Y., Xiao, Q., Barker, D.M., Zhang, X., Michalakes, J., Huang, W., Henderson, T., Bray, J., Chen, Y., Ma, Z., et al.: Four-dimensional variational data assimilation for wrf: Formulation and preliminary results.  Monthly Weather Review 137(1), 299–314 (2009)
  * [21] Kalnay, E.: Atmospheric modeling, data assimilation and predictability.  Cambridge university press (2003)
  * [22] Karimi, A., Paul, M.R.: Extensive chaos in the lorenz-96 model.  Chaos: An Interdisciplinary Journal of Nonlinear Science 20(4), 043105 (2010)
  * [23] Knoll, D.A., Keyes, D.E.: Jacobian-free newton–krylov methods: a survey of approaches and applications.  Journal of Computational Physics 193(2), 357–397 (2004)
  * [24] Lahoz, B.K.W., Menard, R.: Data assimilation.  Springer (2010)
  * [25] Levina, E., Rothman, A., Zhu, J., et al.: Sparse estimation of large covariance matrices via a nested lasso penalty.  The Annals of Applied Statistics 2(1), 245–263 (2008)
  * [26] Lewis, R.M., Torczon, V., Trosset, M.W.: Direct search methods: then and now.  Journal of computational and Applied Mathematics 124(1-2), 191–207 (2000)
  * [27] Liu, C., Xiao, Q., Wang, B.: An ensemble-based four-dimensional variational data assimilation scheme. part i: Technical formulation and preliminary test.  Monthly Weather Review 136(9), 3363–3373 (2008)
  * [28] Loke, M.H., Barker, R.: Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-newton method.  Geophysical prospecting 44(1), 131–152 (1996)
  * [29] Lorenc, A.C.: The potential of the ensemble kalman filter for nwp—a comparison with 4d-var.  Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography 129(595), 3183–3203 (2003)
  * [30] Lorenc, A.C., Bowler, N.E., Clayton, A.M., Pring, S.R., Fairbairn, D.: Comparison of hybrid-4denvar and hybrid-4dvar data assimilation methods for global nwp.  Monthly Weather Review 143(1), 212–229 (2015)
  * [31] Lorenz, E.N.: Designing chaotic models.  Journal of the Atmospheric Sciences 62(5), 1574–1587 (2005).  DOI 10.1175/JAS3430.1
  * [32] Miyoshi, T., Kunii, M.: The local ensemble transform kalman filter with the weather research and forecasting model: Experiments with real observations.  Pure and applied geophysics 169(3), 321–333 (2012)
  * [33] Moré, J.J., Sorensen, D.C.: Computing a trust region step.  SIAM Journal on Scientific and Statistical Computing 4(3), 553–572 (1983)
  * [34] Moré, J.J., Thuente, D.J.: Line search algorithms with guaranteed sufficient decrease.  ACM Transactions on Mathematical Software (TOMS) 20(3), 286–307 (1994)
  * [35] Nerger, L., Schulte, S., Bunse-Gerstner, A.: On the influence of model nonlinearity and localization on ensemble kalman smoothing.  Quarterly Journal of the Royal Meteorological Society 140(684), 2249–2259 (2014)
  * [36] Nino-Ruiz, E.D., Guzman-Reyes, L.G., Beltran-Arrieta, R.: An adjoint-free four-dimensional variational data assimilation method via a modified cholesky decomposition and an iterative woodbury matrix formula.  Nonlinear Dynamics 99(3), 2441–2457 (2020)
  * [37] Nino-Ruiz, E.D., Sandu, A., Deng, X.: A parallel implementation of the ensemble kalman filter based on modified cholesky decomposition.  Journal of Computational Science (2017)
  * [38] Nino-Ruiz, E.D., Sandu, A., Deng, X.: An ensemble kalman filter implementation based on modified cholesky decomposition for inverse covariance matrix estimation.  SIAM Journal on Scientific Computing 40(2), A867–A886 (2018)
  * [39] Nocedal, J.: Updating quasi-newton matrices with limited storage.  Mathematics of computation 35(151), 773–782 (1980)
  * [40] Pan, V.Y., Branham, S., Rosholt, R.E., Zheng, A.L.: Newton’s iteration for structured matrices.  In: Fast reliable algorithms for matrices with structure, pp. 189–210. SIAM (1999)
  * [41] Reichle, R.H.: Data assimilation methods in the earth sciences.  Advances in water resources 31(11), 1411–1418 (2008)
  * [42] Rothman, A.J., Levina, E., Zhu, J.: Generalized thresholding of large covariance matrices.  Journal of the American Statistical Association 104(485), 177–186 (2009)
  * [43] Ruiz, E.D.N., Sandu, A.: A derivative-free trust region framework for variational data assimilation.  Journal of Computational and Applied Mathematics 293, 164–179 (2016)
  * [44] Savard, G., Gauvin, J.: The steepest descent direction for the nonlinear bilevel programming problem.  Operations Research Letters 15(5), 265–272 (1994)
  * [45] Shanno, D.F.: Conditioning of quasi-newton methods for function minimization.  Mathematics of computation 24(111), 647–656 (1970)
  * [46] Shi, Z.J.: Convergence of line search methods for unconstrained optimization.  Applied Mathematics and Computation 157(2), 393–405 (2004)
  * [47] Stroud, J.R., Katzfuss, M., Wikle, C.K.: A bayesian adaptive ensemble kalman filter for sequential state and parameter estimation.  Monthly Weather Review 146(1), 373–386 (2018)
  * [48] Tr’emolet, Y.: Accounting for an imperfect model in 4d-var.  Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography 132(621), 2483–2504 (2006)
  * [49] Trémolet, Y.: Incremental 4d-var convergence study.  Tellus A: Dynamic Meteorology and Oceanography 59(5), 706–718 (2007)
  * [50] Uschmajew, A., Vandereycken, B.: Line-search methods and rank increase on low-rank matrix varieties.  In: Proceedings of the 2014 International Symposium on Nonlinear Theory and its Applications (NOLTA2014), pp. 52–55 (2014)
  * [51] Van Leeuwen, P.J., Cheng, Y., Reich, S.: Nonlinear data assimilation, vol. 2.  Springer (2015)
  * [52] Vanderplaats, G.N.: Numerical optimization techniques for engineering design: with applications, vol. 1.  McGraw-Hill New York (1984)
  * [53] Wang, X., Hamill, T.M., Whitaker, J.S., Bishop, C.H.: A comparison of hybrid ensemble transform kalman filter–optimum interpolation and ensemble square root filter analysis schemes.  Monthly weather review 135(3), 1055–1076 (2007)
  * [54] Wright, S., Nocedal, J.: Numerical optimization.  Springer Science 35(67-68), 7 (1999)
  * [55] Zhou, W., Akrotirianakis, I., Yektamaram, S., Griffin, J.: A matrix-free line-search algorithm for nonconvex optimization.  Optimization Methods and Software pp. 1–24 (2017)
  * [56] Zupanski, M.: Maximum likelihood ensemble filter: Theoretical aspects.  Monthly Weather Review 133(6), 1710–1726 (2005).  DOI 10.1175/MWR2946.1


# Exploring the Properties and Evolution of Neural Network Eigenspaces during
Training

Mats L. Richter
Department of Cognitive Science
University of Osnabrück
49080, Osnabrück
<EMAIL_ADDRESS>
&Leila Malihi
Department of Cognitive Science
University of Osnabrück
49080, Osnabrück
<EMAIL_ADDRESS>
&Anne-Kathrin Patricia Windler
Department of Cognitive Science
University of Osnabrück
49080, Osnabrück
<EMAIL_ADDRESS>
&Ulf Krumnack
Department of Cognitive Science
University of Osnabrück
49080, Osnabrück
<EMAIL_ADDRESS>

###### Abstract

In this work we explore the information processing inside neural networks
using logistic regression probes Alain and Bengio (2018) and the saturation
metric Richter et al. (2020). We show that problem difficulty and neural
network capacity affect the predictive performance in an antagonistic manner,
opening the possibility of detecting over- and under-parameterization of
neural networks for a given task. We further show that the observed effects
are independent from previously reported pathological patterns like the “tail
pattern” described in Richter et al. (2020). Finally we are able to show that
saturation patterns converge early during training, allowing for a quicker
cycle time during analysis

_Keywords_ convolutional neural networks $\cdot$ logistic regression probes
$\cdot$ saturation

##  1 Introduction

Deep Convolutional Neural Networks are opaque machine learning solutions.
Opaque in the sense that the model state itself is neither self evident nor
human interpretable. This has lead to a primarily trial and error driven
approach for development, that relies on the comparison of abstract, model
agnostic performance metrics like accuracy as a measure for classification
performance, the number of parameters as a measure of capacity and the FLOPs
per forward pass as a metric for computational efficiency Tan et al. (2018);
Tan and Le (2019); He et al. (2015). To move towards a more efficient, less
trial and error based design process, a deeper understanding of the model’s
state is required. This understanding does not have to be necessarily complete
with regards to fully understanding the relation of the input and output of
the model. The comparative analysis methods based on SVCCA Raghu et al. (2017)
are good examples of such a non-holistic approaches. The information extracted
from the model by using SVCCA is highly aggregated but allows for useful
insights into the converged model. Logistic Regression Probes Alain and Bengio
(2018) and saturation Richter et al. (2020) aggregate a single layer to a
number, which allows for easy and intuitive analysis, similar to measuring
with a thermometer. While logistic regression probes measure the intermediate
solution quality very directly by training logistic regressions on the output
of a layer, saturation is more task agnostic. The authors show that the size
of the subspace of the feature space 111We refer to feature space as the space
in which the data exists in a specific layer‘s output. responsible for data
processing varies significantly depending on the input resolution, leading to
inefficiencies. In this work we explore the properties of saturation further
by answering the following questions:

  * •

How is saturation influenced by model capacity and problem difficulty?222We
define capacity of a network by the number of trainable parameters, which is
also common in other works Tan and Le (2019); Real et al. (2019); Huang et al.
(2019).? Answer: Yes, we show that problem difficulty and model capacity
behave inversely proportional to each other in their influence on saturation.

  * •

Does the capacity (number of filters) of individual layers influence how the
inference process is distributed over the network? Answer: No, the tested
models seem to be unable to shift processing to other layers.

  * •

How do saturation patterns evolve during training? Answer: They converge at a
similar pace as the loss of the model, with saturation increasing during
training. The way saturation evolves also gives hints on the properties of the
dataset, but it is not influenced by the model over-fitting.

##  2 Background

###  2.1 Logistic Regression Probes

For the analysis of trained models, we use probe classifiers. Probe
classifiers are a tool for analyzing how the solution quality progresses while
the data is propagated through the network’s structure Alain and Bengio
(2018). To do this, logistic regression "probes" are trained on the same task
as the model, using the output of individual layers as input. Since the
softmax layer and the probes effectively solve the same task, the probes can
be used to judge the quality of the intermediate solutions. Typically, probe
performance increases layer by layer, approaching the model’s performance (as
figure 1 demonstrates). We refer to the test accuracy of a probe computed on
the output of layer $l$ as $p_{l}$

###  2.2 Saturation

The saturation $s_{l}$ of a layer $l$ is a simple scalar metric that can be
applied to any layer in a neural network and was first introduced by Shenk et
al. (2019) an explored in more detail by Richter et al. (2020). In simple
terms, saturation is a percentage value that measures how many of the
available dimensions in the output space $Z_{l}$ of the layer $l$ are relevant
for the inference process.

$s_{l}=\frac{\operatorname{dim}E^{k}_{l}}{\operatorname{dim}Z_{l}}$ (1)

Saturation is computed by approximating the ratio of the dimensionality of the
relevant eigenspace $\operatorname{dim}E^{k}_{l}$ of layer $l$ and the
extrinsic dimensionality of layer $l$s activation values
$\operatorname{dim}Z_{l}$. The relevant eigenspace $E^{k}_{l}$ is a subspace
of $Z_{l}$ in which the information is processed. We refer to this space as
"relevant", since a projection of the data into the relevant eigenspace after
the layers output will not result in a loss of predictive performance Richter
et al. (2020). The relevant eigenspace can be considered the subspace in which
the information processing is happening. The approximation of $E^{k}_{l}$ is
done using PCA, where the largest eigendirections are kept in order to explain
$99\%$ of the data‘s variance in the output of layer $l$. Because of this
computation strategy, saturation can be computed on-line during training. This
is in contrast to logistic regression probes, which generally require the
extraction of all latent representations and training of logistic regression
probes, which may take longer than the actual training itself.

###  2.3 The Semantics of Saturation

A sequence of low saturated layers ($<$ 50% of the average saturation of all
other layers) is referred to as a “tail pattern” and indicates that these
layers are not contributing qualitatively to the prediction (see figure 1).
This is also visible when observing the performance of the probes, which
stagnates in the tail-layers.

Figure 1: VGG16 exhibits a tail pattern starting from Conv8. In the tail,
saturation is low and probe accuracy is stagnating at the same level as the
model output.

From these results, we can see that solving a problem saturates the layer more
than simply passing through information in a given layer. However, this does
not mean that the absolute saturation value is indicative of the activity
within a layer. So far, saturation has been only explored as a quicker on-line
computable analogy for logistic regression probes. As such saturation has been
viewed always relative to other layers within a neural network. In this work,
we will explore how the absolute saturation value changes in different
scenarios. We further explore how the saturation evolves during training and
we will explain for the low dimensionality of the tail pattern.

##  3 Capacity and Problem Difficulty behave proportionally

In this section we analyze the relationship between problem difficulty and
model capacity by addressing two questions exploring how this relationship is
reflected in the saturation values. In order to do so, we train the entire
VGG-network family (VGG11, 13, 16 and, 19) on Cifar10 Krizhevsky (2009) and
reduce their capacity evenly over the entire architecture to observe how this
behavior affects saturation values. Our first hypothesis states that the
average saturation $s_{\mu}$ increases proportionally with a reduction in
capacity and a decrease in performance. We then move on to investigate how the
problem difficulty of the network changes the saturation emerging in a neural
architecture. Since the relevant eigenspace is generally larger when the layer
is contributing to the quality of the solution Richter et al. (2020), we
further hypothesize that more processing in a layer requires a larger relevant
eigenspace. If this assumption holds true, the overall saturation level should
increase with an increase in the difficulty of the task. If both working
hypotheses are true, we can conclude that the difficulty of the problem and
the capacity of the layers influence saturation in an antagonistic way.

###  3.1 Methodology

We test our working hypotheses by conducting two experiments. We first train
the VGG-family Simonyan and Zisserman (2015) of networks on Cifar10. We
further train 4 additional variants of each model, which have the respective
number of filters (and thus capacity) reduced by a factor of $\frac{1}{2}$,
$\frac{1}{4}$, $\frac{1}{8}$ and $\frac{1}{16}$. We choose Cifar10 for its
manageable size, which allows for a larger number of model training runs to be
conducted with our available resources, which is necessary for this
experiment. We choose the VGG-family of networks for its architectural
simplicity and because we can test different depths of convolutional neural
networks by experimenting on the entire family of networks. Training itself is
conducted using a stochastic gradient descent (SGD) optimizer with a learning
rate of 0.1, which is decaying after 10 epochs with a decay factor of 0.1. The
models are trained on a batch size of 64 for 30 epochs in total.

The second experiment is conducted on ResNet18. We are using a similar setup
to the first set of experiments. However, we are using a standardized input
resolution of $224\times 224$ pixels, to avoid artifacts caused by the input
resolution. We train the model on multiple datasets of different difficulties
(in ascending order of complexity): MNIST, Cifar10, TinyImageNet, and the
ImageNet dataset LeCun and Cortes (2010); Krizhevsky (2009); Le and Yang
(2015); Deng et al. (2009). While it is hard to precisely define the
complexity of the task, we think the selected datasets can be regarded as
increasingly difficult based on the number of classes and the complexity of
the visual information provided as data points to the model. MNIST features
strictly binary images of $28\times 28$ pixel native resolution for a 10-class
classification problem. Cifar10 features $32\times 32$ real world RGB images
belonging to 10 classes. TinyImageNet consists of $64\times 64$ images with
200 classes and ImageNet is made up of RGB images of various sizes belonging
to 1.000 classes.

###  3.2 Results

The average saturation $s_{\mu}$ is increasing and the predictive performance
is reduced when the capacity of the model is reduced. The exponential
reduction in capacity is reflected in a logarithmic relation between the
increasing $s_{\mu}$ and predictive accuracy measured on the test set (see
figure 2). From these observations, we can conclude that reducing the capacity
of the architecture leads to an increase in saturation.

Figure 2: Reducing the model capacity by reducing the number of filters
increases the average saturation and decreases performance (upper plot).
Training a model on more difficult datasets also increases the overall
saturation level. This indicates that saturation can measure the load on a
ResNet18 model (lower plot).

We further observe in figure 2 that saturation also increases with problem
complexity. The saturation levels of all layers increase when the model is
trained on a more difficult problem. The overall shape of the saturation curve
only deviates slightly, with no tail pattern or similar anomalous shapes
emerging. Since we know from the works of Richter et al. (2020) that a
resolution of $224\times 224$ pixels results in an even distribution of the
inference process for the trained model, we can conclude that less processing
is required for less complex problems.

Combined with the insights gained from training the VGG-variants on Cifar10,
we can conclude that for the pairs of dataset and model in our experiments a
saturation “sweet spot” exists between $s_{\mu}$ = 0.2 and $s_{\mu}=0.4$,
which yields the good predictive performance without being too excessively
overparameterized.

##  4 On the emergence of saturation patterns

The tail pattern that we discussed earlier in this work allows for the
identification of inefficiencies caused by mismatches between the neural
architecture and the input resolution. The analysis of saturation patterns is
generally done after training has concluded. However, since saturation can be
computed life during training with little overhead Richter et al. (2020), we
think that it might be interesting to see how these patterns emerge during the
training process.

###  4.1 Methodology

We first investigate how the saturation level evolves in a single layer under
different circumstances. We train a 3 layer fully connected neural network,
the first layer has 256 units, the second layers has 8, 16, 32, 64 and 128
units (a different number for each run). We train these networks using the
ADAM optimizer and a batch size of 128. The training is conducted twice. Once
using 8 epochs, which is enough for all models to converge, whereas the second
experiment is run for 20 epochs, which results in the loss increasing again
due to overfitting. We compute the saturation of the hidden layer of the
3-layer architecture after each epoch to observe the evolution of the
architecture. We further compute the loss of the architecture to observe a
potential relationship between loss and saturation convergence.

Based on these observations, we repeat the experiment on VGG11 and 19 as well
as sparse (low capacity) versions of these models with $\frac{1}{8}$ of the
original number of filters. We do this to observe whether the evolution of
saturation patterns is also depending on the architecture, depth, and capacity
of the network.

###  4.2 Results

In figure 3, we observe that an increase in the number of units in the fully
connected layer will result in an increased saturation during every epoch.
Furthermore, the saturation converges at a similar pace than the loss,
starting from a low-saturation level and increasing epoch by epoch.

(a) Validation loss during training.

(b) Saturation of layer 2 during training.

Figure 3: Saturation and loss converge at a similar pace.

This joint convergence is interrupted when the model starts to overfit. In
figure 4, we can see that the increase of validation loss is not followed by a
drastic change in saturation. Instead, the curve of saturation values flattens
for all sizes of the second layer. The fact that overfitting is not reflected
in saturation values indicates that the changes to the way the data is
processed when the model starts to overfit are subtle and thus are not
reflected in changes to the relevant eigenspace and therefore saturation.

(a) Validation loss during training.

(b) Saturation of layer 2 during training.

Figure 4: When overfitting the saturation pattern keeps converging, indicating
that overfitting is not drastically altering the relevant eigenspace.

In figure 5, we can see that the converging behavior observed with the fully
connected network can also be observed in a fully convolutional network. This
converging behavior is independent of the position of the layer in the
network, the number of layers, and the capacity of the network.

Figure 5: Saturations of convolutional neural networks show a converging
behavior regarding saturation similar to previous observations in figure 3.

Another interesting observation is that tail pattern seems to be observable
rather early during training, which indicates that an on-line analysis during
training allows the data scientist to detect inefficiencies early before
training has concluded.

##  5 Predictability of Tail Patterns Regarding Complexity

In the previous section, we show that the overall saturation increases when
the network’s capacity is low. In this section, we investigate how this
affects the predictability of tail patterns. Richter et al. (2021) show that
the tail patterns in sequential convolutional neural networks can be predicted
by computing the receptive field of all convolutional layers. The receptive
field size is the height and width of the area on the input image that can
influence the value of a single unit in a convolutional activation map. The
first layers ingesting an input with a receptive field size greater than the
input resolution are referred to as the “border layers” and mark the start of
the tail pattern. In section 4, we have seen that the global saturation level
can be altered by changing the number of filters in a convolutional layer.
Since technically the required capacity is no longer present in a single
layer, the network could distribute the processing among additional layers in
the architecture to process the data. However, if the receptive field
expansion is determining the number of unproductive layers (like the authors
of Richter et al. (2021) suggest), we will observe a tail pattern of
unproductive layers starting at the border layer, which is the same layer as
in the full capacity networks.

###  5.1 Methodology

We test the hypothesis by repeating the experiments conducted by Richter et
al. (2021) regarding the prediction of unproductive layers. The authors of
Richter et al. (2021) were able to predict unproductive layers by computing
the border layer for VGG11, 13, 16, and 19 on Cifar10. We reduce the capacity
of these models by reducing the filter size to $\frac{1}{8}$ of the original
size to see whether a drastic loss in capacity changes how the inference is
distributed. The models are trained for 30 epochs using the SGD-optimizer with
a learning rate of 0.1, decaying by a factor of 0.1 every 10 epochs. The batch
size is 64, each batch is channel-wise normalized, each image is randomly
cropped during inference time as well as randomly horizontally flipped with a
probability of 50%. The receptive field and the border layer are computed
using the formulas provided by Richter et al. (2021).

###  5.2 Results

(a) VGG11

(b) VGG13

(c) VGG16

(d) VGG19

Figure 6: Performance improvements of the logistic regression probes past the
border layer are miniscule, even though the capacity of each layer is reduced
to $\frac{1}{8}$ of the original capacity. This indicates that the networks
are unable to shift processing to otherwise unused layers even if the capacity
is limited. This is consistent with observation made by Richter et al. (2021).

Even though the capacity of the networks has been significantly reduced in
every layer, the networks do not spread the inference process among
significantly more layers.

Based on these results we conclude that the inference dynamics of the tested
networks did not change substantially by reducing their capacity. This means
that the capacity of layers primarily interacts with the difficulty of the
problem, while the presence and absence of tail patterns interact with the
receptive field as exemplified by Richter et al. (2021).

##  6 Conclusion

In this work, we explore the properties of saturation deeper and integrate
this knowledge into insights by the authors of Richter et al. (2020) and
Richter et al. (2021). We find that saturation is influenced in two ways.
First, while the saturation of layers relative to other layers in the network
is indicative of mismatches between input resolution and architecture, the
average saturation level provides insights into the interaction of the
network’s capacity and the complexity of the problem. This observation opens
another way to detect inefficiencies based on over- or under-saturation of the
network. Second, we find that these axes of analysis are independent, meaning
that low capacity networks do not distribute their inference significantly
different among layers than high capacity networks of the same architecture.
The observed saturation patterns converge similar to the loss of network and
can thus be detected early during training allowing for quicker design
iterations. In summary, saturation can be useful as a tool for both,
optimizing neural architectures for specific problems and for performing more
principled investigation of information processing in deep neural networks.

## References

  * Alain and Bengio (2018) Guillaume Alain and Yoshua Bengio.  Understanding intermediate layers using linear classifier probes, 2018\.
  * Richter et al. (2020) Mats L. Richter, Justin Shenk, Wolf Byttner, Anders Arpteg, and Mikael Huss.  Feature space saturation during training, 2020.
  * Tan et al. (2018) Mingxing Tan, Bo Chen, Ruoming Pang, Vijay Vasudevan, and Quoc V. Le.  Mnasnet: Platform-aware neural architecture search for mobile.  _CoRR_ , abs/1807.11626, 2018.  URL http://arxiv.org/abs/1807.11626.
  * Tan and Le (2019) Mingxing Tan and Quoc Le.  EfficientNet: Rethinking model scaling for convolutional neural networks.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pages 6105–6114. PMLR, 09–15 Jun 2019.  URL http://proceedings.mlr.press/v97/tan19a.html.
  * He et al. (2015) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  _CoRR_ , abs/1512.03385, 2015.  URL http://arxiv.org/abs/1512.03385.
  * Raghu et al. (2017) Maithra Raghu, Justin Gilmer, Jason Yosinski, and Jascha Sohl-Dickstein.  SVCCA: Singular vector canonical correlation analysis for deep learning dynamics and interpretability, 2017.
  * Real et al. (2019) Esteban Real, Alok Aggarwal, Yanping Huang, and Quoc V Le.  Regularized evolution for image classifier architecture search, 2019.
  * Huang et al. (2019) Yanping Huang, Youlong Cheng, Ankur Bapna, Orhan Firat, Mia Xu Chen, Dehao Chen, HyoukJoong Lee, Jiquan Ngiam, Quoc V. Le, Yonghui Wu, and Zhifeng Chen.  Gpipe: Efficient training of giant neural networks using pipeline parallelism, 2019.
  * Krizhevsky (2009) Alex Krizhevsky.  Learning multiple layers of features from tiny images.  Technical report, MIT and NYU, 2009.
  * Simonyan and Zisserman (2015) Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition, 2015\.
  * LeCun and Cortes (2010) Yann LeCun and Corinna Cortes.  MNIST handwritten digit database.  http://yann.lecun.com/exdb/mnist/, 2010.  URL http://yann.lecun.com/exdb/mnist/.
  * Le and Yang (2015) Ya Le and Xuan Yang.  Tiny ImageNet Visual Recognition Challenge, 2015.
  * Deng et al. (2009) J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei.  ImageNet: A Large-Scale Hierarchical Image Database.  In _CVPR09_ , 2009.
  * Richter et al. (2021) Mats L. Richter, Wolf Byttner, Ulf Krumnack, Ludwdig Schallner, and Justin Shenk.  Size matters, 2021.
  * Shenk et al. (2019) Mats L. Richter, Wolf Byttner, Ulf Krumnack, Ludwdig Schallner, and Justin Shenk.  Spectral Analysis of Latent Representations, 2021.