# Convergent least-squares optimisation methods for variational data
assimilation

C. Cartis  Mathematical Institute, University of Oxford, UK  M. H. Kaouri
Department of Mathematics and Statistics, University of Reading, UK
Corresponding author. Address: Department of Mathematics and Statistics,
University of Reading, PO Box 220, Reading, RG6 6AX, UK. Email:
<EMAIL_ADDRESS>A. S. Lawless  Department of Mathematics and Statistics,
University of Reading, UK Department of Meteorology, University of Reading, UK
National Centre for Earth Observation, Reading, UK  N. K. Nichols  Department
of Mathematics and Statistics, University of Reading, UK Department of
Meteorology, University of Reading, UK  National Centre for Earth Observation,
Reading, UK

###### Abstract

Data assimilation combines prior (or background) information with observations
to estimate the initial state of a dynamical system over a given time-window.
A common application is in numerical weather prediction where a previous
forecast and atmospheric observations are used to obtain the initial
conditions for a numerical weather forecast. In four-dimensional variational
data assimilation (4D-Var), the problem is formulated as a nonlinear least-
squares problem, usually solved using a variant of the classical Gauss-Newton
(GN) method. However, we show that GN may not converge if poorly initialised.
In particular, we show that this may occur when there is greater uncertainty
in the background information compared to the observations, or when a long
time-window is used in 4D-Var allowing more observations. The difficulties GN
encounters may lead to inaccurate initial state conditions for subsequent
forecasts. To overcome this, we apply two convergent GN variants (line search
and regularisation) to the long time-window 4D-Var problem and investigate the
cases where they locate a more accurate estimate compared to GN within a given
budget of computational time and cost. We show that these methods are able to
improve the estimate of the initial state, which may lead to a more accurate
forecast.

Keywords: Data assimilation, Gauss-Newton, least squares, line search,
optimisation, regularisation

Highlights:

  * •

Poor initialisation of Gauss-Newton method may result in failure to converge.

  * •

Safeguarded Gauss-Newton improves initial state estimate within limited
time/cost.

  * •

Results using twin experiments with long time-window and chaotic Lorenz
models.

  * •

Apply state of the art least-squares convergence theory to data assimilation.

  * •

Improvements to initial state estimate may lead to a more accurate forecast.

##  1 Introduction

Four-dimensional variational data assimilation (4D-Var) aims to solve a
nonlinear least-squares problem that minimizes the error in a prior estimate
of the initial state of a dynamical system together with the errors between
observations and numerical model estimates of the states of the system over
time. In Numerical Weather Prediction (NWP), 4D-Var is used to estimate the
initial conditions for a weather forecast [24]. The 4D-Var scheme is able to
incorporate information from a previous forecast along with observations over
both temporal and spatial domains, weighted by the uncertainties in the prior
and the observations. From a Bayesian point of view the solution is the
maximum a posteriori estimate of the initial state [35]. The nonlinear least-
squares objective function is minimized using an iterative method. The quality
of the estimate and the subsequent forecast depends on how accurately the
4D-Var problem is solved within the time and computational cost available.

In this paper, we investigate the application of globally convergent
optimisation methods to the 4D-Var problem; such methods use safeguards to
guarantee convergence from an arbitrary initial estimate by ensuring a
sufficient, monotonic/strict decrease in the objective function at each
iteration. We focus on the strong-constraint 4D-Var problem where we assume
that the numerical model of the system perfectly represents the true dynamics
of the system or the model errors are small enough to be neglected. This
results in the formulation of variational data assimilation as an
unconstrained nonlinear least-squares problem and is employed by many
operational meteorological centres [37], including the Meteorological Service
of Canada [14], the European Centre for Medium-range Weather Forecasting
(ECMWF) [8, 38] and the Met Office [39].

Ideally in large-scale unconstrained optimisation, we seek a fast rate of
convergence, which can be achieved in nondegenerate cases using a Newton-type
method. However, these methods require the use of second order derivatives of
the objective function, which are too costly to compute and store
operationally. Therefore, optimisation methods that approximate the high order
terms, such as limited memory Quasi-Newton [15, 27, 40, 44], Inexact Newton
[10], Truncated Newton [23, 42], Adjoint Newton [41], Hessian-free Newton [9],
Gauss-Newton [11, 36] and Approximate Gauss-Newton [18] methods have been
considered. To compute efficiently the first derivatives of the objective
function required by these techniques, the adjoint of the numerical model is
generally used [24]. More recently, optimisation methods that do not require
the first derivatives of the objective function are being examined to avoid
the development and maintenance costs associated with using the adjoint [17].
Alternative data assimilation techniques that use ensemble methods to
approximate the objective function gradients, rather than using the adjoint,
are also being investigated [2, 26].

The incremental 4D-Var technique, used commonly in operational centres,
approximately solves a sequence of linear least-squares problems and has been
shown to be equivalent to the Gauss-Newton (GN) method under standard
conditions [22]. In the GN (or incremental) method the linearized problem is
solved in an inner loop of the algorithm; the solution to the nonlinear
problem is then updated in an outer loop and the problem is re-linearized. The
accuracy with which the inner loop is solved is known to affect the
convergence of the outer loop [21, 22]. In our work, we focus on the
convergence of the outer loop, where the exact gradient is used (as is the
case when an adjoint is available) and we assume that the inner loop linear
least-squares problem is solved exactly. Furthermore, we use a variable
transformation usually applied in operational 4D-Var to precondition the
optimisation problem, see [3].

A general drawback of the GN method is that given a poor initialisation, it is
not guaranteed to converge to a solution, known as the ‘analysis’ state, of
the 4D-Var problem [11]. In NWP, the initial guess for the minimisation is
generally chosen to be the predicted initial state from a previous forecast,
known as the ‘prior’ or ‘background’ state. However, for some applications of
4D-Var this choice may not be a good enough estimate of the analysis. We show
that in such cases, the GN minimisation procedure may fail to converge. There
are three main strategies that safeguard GN and make it convergent from an
arbitrary initial guess: line-search, regularisation and trust-region [7, 36].
These can all be regarded as variants of the original Levenberg-Marquardt
algorithm for solving nonlinear least-squares problems [25, 30]. In this work
we investigate GN method with regularisation (REG) and compare its performance
to GN with backtracking Armijo line-search (LS) and GN alone, applied to the
preconditioned 4D-Var problem when there is limited computational time and
evaluations available, such as in NWP.

In previous work, the use of a line-search strategy in combination with a
Quasi-Newton approach was implemented in the ECMWF NWP system to solve the
4D-Var problem and was found to improve the minimisation of the objective
function [15, 38]. This method uses the Wolfe line-search conditions [43] to
safeguard the convergence. The Wolfe conditions require the use of additional
evaluations of the objective function and its gradient, however, which is
computationally costly. Here we instead use the Armijo condition [1], which
only requires additional evaluations of the objective function and not the
gradient. We pair GN with backtracking Armijo line-search and use a fixed
number of computational evaluations to guarantee a reduction in the outer loop
objective function (assuming the inner loop is solved to a high accuracy). We
compare this method to the GN method and to the GN method safeguarded by
quadratic regularisation (REG), using a simple, inexpensive updating strategy.

Using two test models within the 4D-Var framework, we show that where there is
more uncertainty in the background information compared to the observations,
the GN method may fail to converge, yet the convergent methods, LS and REG,
are able to improve the estimate of the analysis. Assimilation over long time
windows is of particular interest. We use accuracy profiles to show
numerically that in the long time-window case and in cases where there is
higher uncertainty in the background information versus the observations, the
globally convergent methods are able to solve more problems than GN in the
limited cost available. By ‘solve’ we mean satisfying a criterion requiring a
reduction in the objective function within a set number of evaluations. We
also show the effect that poor background information has on the quality of
the estimate obtained. We consider the case where the background information
is highly inaccurate compared to the observations and find that the
convergence of all three methods is improved when more observations are
included along the time-window. Finally, for the case where GN performs well,
we recommend further research into the parameter updating strategies used
within the globally convergent methods.

The structure of this paper is organised as follows. In Section 2 we outline
the strong-constraint 4D-Var problem as a nonlinear least-squares problem and
the GN method that is frequently used to solve it. In Section 3 we outline the
globally convergent methods used within this paper. In Section 4 we describe
the experimental design including the dynamical models used. In Section 5 we
present the numerical results obtained when applying GN and the globally
convergent methods to the 4D-Var problem with different features. Finally, we
conclude our findings in Section 6. In an appendix we detail the proofs of
convergence for the REG and LS methods.

##  2 Variational data assimilation

###  2.1 4D-Var: least-squares formulation

In four-dimensional variational data assimilation (4D-Var), the analysis
$\mathbf{x}_{0}^{a}\in\mathbb{R}^{n}$ is obtained by minimising a objective
function consisting of two terms: the background term and the observation
term, namely;

$\mathcal{J}(\mathbf{x}_{0})=\frac{1}{2}(\mathbf{x}_{0}-\mathbf{x}_{0}^{b})^{T}\mathbf{B}_{0}^{-1}(\mathbf{x}_{0}-\mathbf{x}_{0}^{b})+\frac{1}{2}\sum_{i=0}^{N}(\mathbf{y}_{i}-\mathcal{H}_{i}(\mathbf{x}_{i}))^{T}\mathbf{R}_{i}^{-1}(\mathbf{y}_{i}-\mathcal{H}_{i}(\mathbf{x}_{i})).$
(1)

The background term measures the difference between the initial state of the
system and the background state vector $\mathbf{x}_{0}^{b}\in\mathbb{R}^{n}$,
which contains prior information. The observation term measures the difference
between information from observations at times $t_{i}$ in the observation
vector $\mathbf{y}_{i}\in\mathbb{R}^{p_{i}}$ and the model state vector
$\mathbf{x}_{i}\in\mathbb{R}^{n}$ at the same time through use of the
observation operator
$\mathcal{H}_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{p_{i}}$ that maps from
the model state space to the observation space. Both terms are weighted by
their corresponding covariance matrices to represent the uncertainty in the
respective measures, the background error covariance matrix
$\mathbf{B}\in\mathbb{R}^{n\times n}$ and the observation error covariance
matrices at times $t_{i}$, $\mathbf{R}_{i}\in\mathbb{R}^{p_{i}\times p_{i}}$,
which are assumed to be symmetric positive definite. We note that observations
are distributed both in time and space and there are usually fewer
observations available than there are state variables so $p<n$, where
$p=\sum_{i=0}^{N}p_{i}$. The 4D-Var objective function (1) is subject to the
nonlinear dynamical model equations which contain the physics of the system

$\mathbf{x}_{i}=\mathcal{M}_{0,i}(\mathbf{x}_{0}),$ (2)

where the nonlinear model
$\mathcal{M}_{0,i}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ evolves the state
vector from the initial time point $t_{0}$ to the time point $t_{i}$.

We precondition the 4D-Var problem using a variable transform, which has been
shown to improve the conditioning of the variational optimisation problem [19,
20]. To be able to use the negative square root of $\mathbf{B}$ in our
variable transformation, we first require the assumption that the matrix
$\mathbf{B}$ is full rank. This assumption is satisfied for our choices of
$\mathbf{B}$ in Section 5. We define a new variable $\mathbf{v}$ to be,

$\mathbf{v}=\mathbf{B}^{-1/2}(\mathbf{x}_{0}-\mathbf{x}_{0}^{b}).$ (3)

The 4D-Var objective function can then be written in terms of $\mathbf{v}$,
known as the control variable in data assimilation (DA), and minimised with
respect to this instead. Furthermore, by including the model information
within the objective function, we are able to write the constrained
optimisation problem (1)-(2) in the form of an unconstrained optimisation
problem and apply the minimisation methods described later in this paper. The
preconditioned 4D-Var objective function is given by

$\mathcal{J}(\mathbf{v})=\frac{1}{2}\mathbf{v}^{T}\mathbf{v}+\frac{1}{2}\sum_{i=0}^{N}(\mathbf{y}_{i}-\mathcal{H}_{i}(\mathcal{M}_{0,i}(\mathbf{B}^{1/2}\mathbf{v}+\mathbf{x}_{0}^{b})))^{T}\mathbf{R}_{i}^{-1}(\mathbf{y}_{i}-\mathcal{H}_{i}(\mathcal{M}_{0,i}(\mathbf{B}^{1/2}\mathbf{v}+\mathbf{x}_{0}^{b}))).$
(4)

We note that the function (4) is continuously differentiable if the operators
$\mathcal{H}_{i}$ and $\mathcal{M}_{0,i}$ are continuously differentiable. To
save both computational cost and time in 4D-Var, tangent linear approximations
of the nonlinear operators in (4) are use in the inner loop [8]. The tangent
linear model (TLM) and tangent linear observation operator are usually derived
by linearising the discrete nonlinear model equations.

In a nonlinear least-squares problem, the function
$\mathcal{J}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ has a special form, as
defined in the following,

$\min_{\mathbf{v}}\mathcal{J}(\mathbf{v})=\frac{1}{2}\|\mathbf{r}(\mathbf{v})\|_{2}^{2},$
(5)

where $\mathbf{r}(\mathbf{v})=[r_{1}(\mathbf{v}),...,r_{n+p}(\mathbf{v})]^{T}$
and each $r_{j}:\mathbb{R}^{n}\rightarrow\mathbb{R}$, for $j=1,2,\dots,n+p$,
is referred to as a residual. In (5), $\|.\|_{2}$ denotes the $l_{2}$-norm,
which will be used throughout this paper. Equation (4) is equivalent to (5)
where the residual vector $\mathbf{r}(\mathbf{v})\in\mathbb{R}^{(n+p)}$ and
its Jacobian $\mathbf{J}(\mathbf{v})$ are given by

$\mathbf{r}(\mathbf{v})=\begin{pmatrix}\mathbf{v}\\\
\mathbf{R}_{0}^{-1/2}(\mathbf{y}_{0}-\mathcal{H}_{0}(\mathbf{B}^{1/2}\mathbf{v}+\mathbf{x}_{0}^{b}))\\\
\mathbf{R}_{1}^{-1/2}(\mathbf{y}_{1}-\mathcal{H}_{1}(\mathcal{M}_{0,1}(\mathbf{B}^{1/2}\mathbf{v}+\mathbf{x}_{0}^{b})))\\\
\vdots\\\
\mathbf{R}_{N}^{-1/2}(\mathbf{y}_{N}-\mathcal{H}_{N}(\mathcal{M}_{0,N}(\mathbf{B}^{1/2}\mathbf{v}+\mathbf{x}_{0}^{b})))\end{pmatrix}\text{
and }\mathbf{J}(\mathbf{v})=\begin{pmatrix}\mathbf{I}\\\
-\mathbf{R}_{0}^{-1/2}\mathbf{H}_{0}\mathbf{B}^{1/2}\\\
-\mathbf{R}_{1}^{-1/2}\mathbf{H}_{1}\mathbf{M}_{0,1}\mathbf{B}^{1/2}\\\
\vdots\\\
-\mathbf{R}_{N}^{-1/2}\mathbf{H}_{N}\mathbf{M}_{0,N}\mathbf{B}^{1/2}\\\
\end{pmatrix},$ (6)

where

$\mathbf{M}_{0,i}=\frac{\partial\mathcal{M}_{0,i}}{\partial\mathbf{v}}\big{|}_{\mathcal{M}_{0,i}(\mathbf{B}^{1/2}\mathbf{v}+\mathbf{x}_{0}^{b})}\text{
and
}\mathbf{H}_{i}=\frac{\partial\mathcal{H}_{0}}{\partial\mathbf{v}}\big{|}_{\mathcal{M}_{0,i}(\mathbf{B}^{1/2}\mathbf{v}+\mathbf{x}_{0}^{b})}$
(7)

are the Jacobian matrices of the model operator $\mathcal{M}_{0,i}$ and
observation operator $\mathcal{H}_{i}$ respectively, where
$\mathbf{M}_{0,i}\in\mathbb{R}^{n\times n}$ is the tangent linear of
$\mathcal{M}_{0,i}$ and $\mathbf{H}_{i}\in\mathbb{R}^{p_{i}\times n}$ is the
tangent linear of $\mathcal{H}_{i}$ [35]. In practice, an adjoint method is
used to calculate the gradient of (4), defined as

$\nabla\mathcal{J}(\mathbf{v})=\mathbf{J}(\mathbf{v})^{T}\mathbf{r}(\mathbf{v}).$
(8)

The Hessian is the matrix of second-order partial derivatives of (4),

$\nabla^{2}\mathcal{J}(\mathbf{v})=\mathbf{J}(\mathbf{v})^{T}\mathbf{J}(\mathbf{v})+\sum\limits_{j=1}^{n+p}r_{j}(\mathbf{v})\nabla^{2}r_{j}(\mathbf{v}).$
(9)

In data assimilation, the second-order terms in (9) are often difficult to
calculate in the time and cost available and too large to store, and so one
cannot easily use Newton-type methods for 4D-Var. Therefore, a first-order
approximation to the Hessian of the objective function (4) is used, resulting
in a GN method, and is given by

$\mathbf{S}=\mathbf{J}(\mathbf{v})^{T}\mathbf{J}(\mathbf{v})=\mathbf{I}+\sum_{i=0}^{N}\mathbf{B}^{1/2}\mathbf{M}_{0,i}^{T}\mathbf{H}_{i}^{T}\mathbf{R}_{i}^{-1}\mathbf{H}_{i}\mathbf{M}_{0,i}\mathbf{B}^{1/2},$
(10)

which is, by construction, full rank and symmetric positive definite. The
condition number in the $l_{2}$-norm of (10), $\kappa(\mathbf{S})$, is the
ratio of its largest and smallest eigenvalues and is related to the number of
iterations used for the linear minimisation problems in 4D-Var and how
sensitive the estimate of the initial state is to perturbations of the data.
We can use $\kappa(\mathbf{S})$ to indicate how quickly and accurately the
optimisation problem can be solved [16].

###  2.2 4D-Var implementation

The incremental 4D-Var method, which was first proposed for practical
implementation of the NWP problem in [8], has been shown to be equivalent to
the GN method when an exact TLM is used in the inner loop. When an approximate
TLM is used, the method is equivalent to an inexact GN method [18, 22]. A
summary of the GN method is given next.

Algorithm 2.1: GN algorithm applied to (5) [11].
Step $0$: Initialisation. Given $\mathbf{v}^{(0)}\in\mathbb{R}^{n}$ and some
stopping criteria. Set $k=0$. Step $1$: Check stopping criteria. While the
stopping criteria are not satisfied, do: Step $2$: Step computation. Compute a
step $\mathbf{s}^{(k)}$ that satisfies
$\mathbf{J}(\mathbf{v}^{(k)})^{T}\mathbf{J}(\mathbf{v}^{(k)})\mathbf{s}^{(k)}=-\mathbf{J}(\mathbf{v}^{(k)})^{T}\mathbf{r}(\mathbf{v}^{(k)}).$
(11) Step $3$: Iterate update. Set
$\mathbf{v}^{(k+1)}=\mathbf{v}^{(k)}+\mathbf{s}^{(k)}$, $k:=k+1$ and go to
Step 1.

In Algorithm 2.2, the updated control variable $\mathbf{v}^{(k+1)}$ is
computed by finding a step $\mathbf{s}^{(k)}$ that satisfies (11), which is
known as the preconditioned linearised subproblem. By substituting
$\mathbf{v}^{(k+1)}$ into (3) and rearranging, we obtain the current estimate
$\mathbf{x}_{0}^{(k+1)}$ of the initial state to the original nonlinear 4D-Var
problem.

To reduce the computational cost in large DA systems and to solve the DA
problem in real time, the series of problems (11) can be solved approximately
in the inner loop using iterative optimisation methods such as Conjugate
Gradient (CG) where a limited number of CG iterations are allowed and an exact
or approximate $\mathbf{J}$ is used [18]. We do not focus on this here and
assume that (11) is solved exactly.

We note that the step calculation (11) uniquely defines $\mathbf{s}^{(k)}$,
and $\mathbf{s}^{(k)}$ is a descent direction when $\mathbf{J}(\mathbf{v})$ is
full column rank. This is the case in 4D-Var as the Jacobian,
$\mathbf{J}(\mathbf{v})$ in (6) is full column rank due to the presence of the
identity matrix, thus ensuring that $\mathbf{s}^{(k)}$ is a descent direction.

The definitions of two solution types, namely, local and global minima, are
stated in Appendix A, along with a brief explanation of the local convergence
property of GN. Although the GN method benefits from local convergence
properties, convergence can only be guaranteed if the initial guess
$\mathbf{v}^{(0)}$ of the algorithm is in some neighbourhood around an
(unknown) local solution $\mathbf{v}^{*}$, that is, convergence from an
arbitrary initial guess is not guaranteed [11]. Even if the GN method does
converge, it may not necessarily converge to the global minimum due to the
fact that multiple local minima of a nonlinear least-squares objective
function may exist.

GN has no way of adjusting the length of the step $\mathbf{s}^{(k)}$ and
hence, may take steps that are too long and fail to decrease the objective
function value and thus to converge, see Example 10.2.5 in [11] and later in
Section 5 where the poor performance of GN is demonstrated. As GN only
guarantees local convergence, we are interested in investigating methods that
converge when $\mathbf{v}^{(0)}$ is far away from a local minimiser
$\mathbf{v}^{*}$. We refer to these methods as ‘globally convergent’.
Mathematical theory on global strategies can be found in [36] and [11]. Two
globally convergent methods are GN with line search and GN with quadratic
regularisation, which use a strategy within the GN framework to achieve
convergence to a stationary point given an arbitrary initial guess by
adjusting the length of the step. These methods will be presented in the next
section.

##  3 Globally convergent methods

Within this section, we outline the two globally convergent algorithms that we
apply in Section 5 to the preconditioned 4D-Var problem.

###  3.1 Gauss-Newton with line search (LS)

A line search method aims to restrict the step $\mathbf{s}^{(k)}$ in (11) so
as to guarantee a decrease in the value of $\mathcal{J}$. Within our work, an
inexact line search method known as the backtracking-Armijo (bArmijo)
algorithm is used within the inner loop of GN to find a step length $\alpha>0$
that satisfies the Armijo condition [1]. The Gauss-Newton with backtracking-
Armijo line search (LS) method is as follows.

Algorithm 3.1: LS algorithm applied to (5) [36].
Step $0$: Initialisation. Given $\mathbf{v}^{(0)}\in\mathbb{R}^{n}$,
$\tau\in(0,1)$ and $\beta\in(0,1)$ and $\alpha_{0}>0$ and some stopping
criteria. Set $k=0$. Step $1$: Check stopping criteria. While the stopping
criteria are not satisfied, do: Step $2$: Step computation. Compute a step
$\mathbf{s}^{(k)}$ that satisfies
$\mathbf{J}(\mathbf{v}^{(k)})^{T}\mathbf{J}(\mathbf{v}^{(k)})\mathbf{s}^{(k)}=-\mathbf{J}(\mathbf{v}^{(k)})^{T}\mathbf{r}(\mathbf{v}^{(k)})$
(12) and set $\alpha^{(k)}=\alpha_{0}$. Step $3$: Check Armijo condition.
While the following (Armijo) condition is not satisfied
$\mathcal{J}(\mathbf{v}^{(k)}+\alpha^{(k)}\mathbf{s}^{(k)})\leq\mathcal{J}(\mathbf{v}^{(k)})+\beta\alpha^{(k)}{\mathbf{s}^{(k)}}^{T}\nabla\mathcal{J}(\mathbf{v}^{(k)}),$
(13) do: Step $4$: Shrink stepsize. Set $\alpha^{(k)}:=\tau\alpha^{(k)}$ and
go to Step 3. Step $5$: Iterate update. Set
$\mathbf{v}^{(k+1)}=\mathbf{v}^{(k)}+\alpha^{(k)}\mathbf{s}^{(k)}$, $k:=k+1$
and go to Step 1.

In Algorithm 3.1, the control parameter $\beta$ in (13) is typically chosen to
be small (see [36]). The step equation (12) is the same as the GN step
equation (11); thus when $\alpha^{(k)}=1$, the GN and LS iterates coincide at
(the same) point $\mathbf{v}^{(k)}$. The use of condition (13) in this method
ensures that the accepted steps produce a sequence of strictly decreasing
function values given
$\triangledown\mathcal{J}(\mathbf{v}^{(k)})^{T}\mathbf{s}^{(k)}<0$. This
latter condition is satisfied by $\mathbf{s}^{(k)}$ defined in (12) whenever
$\mathcal{J}(\mathbf{v}^{(k)})$ is full column rank (which is the case here)
as mentioned in Section 2 [36].

Despite its global convergence property (see Appendix A.1), the LS method has
some disadvantages. We remark that the use of the step length $\alpha^{(k)}$
may sometimes unnecessarily shorten the step $\mathbf{s}^{(k)}$, slowing down
the convergence. Furthermore, LS may be computationally costly due to the need
to calculate the value of the function $\mathcal{J}$ each time $\alpha^{(k)}$
is adjusted, although more sophisticated updating strategies for $\alpha$ may
be used to try to reduce this effect.

Other line search strategies are possible such as Wolfe, Goldstein-Armijo and
more [36], but they are more involved and potentially more computationally
costly. As LS requires the re-evaluation of the outer loop objective function
each time it adjusts its line search parameter, its applicability to real
systems has been in doubt due to the computational cost limitations in 4D-Var
[38]. In Section 5, we show that given the same cost as the GN method, the LS
method can in some cases, better minimise the preconditioned 4D-Var objective
function.

###  3.2 Gauss-Newton with regularisation (REG)

The GN method may also be equipped with a globalisation strategy by including
a regularisation term $\gamma^{(k)}\mathbf{s}^{(k)}$ in the step calculation
(11) of Algorithm 2.2. This ensures that the accepted steps produce a sequence
of monotonically decreasing function values. This is a common variation of the
GN method known as the Levenberg-Marquardt method, proposed in [25] and [30].
The effect of the regularisation parameter $\gamma^{(k)}$ is to implicitly
control the length of the step $\mathbf{s}^{(k)}$. Increasing $\gamma^{(k)}$
shortens the steps, thus increasing the possibility that the procedure will
decrease the objective function in the next iteration. The REG method can be
summarised as follows.

Algorithm 3.2: REG algorithm applied to (5) [33].
Step $0$: Initialisation. Given $\mathbf{x}^{(0)}\in\mathbb{R}^{n}$,
$1>\eta_{2}\geq\eta_{1}>0$, $\gamma^{(0)}>0$ and some stopping criteria. Set
$k=0$. Step $1$: Check stopping criteria. While the stopping criteria are not
satisfied, do: Step $2$: Step computation. Compute a step $\mathbf{s}^{(k)}$
that satisfies
$\left(\mathbf{J}(\mathbf{v}^{(k)})^{T}\mathbf{J}(\mathbf{v}^{(k)})+\gamma^{(k)}\mathbf{I}\right)\mathbf{s}^{(k)}=-\mathbf{J}(\mathbf{v}^{(k)})^{T}\mathbf{r}(\mathbf{v}^{(k)}).$
(14) Step $3$: Iterate update. Compute the ratio
$\rho^{(k)}=\frac{\mathcal{J}(\mathbf{v}^{(k)})-\mathcal{J}(\mathbf{v}^{(k)}+\mathbf{s}^{(k)})}{\mathcal{J}(\mathbf{v}^{(k)})-m(\mathbf{s}^{(k)})},$
(15) where
$m(\mathbf{s}^{(k)})=\frac{1}{2}\|\mathbf{J}(\mathbf{v}^{(k)})\mathbf{s}^{(k)}+\mathbf{r}(\mathbf{v}^{(k)})\|_{2}^{2}+\frac{1}{2}\gamma^{(k)}\|\mathbf{s}^{(k)}\|_{2}^{2}.$
(16) Set
$\mathbf{v}^{(k+1)}=\begin{cases}\mathbf{v}^{(k)}+\mathbf{s}^{(k)},&\text{if
$\rho^{(k)}\geq\eta_{1}$}\\\ \mathbf{v}^{(k)},&\text{otherwise}.\end{cases}$
(17) Step $4$: Regularisation parameters update. Set
$\gamma^{(k+1)}=\begin{cases}\frac{1}{2}\gamma^{(k)},&\text{if
$\rho^{(k)}\geq\eta_{2}$}\text{ (very successful iteration)}\\\
\gamma^{(k)},&\text{if $\eta_{1}\leq\rho^{(k)}<\eta_{2}$}\text{ (successful
iteration)}\\\ 2\gamma^{(k)},&\text{otherwise,}\text{ (unsuccessful
iteration)}\end{cases}$ (18) Let $k:=k+1$ and go to Step 1.

As in Algorithms 2.2 and 3.1, the step equation (14) is solved directly in the
numerical experiments in Section 5. We note that when $\gamma^{(k)}=0$ in
(14), the REG step in (14) is the same as the GN step in (11). By comparing
(14) with (11), we are able to see how the REG step differs from the GN step.
The diagonal entries of the Hessian of the 4D-Var objective function (4) are
increased by the regularisation parameter $\gamma^{(k)}$ at each iteration of
the REG method. The method is able to vary its step between a GN and a
gradient descent step by adjusting $\gamma^{(k)}$ (see [36]) but may be costly
due to the need to calculate the value of the function $\mathcal{J}$ to assess
the step. Note that other choices of the factors $\frac{1}{2}$ and $2$ for
updating $\gamma^{(k)}$ in (18) are possible and even more sophisticated
variants for choosing $\gamma^{(k)}$ have been proposed. The proof of global
convergence of the REG method is presented in Appendix A.2.

##  4 Experimental design

Before evaluating the GN, LS and REG methods numerically, we first explain the
experimental design.

Twin experiments are commonly used to test DA methods. They use error
statistics that satisfy the DA assumptions as well as synthetic observations
generated by running the nonlinear model forward in time to produce a
reference state (not generally a local minimum of (5)). Within this section,
we define our choices for the twin experimental design. We begin by briefly
outlining two commonly used dynamical models, which are sensitive to initial
conditions (chaotic nature), a property shared with NWP models.

###  4.1 Models

Lorenz 1963 model (L63) Proposed in [28], the Lorenz 63 model (L63) is a
popular experimental dynamical system that represents meteorological processes
using a simple model. The model consists of three nonlinear, ordinary
differential equations given as

$\displaystyle\frac{dx}{dt}$ $\displaystyle=\sigma(y-x),$ (19)
$\displaystyle\frac{dy}{dt}$ $\displaystyle=x(\rho-z)-y,$
$\displaystyle\frac{dz}{dt}$ $\displaystyle=xy-\beta z,$

where the state vector consists of $n=3$ time-dependent variables
$\mathbf{x}=[x(t),y(t),z(t)]^{T}\in\mathbb{R}^{3}$. The scalar parameters are
chosen to be $\sigma=10$, $\rho=\frac{8}{3}$ and $\beta=28$, making the system
chaotic. A second-order Runge-Kutta method is used to discretise the model
equations using a time step $\Delta t=0.025$.

Lorenz 1996 model (L96) Another popular experimental system is the atmospheric
Lorenz 96 model (L96) [29] given by the following $n$ equations,

$\frac{dx_{j}}{dt}=-x_{j-2}x_{j-1}+x_{j-1}x_{j+1}-x_{j}+F,$ (20)

where $j=1,2,\ldots,n$ is a spatial coordinate. For a forcing term $F=8$ and
$n=40$ state variables, the system is chaotic [29]. The variables are evenly
distributed over a circle of latitude of the Earth with $n$ points with a
cyclic domain and a single time unit is equivalent to approximately 5
atmospheric days. A fourth-order Runge-Kutta method is used to discretise the
model equations using a time step $\Delta t=0.025$ (approximately 3 hours).

For both the L63 and L96 models, the time-window length $t_{a}$ is varied in
the numerical experiments in Section 5.1. We will now outline how we formulate
the twin experiments, beginning with generating the reference state.

###  4.2 Twin experiments

The reference state at time $t_{0}$, $\mathbf{x}^{ref}_{0}$ is used as the
basis of a twin experiment in the definition of the background state (the
initial guess for the optimisation algorithms) as well as to generate the
observations using a nonlinear model run called the ‘nature’ run. We begin by
explaining how we obtain $\mathbf{x}^{ref}_{0}$.

Reference state A vector of length $n$ is drawn from the uniform distribution
and used as the initial vector of state variables $\mathbf{x}^{rand}$. For the
L63 model, $\mathbf{x}^{rand}$ is integrated forward using a second-order
Runge-Kutta method, which is spun-up over 1000 time steps to obtain the
reference state on the model attractor for the L63 twin experiments,
$\mathbf{x}^{ref}_{0}\in\mathbb{R}^{3}$. This is the same for the L96 model
except a fourth-order Runge-Kutta method is used to obtain
$\mathbf{x}^{ref}_{0}\in\mathbb{R}^{40}$. The reference state at time $t_{0}$,
$\mathbf{x}^{ref}_{0}$ can then be used to obtain the full nonlinear model
trajectory.

We next explain how we obtain the background state vector used within our twin
experiments to be used as the initial guess for the optimisation algorithms.

Background In 4D-Var, the initial guess for the optimisation algorithm is
taken to be the background state at time $t_{0}$, $\mathbf{x}_{0}^{b}$, which
incorporates information from previous forecasts. In our experiments, the
background state vector $\mathbf{x}^{b}_{0}$ is generated by adding Gaussian
noise

$\mathbf{\varepsilon_{b}}\sim\mathcal{N}(0,\mathbf{B}),$ (21)

to the reference state at time $t_{0}$, $\mathbf{x}^{ref}_{0}$. For the
background error covariance matrix, we choose
$\mathbf{B}=\sigma_{b}^{2}\mathbf{I}_{n}$ where $\sigma_{b}^{2}$ is the
background error variance. The standard deviations of the errors from the
reference state for each model are based on the average order of magnitude of
the entries of $\mathbf{x}^{ref}_{0}$. For the L63 experiments,
$\sigma_{b}^{2}=0.25,1,6.25$ and $25$ represent a $5\%,10\%,25\%$ and $50\%$
error respectively. Similarly for the L96 experiments we set
$\sigma_{b}^{2}=0.0625,0.25,1.5625$ and $6.25$.

As previously mentioned, we generate synthetic observations from a nonlinear
model run using the reference state at time $t_{0}$, $\mathbf{x}^{ref}_{0}$.
We next describe the choices we made when specifying these observations.

Observations We consider both the spatial and temporal locations of the
observations. We assume that for both models observations of single state
variables are taken and $\mathbf{H}_{i}$ are the exact observation operators
at times $t_{i}$ used to map to observation space. For the L63 model, we
consider where we have $p=2$ observations, one of $x$ and one of $z$ per
observation location in time. For the L96 model, we consider where we have an
observation of the first half of the state variables per observation location
in time. This choice mimics what we may expect in reality where we have more
observations concentrated in one part of the globe. For both models, we first
consider where there is only one set of observations at time $N$ (Nobs1) and
then show the effect of using more observations along the time-window in later
experiments. We use imperfect observations where the observations
$\mathbf{y}_{i}$ are generated by adding Gaussian noise

$\mathbf{\varepsilon_{o}}\sim\mathcal{N}(0,\mathbf{R}_{i}),$ (22)

to $\mathbf{H}_{i}\mathbf{x}^{ref}_{i}$ for each observation location in time.
For the observation error covariance matrix we choose
$\mathbf{R}_{i}=\sigma_{o}^{2}\mathbf{I}_{p}$ where $\sigma_{o}^{2}$ is the
observation error variance. We expect the problem (4) to be more ill-
conditioned, thus difficult to solve accurately, when the ratio

$\frac{\sigma_{b}}{\sigma_{o}}$ (23)

is large [19, 20]. The ratio (23) controls the influence of the observation
term in the preconditioned objective function (4). For all experiments, we set
the standard deviation of the observation error to be $10\%$ of the average
order of magnitude of the entries of $\mathcal{H}(\mathbf{x}^{ref}_{i})$ for
both models. For the L63 model, this is $\sigma_{o}^{2}=1$ and for the L96
model, this is $\sigma_{o}^{2}=0.25$. We vary the background error variance
$\sigma_{b}^{2}$ above and below $\sigma_{o}^{2}$ such that the ratio (23)
varies. This can be thought of as having more confidence in the observations
compared to background when $\sigma_{b}>\sigma_{o}$ and vice versa.
Furthermore, as the initial guess is set to be the background state vector,
which is dependent on the value of $\sigma_{b}$, by varying $\sigma_{b}^{2}$
we are essentially varying the initial guess of the algorithms, thus
eliminating starting point bias from our results [4]. It is important to
recall here that under certain conditions, the GN method is known for its fast
convergence properties when in close vicinity to a local minimum, see [11]. By
choosing a small value of $\sigma_{b}^{2}$, we expect the performance of GN to
beat that of both LS and REG as it does not require the adjustment of the
additional parameters $\alpha^{(k)}$ and $\gamma^{(k)}$. Also, when assuming
that the observations are more accurate than the background, the use of more
observation locations in time means that we are constraining the estimate of
the initial state more tightly to the reference state in the twin experiment
design. The effect this has on the convergence of the optimisation methods
will be investigated. We next outline the algorithmic choices we have made.

###  4.3 Algorithmic choices

Stopping criteria We now outline the criteria used to terminate Algorithms
2.2, 3.1 and 3.2. Due to the limited time and computational cost available in
practice, the GN method is not necessarily run to convergence and a stopping
criterion is used to limit the number of iterations. Each calculation of the
residual vector $\mathbf{r}(\mathbf{v})$ requires the non-linear model to be
run forward to obtain the state at each observation location in time. This can
then be used to calculate the value of the objective function. Furthermore,
one run of the adjoint model is required to calculate the gradient.

To reduce computational cost in practical implementations of 4D-Var, the inner
loop problem is solved at a lower resolution than the outer loop problem [13].
However, as the dimension of the problems used within this paper are
relatively small compared to DA systems in practice, we solve the full
resolution inner loop problem using the full resolution residual and Jacobian
given in (6) and solve the inner loop problem using MATLAB’s backslash
operator where an appropriate solver is chosen according to the properties of
the Hessian matrix $\nabla^{2}\mathcal{J}(\mathbf{v})$ (see [31] for more
details). The limit on the total number of function and Jacobian evaluations
is achieved by using the following criterion

$k_{J}+l\leq\tau_{e},$ (24)

where $k_{J}$ is the total number of Jacobian evaluations (which is equivalent
to the number of outer iterations $k$ in 4D-Var), $l$ is the total number of
function evaluations and $\tau_{e}$ is the tolerance. The tolerance $\tau_{e}$
can be chosen according to the maximum number of evaluations desired. We note
that for GN, $k_{J}=l$ as the method requires as many Jacobian evaluations as
function evaluations. However, for both LS and REG there could be more than
one function evaluation per Jacobian evaluation since for unsuccessful steps,
the Jacobian is not updated so $k_{J}\leq l$.

To ensure that the algorithms are stopped before the function values stagnate,
the following common termination criterion based on the relative change in the
function at each iteration is also used

$\frac{|\mathcal{J}(\mathbf{v}^{(k-1)})-\mathcal{J}(\mathbf{v}^{(k)})|}{1+\mathcal{J}(\mathbf{v}^{(k)})}\leq\tau_{s},$
(25)

for $k\geq 1$, where $\tau_{s}$ is the tolerance, chosen to be $10^{-5}$. The
criterion (25) is used throughout Section 5 unless indicated otherwise.

We expect the norm of the gradient of the objective function,
$\|\nabla\mathcal{J}(\mathbf{v}^{(k)})\|$ to be close to zero at a stationary
point. The following termination criterion will be used in Section 5.2 to
identify whether or not a given method has located a stationary point

$\|\nabla\mathcal{J}(\mathbf{v}^{(k)})\|\leq 10^{-5}.$ (26)

Parameter choices For the LS method, we choose $\alpha_{0}=1$ so that the
first step assessed by the bArmijo rule is the GN step. We set $\beta=0.1$ and
to adjust the step length, $\tau=0.5$.

For the REG method, we select the initial regularisation parameter to be
$\gamma^{(0)}=1$ so that the condition in Algorithm 3.2, $\gamma^{(0)}>0$, is
satisfied and the REG step differs from the GN step. Furthermore, we choose
$\eta_{1}=0.1$ and $\eta_{2}=0.9$ to assess how well the model (16)
approximates the true function value at the next iteration.

For all three optimisation methods, we set $\tau_{e}=8,100$ or $1000$
depending on the experiment. The choice of $\tau_{e}=8$ comes from that which
is used operationally in the ECMWF Integrated Forecasting System [12], whereas
the choice of $\tau_{e}=100$ or $1000$ is used to measure the performance of
the optimisation methods when closer to convergence.

In order to best present our results, we use accuracy profiling described as
follows.

Accuracy profiles An accuracy profile shows the proportion of problems a given
method can solve within a fixed amount of work ($\tau_{e}$) and a given
tolerance ($\tau_{f}$) of the change in the function value [34]. To ensure the
robustness of our results, we apply the three optimisation methods to a series
of $n_{r}$ randomly generated problems, where the randomness occurs through
the background and observation error vectors, $\mathbf{\varepsilon_{b}}$ and
$\mathbf{\varepsilon_{o}}$. For each realisation, a new
$\mathbf{\varepsilon_{b}}$ and $\mathbf{\varepsilon_{o}}$ are generated from
their respective distributions, (21) and (22). The following criterion
proposed in [34] is used to flag that an estimate of the initial state has
been obtained by an optimisation method

$\frac{\mathcal{J}(\mathbf{v}_{0}^{(l)})-\mathcal{J}(\mathbf{v}_{0}^{t})}{\mathcal{J}(\mathbf{v}_{0}^{(0)})-\mathcal{J}(\mathbf{v}_{0}^{t})}\leq\tau_{f},$
(27)

where $\mathbf{v}_{0}^{t}$ is a solution of (4) referred to as the ‘truth’ and
$\tau_{f}$ is the tolerance. The measure (27) compares the optimality gap
$\mathcal{J}(\mathbf{v}_{0}^{(l)})-\mathcal{J}(\mathbf{v}_{0}^{t})$ relative
to the best reduction
$\mathcal{J}(\mathbf{v}_{0}^{(0)})-\mathcal{J}(\mathbf{v}_{0}^{t})$ [34]. This
ensures that the 4D-Var problem is only flagged as solved by the optimisation
method once the value of the objective function is within some error
($\tau_{f}$) of the truth.
For our problems, the truth is unknown. We only know that, due to the
nonlinearity of the 4D-Var problem, there may exist many values of
$\mathbf{v}_{0}$ that could minimise (4) locally. We are interested in the
estimate $\mathbf{v}_{0}^{t}$ that gives the greatest reduction in (4) that
any of the three methods can obtain. Therefore, we set the truth to be the
$\mathbf{v}_{0}^{(l)}$ obtained by any of the three methods that gives the
smallest function value within the given number of evaluations. Using this
criterion allows us to benchmark the methods against each other using accuracy
profiles.

For each experiment, we plot the proportion of the same $n_{r}=100$
realisations solved by each method against the relative accuracy obtained,
$\tau_{f}$. The relative accuracy obtained is varied using $\tau_{f}=10^{-i}$,
where $i=0,0.01,0.02,\dots,5$.

##  5 Numerical results

In this section, we present the results when applying GN, LS and REG using the
experimental design described in the previous section. We begin by conducting
experiments showing the effect of the length of the assimilation time-window
on the convergence of the three methods.

###  5.1 Effect of time-window length

We produce accuracy profiles for different time-window lengths to understand
the effect this has on the convergence of each method while limiting the
number of function and Jacobian evaluations to $\tau_{e}=8$. We choose a
background error of $50\%$ and an observation error of $10\%$ so that the
ratio (23) is large relative to the other cases we consider. For both the L63
and L96 models, we consider both short and long time-window lengths of 6 hours
($t_{a}=0.05$), 12 hours ($t_{a}=0.1$), 1 day ($t_{a}=0.2$) and 5 days
($t_{a}=1$) with the results shown in Figure 1.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 1: Accuracy profiles for the GN (black), LS (red) and REG (blue)
methods applied to the L63 and L96 problems using different time-window
lengths $t_{a}$. These show the proportion of $n_{r}=100$ problems solved by
each of the methods against the specified accuracy $-\log(\tau_{f})$ when
$\tau_{e}=8$. The GN line is below the LS line in (a), (b), (e), (f) and (g).

From Figure 1, we see that as the length of the time-window of both the L63
and L96 problems is increased, the performance of the GN, LS and REG methods
suffers.

For the L63 problems, Figures 1(a) and 1(b) show that GN and LS perform
similarly and solve more problems to the highest accuracy than REG. However,
as $\tau_{f}$ is increased, REG is solving all of the problems, so the REG
estimate must be close to that of GN and LS. In Figure 1(c), both LS and REG
solve fewer problems compared to GN, even for relatively large choices of
$\tau_{f}$. However, there is a choice of $\tau_{f}$ where all three methods
are solving all problems, again indicating that the LS and REG estimates are
close to the GN estimate. The initial guess for the three methods (the
background) appears to be close enough to the solution and so the GN step is
able to attain a sufficient decrease in the objective function as predicted by
its local convergence properties. LS and REG are inadvertently shortening the
GN step, which is a good step in the short time-window case. As we know, LS
and REG need to adjust their respective parameters, $\alpha^{(k)}$ and
$\gamma^{(k)}$ to attain GN’s fast local convergence, so LS and REG are
requiring more evaluations than GN to achieve the same result. For the L96
short time-window results in Figures 1(e), 1(f) and 1(g), this is not the
case. In fact, REG is outperforming GN and LS and it appears that LS is
mimicking the behaviour of GN quite closely as the GN step is attaining a
sufficient decrease in the objective function. However the decrease that the
REG step is achieving appears to be much greater for the L96 problems.
Therefore, REG is able to solve a greater number of problems within a higher
level of accuracy, which explains the difference between the L63 results in
Figures 1(a), 1(b) and the L96 results in 1(e) and 1(f).

The long time-window results for the L63 and L96 problems are shown in Figures
1(d) and 1(h), respectively. In both figures, LS is outperforming GN. For the
L63 problems, the performance of GN does not differ much from the performance
of REG. However, comparing the performance of GN in 1(c) with 1(d), we can see
that performance of GN has deteriorated greatly when increasing the length of
the time-window. In fact, in the results where even longer time-windows are
used (not included here), LS and REG outperform the GN method for the L63
problems, as in 1(h).

For the remainder of our experiments, we set $t_{a}=1$ in order to consider a
long time-window case only, as this is where we expect to see the greatest
benefit from the globally convergent methods.

###  5.2 Behaviour of methods and stagnation of GN

In order to gain an understanding of how the globally convergent methods, LS
and REG, contend with GN, we next demonstrate the behaviour of GN, LS and REG
when applied to typical preconditioned 4D-Var L63 and L96 problems, where the
background error is large and the time-window length is long.

Figure 2 shows the convergence plots for two typical realisations when using
the GN, LS and REG methods to obtain a solution to the preconditioned 4D-Var
problem with the L63 and L96 models. In this figure, the total number of
function and Jacobian evaluations allowed is set to $\tau_{e}=100$ for both
the L63 and the L96 problems to see if any progress is made beyond the number
of evaluations allowed in practice. We recall that GN updates the gradient (8)
when the function (4) is updated, so there are as many function evaluations as
Jacobian evaluations. However, both LS and REG only update the Jacobian on
successful iterations when there is a reduction in the objective function.
Therefore, the total number of evaluations used by each of the methods could
consist of a different combination of function and Jacobian evaluations. As in
Section 5.1, we set the ratio (23) to be large. It is in this case that we are
able to best demonstrate the benefit of the globally convergent methods, LS
and REG. In Figure 2, we set $\tau_{s}=10^{-3}$ to ensure that the methods
stop before the function values stagnate. As Figure 2 includes function
evaluations for both successful and unsuccessful step calculations, it is
natural to see jumps in the function values of LS and REG while their
parameters, $\alpha^{(k)}$ and $\gamma^{(k)}$ are being adjusted to guarantee
a reduction in the function.

(a)

(b)

Figure 2: Convergence plots showing the value of the objective function at
each iteration (including unsuccessful iterations) of the GN (black), LS (red)
and REG (blue) methods when applied to a L63 problem (a) and a L96 problem
(b).

For the L63 problems (Figure 2(a)), all three methods stop when the relative
change in the function criterion (25) is satisfied and before the limit on the
total number of function and Jacobian evaluations (24) is met. Table 1
supports this figure by showing the algorithmic output for each of the GN, LS
and REG methods when two different stopping criteria are used. From these
results, we see that both LS and REG stop at the same function value, although
REG requires fewer evaluations to do so, and that GN is converging towards a
larger value of the objective function (4) than LS and REG. By instead
stopping on the criterion (26) and setting $\tau_{e}=1000$, we see in Table 1
that all three methods are still making progress on the gradient and iterate
level, indicating that the methods are in fact locating stationary points
despite a small change in the function value beyond those shown in Figure 2.

Table 1: Table of algorithmic output when applying, GN, LS and REG to a typical realisation of the L63 problems, corresponding to Figure 2(a). Criteria | Method | $l$ | $k_{J}$ | $\mathcal{J}(\mathbf{v}^{(k_{J})})$ | $\|\mathbf{v}^{(k_{J})}-\mathbf{v}^{(k_{J}-1)}\|$ | $\|\nabla\mathcal{J}(\mathbf{v}^{(k_{J})})\|$
---|---|---|---|---|---|---
| GN | 20 | 20 | 81.55 | 0.42 | 86.35
(25) | LS | 27 | 14 | 8.69 | 0.03 | 5.18
| REG | 14 | 14 | 8.69 | 0.05 | 1.00
| GN | 101 | 101 | 78.87 | $3.54^{-8}$ | $8.47^{-6}$
(26) | LS | 43 | 27 | 8.69 | $8.21^{-7}$ | $8.31^{-6}$
| REG | 66 | 66 | 8.69 | $7.34^{-7}$ | $9.24^{-6}$

For the L96 problems (Figure 2(b)), LS and REG stop when (25) is satisfied and
before (24) is satisfied, whereas GN only satisfies (24). Table 2 supports
this figure by showing the algorithmic output for each of the GN, LS and REG
methods when two different stopping criteria are used. From these results, we
see that both GN and LS are stopping at a larger value of the objective
function (4) than REG. Recall that the norm of the gradient criterion (26) can
be used to identify whether or not a given method has located a stationary
point. The values of $\|\nabla\mathcal{J}(\mathbf{v}^{(k_{J})})\|$ for LS and
REG when the relative change in the function criterion (25) is used are much
smaller than that of GN. However, when we instead use the norm of the gradient
criterion (26) and limit the number of iterations to $\tau_{e}=1000$, the
methods stop on the limit of the number of iterations. Therefore, our results
do not indicate that the estimates of LS and REG may indeed be stationary
points of the objective function as they did for the L63 problems. However, LS
and REG are are able to make some improvement (REG more so than LS) on the
gradient norm level, unlike GN, which appears to fluctuate at gradient level,
even after $\tau_{e}=1000$ evaluations.

Table 2: Table of algorithmic output when applying, GN, LS and REG to a typical realisation of the L96 problems, corresponding to Figure 2(b). Criteria | Method | $l$ | $k_{J}$ | $\mathcal{J}(\mathbf{x}^{(k_{J})})$ | $\|\mathbf{v}^{(k_{J})}-\mathbf{v}^{(k_{J}-1)}\|$ | $\|\nabla\mathcal{J}(\mathbf{v}^{(k_{J})})\|$
---|---|---|---|---|---|---
| GN | 50 | 50 | 1728.99 | 20.02 | 5758.47
(25) | LS | 24 | 14 | 12.72 | 0.07 | 10.09
| REG | 19 | 16 | 5.52 | 0.08 | 1.89
| GN | 500 | 500 | 960.32 | 15.88 | 8015.13
(26) | LS | 967 | 32 | 12.71 | 0 | 10.09
| REG | 967 | 32 | 5.51 | 0 | 0.03

Table 2 shows that as LS and REG iterate beyond what is shown in Figure 2(b),
there is very little change in the value of the cost function, despite making
some change on the iterate and/or gradient level. The effect of rounding error
means that although we see progress made, the function value may remain
stagnant because of limitations in computer precision and because of the
conditioning of the problem. The condition number of the Hessian
$\kappa(\mathbf{S})$ can be used to indicate the accuracy we could be able to
achieve. In our work, both the L63 and L96 problems are well-conditioned.

The observed behaviour in this section is partly due to the fact that there is
no mechanism in GN to force it to converge as there is in LS and REG. The
benefit of these mechanisms is clearly shown in Figure 2(b) where the GN
method is stagnating while the LS and REG methods are converging, further
motivating our investigation of these methods.

###  5.3 Effect of background error variance

In this section, we study the effect on the performance of the three methods
when the uncertainty in the background information is increased whilst the
uncertainty in the observations is fixed. Figure 3 shows the accuracy profiles
used to benchmark the performance of the GN, LS and REG methods as the
tolerance $\tau_{f}$ is reduced, where $\tau_{e}=8$, while Figure 4 allows
$\tau_{e}$ to increase for both models with the increase chosen relative to
the dimension of the models, i.e. a larger increase in $\tau_{e}$ is allowed
for the L63 problems, where $n=3$, than the L96 problems, where $n=40$. From
both these figures, we generally see that as the error in the background is
reduced, the performance of all three methods improves. The conditioning of
the problem has been shown to depend on the ratio of the standard deviations
of the background and observation errors (23) [19, 20]. Therefore, the
estimate obtained by any of the optimisation methods may not be accurate
enough to produce a reliable forecast if the ratio (23) is large. The accuracy
of the estimate obtained by each method will be investigated further later on
in the paper.

Figures 3(a) and 3(e) show that a globally convergent method is able to find a
smaller function value than GN. As the ratio (23) is reduced, from Figures
3(b), 3(c), 3(f) and 3(g) we see that the REG method is able to solve the most
problems at the highest level of accuracy. When there is less uncertainty in
the background versus the observations, Figure 3(d) shows that for the L63
problems, all three methods are solving close to all of the problems within a
high level of accuracy. This is because the three methods are able to solve a
large portion of the cases when the problem is well-conditioned, which could
explain this result. However, for the L96 problems Figure 3(h) shows that the
GN and LS methods are solving the majority of the problems and REG is not
performing as well at higher levels of background accuracy. We can see the
performance of REG improving for the L96 problems when more evaluations are
allowed in Figure 4(h).

In Figure 4, where more evaluations are allowed than in Figure 3, we see a
much greater difference between the globally convergent methods and GN when
the background error is larger than the observation error. In Figures 4(a),
4(b), 4(e) and 4(f), it appears that when more evaluations are allowed, the
performance of GN worsens relative to LS and REG in the case when $\sigma_{b}$
is large. The globally convergent methods are able to locate estimates of the
initial states for the preconditioned 4D-Var problem, which when compared to
GN, better minimise the objective function (4). When the background error is
the same as the observation error in Figure 4(c), it is GN that is performing
better than LS and REG for the L63 problems. For LS, this could be because LS
is unnecessarily shortening the GN step, causing slower convergence. For the
REG method, the regularisation parameter must be shrunk and therefore, REG
requires more iterations to benefit from GN’s fast convergence property.

In Figure 4(d), all three methods are solving essentially the same number of
problems, with a slight decrease in success for REG, that again could be due
to the need to adjust the regularisation parameter. For the L96 problems, we
see a slightly different result. Figures 4(g) and 4(h) show that a globally
convergent method is solving more problems, more accurately than GN despite
the background error being at most equal to the observation error. This is an
interesting result for this higher-dimensional model as we would expect GN to
locally converge at a faster rate than the globally convergent methods due to
the fact that GN does not need to adjust any parameters; however, we find this
not to be the case.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 3: Accuracy profiles for the GN (black), LS (red) and REG (blue)
methods applied to the L63 problems in (a)-(d) and the L96 problems in (e)-(h)
where $n_{r}=100$, $\tau_{e}=8$ and where there is one observation at the end
of the time-window. The observation error is $10\%$ and the background error
is varied above and below this, as indicated in the plot captions. The GN line
is below the LS line in (c), (d), (g) and (h).

In DA, we are interested in knowing the accuracy of the estimate obtained as
in applications such as NWP, the estimate is used as the initial conditions
for a forecast and so the quality of this forecast will depend on the errors
in the estimate. In the following section, we quantify and compare the errors
in the estimates obtained by each method.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 4: Accuracy profiles for the GN (black), LS (red) and REG (blue)
methods applied to the L63 problems where $\tau_{e}=1000$ in (a)-(d) and the
L96 problems where $\tau_{e}=100$ in (e)-(h). We set $n_{r}=100$ and there is
one observation at the end of the time-window. The observation error is $10\%$
and the background error is varied above and below this, as indicated in the
plot captions.

###  5.4 Quality of the analysis

We recall that the initial guess of the algorithms is the reference state
$\mathbf{x}_{0}^{ref}$ perturbed by the background error
$\mathbf{\varepsilon_{b}}$. In order to compare the quality of the estimate
obtained by each method, we compare their estimate to the reference state
$\mathbf{x}_{0}^{ref}$ to understand how far the estimates obtained by the
methods have deviated from this. The analysis error for each state variable is
given by $\varepsilon^{a}_{i}=x^{a}_{i}-x^{ref}_{i}$. For each realisation, we
calculate the root mean square error (RMSE) of the analysis error, which is
the difference between the reference state and the estimate obtained by each
method,

$RMSE=\frac{1}{\sqrt{n}}\|\varepsilon^{a}\|_{2}.$ (28)

For each method, we plot the percentage of problems solved (according to the
criterion (27) where $\tau_{f}=10^{-3}$) within a specified tolerance of the
RMSE (28). We acknowledge in this work that the code for the RMSE profiles has
been adapted from the code for the data profiles used in [34].

The results for the L63 and L96 problems are in Figure 5, which coincides with
the case shown in Figure 3 where $\tau_{f}=10^{-3}$. From this, we see that
the GN method solves fewer problems within the same level of RMSE accuracy as
LS and REG when the background error is large in Figures 5(a), 5(b), 5(e) and
5(f). Furthermore, we see how the RMSE of the analyses successfully found by
each method reduces as the background error variance is reduced. This can be
seen in the scale of the x axis in Figures 5(a), 5(b), 5(c) and 5(d) for the
L63 problems and Figures 5(e), 5(f), 5(g) and 5(h) for the L96 problems. For
both models, the concentration of points in Figures 5(a) and 5(e) shows us
that the LS method is solving more problems than GN and REG within the same
RMSE tolerance. A similar result can be seen for REG in Figures 5(b), 5(c),
5(f) and 5(g). In Figures 5(d) and 5(h), we see that all three methods are
performing similarly, the RMSE errors for each of the analyses are very close
together.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 5: RMSE plots for the GN (black), LS (red) and REG (blue) methods
applied to the L63 problems in (a)-(d) and the L96 problems in (e)-(h) where
$n_{r}=100$, $\tau_{e}=8$, $\tau_{f}=10^{-3}$ and where there is one
observation at the end of the time-window. The observation error is $10\%$ and
the background error is varied above and below this, as indicated in the plot
captions.

Including more observations constrains the solution to be closer to the
reference state when the observation error is small. We next show the effect
on the performance of the methods as we include more observations and see if
this gives any improvement in the performance of the methods when the
background error is much larger than the observation error.

###  5.5 Effect of observations

Within this section, we show how the use of more observation locations in time
affects the performance of the three methods. We take the worst case for the
three methods when there is a $50\%$ error in the background and see if
including more observations in time with a $10\%$ error affects the
performance of the methods. For both models, we consider only equally spaced
observations in time, one set of observations at time $N$ (Nobs1), times $N/2$
and $N$ (Nobs2), times $N/4,N/2,3N/4$ and $N$ (Nobs3) and the even time points
(Nobs4), where $N=40$. For the Nobs1 case, observations are based on the
reference state at the end of the time-window and more observations are
included over time in the Nobs2, Nobs3 and Nobs4 cases. This not only
increases the condition number of the problem but also constrains the estimate
more tightly to the reference state.

For the L63 problems from Figures 6(a), 6(b), 6(c) and 6(d), we see that as
the number of observation locations in time is increased, all three methods
are solving more problems at a higher level of accuracy. This is more apparent
when more evaluations are allowed as shown in Figure 7(a), 7(b), 7(c) and
7(d). Here, the performance of GN improves drastically between the Nobs1 and
Nobs2 cases (Figures 6(a) and 6(b)) while there is less significant change in
the behaviour of LS and REG. In Figure 6(d), we see that GN is able to solve
more problems than LS and REG. Again, this could be because the LS and REG
methods require more iterations to converge when GN is performing well due to
the need to adjust their parameters. This argument coincides with Figure 7(d)
where more evaluations are allowed and the LS and REG methods are able to
perform as well as or better than GN. For the L96 problems, we see a different
result. From Figure 6, we only see a significant improvement in the
performance of GN in the Nobs4 case (Figure 6(h)). Otherwise, there is little
effect. This conclusion can also be drawn from Figure 7(g) and 7(h) where more
evaluations are allowed.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 6: Accuracy profiles where $n_{r}=100$ and $\tau_{e}=8$ for the L63
problems in (a)-(d) and the L96 problems in (e)-(h) for different observation
locations in time, as indicated in the plot captions, where the background
error is $50\%$ and the observation error is $10\%$.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 7: Accuracy profiles where $n_{r}=100$ for the L63 problems where
$\tau_{e}=1000$ in (a)-(d) and the L96 problems where $\tau_{e}=100$ in
(e)-(h) for different observation locations in time, as indicated in the plot
captions, where the background error is $50\%$ and the observation error is
$10\%$. The GN line is below the LS line in (d).

Similar studies were carried out on the performance of GN, LS and REG when
applied to the preconditioned 4D-Var problem where we instead choose
$\mathbf{B}=\sigma_{b}^{2}\mathbf{C}_{B}$, where $\mathbf{C}_{B}$ is a
correlation matrix; similar conclusions are drawn but due to space
constraints, are not included within this paper.

##  6 Conclusion

We have shown that the globally convergent methods, LS and REG, have the
capacity to improve current estimates of the DA analysis within the limited
time and cost available in DA, through the use of safeguards within GN which
guarantee the convergence of the method from any initial guesses.

Using the L63 and L96 models in the preconditioned 4D-Var framework, we have
shown that when there is more uncertainty in the background information
compared to the observations, the GN method may fail to converge in the long
time-window case yet the globally convergent methods LS and REG are able to
improve the estimate of the initial state. We compare the quality of the
estimate obtained using the RMSE of the analysis and show that even in the
case where the background is highly inaccurate compared to the observations,
the globally convergent methods find estimates with an RMSE less than or equal
to the RMSE of the estimates GN obtains. We take the case where the background
is highly inaccurate compared to the observations and find that the
convergence of all three methods is improved when more observations are
included along the time-window. In addition to the numerical results, the
assumptions made in the global convergence theorems of both LS and REG when
applied to a general nonlinear least-squares problem and a discussion as to
whether these assumptions are satisfied in DA is presented in the appendix. We
note that preconditioning the second derivative matrix is not necessary for
these results to hold, although this is the case we have focused on within our
work.

Our findings are important in DA as they show that in cases where the accuracy
of the prior information is poor and when there is limited computational
budget, the globally convergent methods are able to minimise the 4D-Var
objective function, unlike GN. We recommend that these methods are tested on
DA problems with realistic models and for different applications to understand
if these conclusions continue to hold. In particular, one should consider such
problems where an accurate initial guess for the algorithms is unavailable and
a long assimilation time-window is used, as we found that it is in this case
that LS and REG have an advantage over GN.

Within this paper, the 4D-Var inner loop problem is solved exactly. In
practice this must be solved inexactly, due to the size of the control vector,
and by the use of approximations to meet the computational and time
constraints. This is a common area of research in the DA community in order to
improve the quality of the assimilation analysis as well as the speed of
convergence of the algorithms. Furthermore, in the case where GN performs
better than LS and REG, further research is needed on updating the
globalisation parameters (stepsize $\alpha^{(k)}$ and regularisation parameter
$\gamma^{(k)}$) to speed up convergence.

Acknowledgements This work has been funded in part by the UK Engineering and
Physical Sciences Research Council Centre for Doctoral Training in Mathematics
of Planet Earth, the University of Reading EPSRC studentship (part of
Grant/Award Number: EP/N509723/1) and by the NERC National Centre for Earth
Observation. We acknowledge in this work that the code for the Lorenz 1996
model was developed by Adam El-Said.

Declarations of interest None.

##  7 Appendix

##  Appendix A Convergence theorems

In this section, we outline some existing global convergence results for the
LS and REG methods and discuss whether the assumptions made hold in DA. We
first state the definitions of a local and global minimum of an optimisation
problem $\min_{\mathbf{v}\in\mathbb{R}^{n}}f(\mathbf{v})$ where
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and $\mathbf{v}\in\mathbb{R}^{n}$.

######  Definition A.1 (Local minimiser [36]).

A point $\mathbf{v}^{*}$ is a local minimiser of $f$ if there is a
neighbourhood $\mathcal{N}$ of $\mathbf{v}^{*}$ such that
$f(\mathbf{v}^{*})\leq f(\mathbf{v})$ for all $\mathbf{v}\in\mathcal{N}$.

######  Definition A.2 (Global minimiser [36]).

A point $\mathbf{v}^{*}$ is a global minimiser of
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ if $f(\mathbf{v}^{*})\leq
f(\mathbf{v})$ for all $\mathbf{v}\in\mathbb{R}^{n}$.

A global solution is difficult to locate in most cases due to the nonlinearity
of the problems. Therefore, a local solution is often sought by algorithms for
nonlinear optimisation.

We focus on nonlinear least-squares optimisation problems of the form (5) for
the remainder of this section. The GN method can only guarantee local
convergence under certain conditions and not necessarily global convergence.
This is dependent on how close the initial guess is from the local minimum the
algorithm locates and whether or not the residual vector $\mathbf{r}$ of (5)
is a zero vector at a solution $\mathbf{v}^{*}$. Furthermore, the region of
local convergence depends on problem constants not known a priori, such as
Lipschitz constants of the gradient.

A local convergence result for the GN method can be found in Theorem 10.2.1 of
[11] where the performance of GN is shown to be dependent on whether or not
the second-order terms in (9) evaluated at the solution $\mathbf{v}^{*}$ are
close to zero. Another local convergence result can be found in Theorem 4 of
[18] where GN is treated as an inexact Newton method. The theorem guarantees
convergence of the GN method if for each iteration $k=0,1,\ldots,$ the norm of
the ratio of $\mathbf{Q}(\mathbf{v}^{(k)})$ and
$\mathbf{J}(\mathbf{v}^{(k)})^{T}\mathbf{J}(\mathbf{v}^{(k)})$, the second and
first terms of (9) respectively, is less than or equal to some constant
$\hat{\eta}$ where $0\leq\hat{\eta}\leq 1$.

It is important to note here that the globally convergent methods we are
concerned with, namely LS and REG, can only guarantee global convergence to a
local minimum under certain conditions and not necessarily to a global
minimum.

Before we list the assumptions for the global convergence theorems, we first
state the definition of the Lipschitz continuity property of a general
function $g$ as this is widely used in the theorems.

######  Definition A.3 (Lipschitz continuous function (see [36] A.42)).

Let $g$ be a function where $g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ for
general $n$ and $m$. The function $g$ is said to be Lipschitz continuous on
some set $\mathcal{N}\subset\mathbb{R}^{n}$ if there exists a constant $L>0$
such that,

$\|g(\mathbf{v})-g(\mathbf{w})\|\leq
L\|\mathbf{v}-\mathbf{w}\|,\qquad\forall\mathbf{v},\mathbf{w}\in\mathcal{N}.$
(29)

The following assumptions are used to prove global convergence of both the LS
and REG methods.

######  A1.

$\mathbf{r}$ is uniformly bounded above by $\omega>0$ such that
$\|\mathbf{r}(\mathbf{v})\|\leq\omega$.

######  A2.

$\mathbf{r}\in\mathcal{C}^{1}(\mathbb{R}^{n})$ is Lipschitz continuous on
$\mathbb{R}^{n}$ with Lipschitz constant $L_{r}>0$.

######  A3.

$\mathbf{J}$ is Lipschitz continuous on $\mathbb{R}^{n}$ with Lipschitz
constant $L_{J}>0$.

We remark that for the LS method, we can weaken assumptions A2 and A3 using
the open set $\mathcal{N}$ containing the level set

$\mathcal{L}=\left\\{\mathbf{v}\in\mathbb{R}^{n}|\mathcal{J}(\mathbf{v})\leq\mathcal{J}(\mathbf{v}^{(0)})\right\\}.$
(30)

In order to achieve the sufficient decrease property of the LS method, the
following assumption must be made.

######  A4.

$\mathbf{J}(\mathbf{v})$ in (6) is uniformly full rank for all
$\mathbf{v}\in\mathbb{R}^{n}$, that is, the singular values of
$\mathbf{J}(\mathbf{v})$ are uniformly bounded away from zero, so there exists
a constant $\nu$ such that
$\|\mathbf{J}(\mathbf{v})\mathbf{z}\|\geq\nu\|\mathbf{z}\|$ for all
$\mathbf{v}$ in a neighbourhood $\mathcal{N}$ of the level set $\mathcal{L}$
where $\mathbf{z}\in\mathbb{R}^{n}$.

In 4D-Var practice, it is reasonable to assume that the physical quantities
are bounded. Therefore, we can say that both $\mathbf{x}_{0}-\mathbf{x}^{b}$
and the innovation vector $\mathbf{y}-\mathcal{H}(\mathbf{x})$ are bounded in
practice, thus satisfying assumption A1. In 4D-Var, we must assume that the
nonlinear model $\mathcal{M}_{0,i}$ is Lipschitz continuous in order for A2 to
hold. As discussed in [32], this is a common assumption made in the
meteorological applications. However, we cannot say that this is necessarily
the case in 4D-Var practice.

In order for the Jacobian $\mathbf{J}$ to be Lipschitz continuous, we require
its derivative to be bounded above by its Lipschitz constant. Therefore, for
assumption A3 to hold, we require $\mathbf{r}$ to be twice continuously
differentiable in practice, which is a common assumption made in 4D-Var, and
also, that these derivatives of $\mathbf{r}$ are bounded above.

As mentioned in Section 2, the preconditioned 4D-Var Hessian (10) is full rank
by construction as it consists of the identity matrix and a non-negative
definite term. Therefore, the Jacobian of the residual of the preconditioned
problem in (6) is full rank and assumption A4 holds. This is also the case for
the standard 4D-Var problem (1), because of the presence of $\mathbf{B}^{1/2}$
in its Jacobian.

We now outline the global convergence theorems for the LS and REG methods,
using these assumptions.

###  A.1 Global convergence of the LS method

Nocedal et al. outline the proof for the GN method with Wolfe line search
conditions in [36], which uses the Zoutendijk condition. This proof can be
adapted to prove the global convergence theorem of the LS method, Algorithm
3.1, given as follows.

######  Theorem A.4 (Global convergence for the Gauss-Newton with bArmijo line
search method, Algorithm 3.1).

Suppose we have a function $\mathcal{J}=\frac{1}{2}\mathbf{r}^{T}\mathbf{r}$
and its gradient $\nabla\mathcal{J}=\mathbf{J}^{T}\mathbf{r}$ where
$\mathbf{r}\in\mathcal{C}^{1}(\mathbb{R}^{n})$ and $\mathbf{J}$ is the
Jacobian of $\mathbf{r}$. Assume A1 \- A4 hold. Then if the iterates
$\\{\mathbf{v}^{(k)}\\}$ are generated by the GN method with stepsizes
$\alpha^{(k)}$ that satisfy the Armijo condition (13), we have

$\lim_{k\rightarrow\infty}\mathbf{J}(\mathbf{v}^{(k)})^{T}\mathbf{r}(\mathbf{v}^{(k)})=0.$
(31)

That is, the gradient norms converge to zero, and so the Gauss-Newton method
with bArmijo line search is globally convergent.

The proof of Theorem A.4 requires the bArmijo chosen stepsizes $\alpha^{(k)}$
to be bounded below, which can be derived using assumptions A1 \- A3. Using
this lower bound, as well as assumption A4, we are able to prove the
Zoutendijk condition (as in [36]) and its variant

$\sum_{k\geq
0}\cos(\theta^{(k)})\|\nabla\mathcal{J}(\mathbf{v}^{(k)})\|_{2}\|\mathbf{s}^{(k)}\|_{2}<\infty$
(32)

hold. Both the Zoutendijk condition and its variant (32) use the angle between
$\mathbf{s}^{(k)}$ (the GN search direction) and
$-\nabla\mathcal{J}(\mathbf{v}^{(k)})$ (the steepest descent direction),
$\theta^{(k)}$, which is given by

$\cos(\theta^{(k)})=\frac{(-\nabla\mathcal{J}(\mathbf{v}^{(k)}))^{T}\mathbf{s}^{(k)}}{\|\nabla\mathcal{J}(\mathbf{v}^{(k)})\|_{2}\|\mathbf{s}^{(k)}\|_{2}}.$
(33)

By showing that the angle is uniformly bounded away from zero with $k$, one
can show that GN with line search is a globally convergent method.

We will next present the global convergence theorem for the REG method. The
REG method has no sufficient decrease condition as in the LS method.
Therefore, the use of the level set (30) is not required. The assumptions for
convergence are similar to the LS method aside from the requirement of
$\mathbf{J}(\mathbf{v})$ being full rank.

###  A.2 Global convergence of the REG method

The global convergence theorem for the GN with quadratic regularisation
method, Algorithm 3.2, is given as follows.

######  Theorem A.5 (Global convergence for the Gauss-Newton with
regularisation method, Algorithm 3.2).

Suppose we have a function $\mathcal{J}=\frac{1}{2}\mathbf{r}^{T}\mathbf{r}$
and its gradient $\nabla\mathcal{J}=\mathbf{J}^{T}\mathbf{r}$ where
$\mathbf{r}\in\mathcal{C}^{1}(\mathbb{R}^{n})$ and $\mathbf{J}$ is the
Jacobian of $\mathbf{r}$. Assume A1 \- A3 hold. Then if the iterates
$\\{\mathbf{v}^{(k)}\\}$ are generated by the Gauss-Newton with regularisation
method, we have that

$\lim_{k\rightarrow\infty}\mathbf{J}(\mathbf{v}^{(k)})^{T}\mathbf{r}(\mathbf{v}^{(k)})=0.$
(34)

That is, the gradient norms converge to zero, and so the Gauss-Newton method
with regularisation is globally convergent.

We first note that some adaptations of the lemmas from the global convergence
proof of the Adaptive Regularisation algorithm using Cubics (ARC method) are
used to prove Theorem A.5, see [5] and [6]. We begin the proof by deriving an
expression for the predicted model decrease in terms of the gradient. We
require the use of an upper bound on $\gamma^{(k)}$, denoted as
$\gamma_{\max}$, which is derived using a property of Lipschitz continuous
gradients. We show that $\gamma^{(k)}\leq\gamma_{\max}$ for all $k\geq 0$ by
first showing that if $\gamma^{(k)}$ is large enough, then we have a
successful step so that $\gamma^{(k)}$ can stop increasing due to unsuccessful
steps in Algorithm 3.2. We use the expression for $\gamma_{\max}$ to prove
global convergence of the REG method under assumptions A1-A3 by showing that
the gradient norms converge to zero as we iterate.

Note that for both the LS and REG, if $\mathbf{r}(\mathbf{v}^{(k)})\rightarrow
0$, i.e. (5) is a zero residual problem, then we have that (31) and (34) hold
as $|\mathcal{J}(\mathbf{v}^{(k)})|$ is uniformly bounded. However, in
practice the variational problem is not usually a zero residual problem.

## References

  * [1] L. Armijo.  Minimization of functions having Lipschitz continuous first partial derivatives.  Pacific Journal of Mathematics, 16(1):1–3, 1966.
  * [2] R. Bannister.  A review of operational methods of variational and ensemble-variational data assimilation.  Quarterly Journal of the Royal Meteorological Society, 143(703):607–633, 2017.
  * [3] R. N. Bannister.  A review of forecast error covariance statistics in atmospheric variational data assimilation. II: Modelling the forecast error covariance statistics.  Quarterly Journal of the Royal Meteorological Society, 134(637):1971–1996, 2008.
  * [4] V. Beiranvand, W. Hare, and Y. Lucet.  Best practices for comparing optimization algorithms.  Optimization and Engineering, 18(4):815–848, 2017.
  * [5] C. Cartis, N. I. Gould, and P. L. Toint.  Adaptive cubic regularisation methods for unconstrained optimization. Part I: motivation, convergence and numerical results.  Mathematical Programming, 127(2):245–295, 2011.
  * [6] C. Cartis, N. I. Gould, and P. L. Toint.  Adaptive cubic regularisation methods for unconstrained optimization. Part II: worst-case function-and derivative-evaluation complexity.  Mathematical Programming, 130(2):295–319, 2011.
  * [7] A. R. Conn, N. I. Gould, and P. L. Toint.  Trust Region Methods.  SIAM, 2000.
  * [8] P. Courtier, J.-N. Thépaut, and A. Hollingsworth.  A strategy for operational implementation of 4D-Var, using an incremental approach.  Quarterly Journal of the Royal Meteorological Society, 120(519):1367–1387, 1994.
  * [9] D. N. Daescu and I. M. Navon.  An analysis of a hybrid optimization method for variational data assimilation.  International Journal of Computational Fluid Dynamics, 17(4):299–306, 2003.
  * [10] R. S. Dembo, S. C. Eisenstat, and T. Steihaug.  Inexact Newton methods.  SIAM Journal on Numerical Analysis, 19(2):400–408, 1982.
  * [11] J. E. Dennis Jr and R. B. Schnabel.  Numerical methods for unconstrained optimization and nonlinear equations.  SIAM, 1996.
  * [12] ECMWF.  ECMWF Newsletter No.158 Winter 2018/19.  (158):21–26, 2019.
  * [13] M. Fisher, J. Nocedal, Y. Trémolet, and S. J. Wright.  Data assimilation in weather forecasting: a case study in PDE-constrained optimization.  Optimization and Engineering, 10(3):409–426, 2009.
  * [14] P. Gauthier, M. Tanguay, S. Laroche, S. Pellerin, and J. Morneau.  Extension of 3DVAR to 4DVAR: Implementation of 4DVAR at the Meteorological Service of Canada.  Monthly Weather Review, 135(6):2339–2354, 2007.
  * [15] J. C. Gilbert and C. Lemaréchal.  Some numerical experiments with variable-storage quasi-Newton algorithms.  Mathematical Programming, 45(1-3):407–435, 1989.
  * [16] G. H. Golub and C. F. Van Loan.  Matrix computations, volume 3.  Johns Hopkins University Press, 2012.
  * [17] S. Gratton, P. Laloyaux, and A. Sartenaer.  Derivative-free optimization for large-scale nonlinear data assimilation problems.  Quarterly Journal of the Royal Meteorological Society, 140(680):943–957, 2014.
  * [18] S. Gratton, A. S. Lawless, and N. K. Nichols.  Approximate Gauss–Newton methods for nonlinear least squares problems.  SIAM Journal on Optimization, 18(1):106–132, 2007.
  * [19] S. A. Haben, A. S. Lawless, and N. K. Nichols.  Conditioning and preconditioning of the variational data assimilation problem.  Computers & Fluids, 46(1):252–256, 2011.
  * [20] S. A. Haben, A. S. Lawless, and N. K. Nichols.  Conditioning of incremental variational data assimilation, with application to the Met Office system.  Tellus A: Dynamic Meteorology and Oceanography, 63(4):782–792, 2011\.
  * [21] S. Laroche and P. Gauthier.  A validation of the incremental formulation of 4D variational data assimilation in a nonlinear barotropic flow.  Tellus A: Dynamic Meteorology and Oceanography, 50(5):557–572, 1998\.
  * [22] A. S. Lawless, S. Gratton, and N. K. Nichols.  An investigation of incremental 4D-Var using non-tangent linear models.  Quarterly Journal of the Royal Meteorological Society, 131(606):459–476, 2005.
  * [23] F.-X. Le Dimet, I. M. Navon, and D. N. Daescu.  Second-order information in data assimilation.  Monthly Weather Review, 130(3):629–648, 2002.
  * [24] F.-X. Le Dimet and O. Talagrand.  Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects.  Tellus A: Dynamic Meteorology and Oceanography, 38(2):97–110, 1986\.
  * [25] K. Levenberg.  A method for the solution of certain non-linear problems in least squares.  Quarterly of Applied Mathematics, 2(2):164–168, 1944.
  * [26] C. Liu, Q. Xiao, and B. Wang.  An ensemble-based four-dimensional variational data assimilation scheme. Part I: Technical formulation and preliminary test.  Monthly Weather Review, 136(9):3363–3373, 2008.
  * [27] D. C. Liu and J. Nocedal.  On the limited memory BFGS method for large scale optimization.  Mathematical Programming, 45(1-3):503–528, 1989.
  * [28] E. N. Lorenz.  Deterministic nonperiodic flow.  Journal of the Atmospheric Sciences, 20(2):130–141, 1963.
  * [29] E. N. Lorenz.  Predictability: A problem partly solved.  In Proc. Seminar on predictability, volume 1, 1996.
  * [30] D. W. Marquardt.  An algorithm for least-squares estimation of nonlinear parameters.  Journal of the Society for Industrial and Applied Mathematics, 11(2):431–441, 1963.
  * [31] MathWorks.  mldivide documentation.  https://www.mathworks.com/help/matlab/ref/mldivide.html, 2021.  [Online; accessed 30-May-2021].
  * [32] A. J. Moodey, A. S. Lawless, R. W. Potthast, and P. J. Van Leeuwen.  Nonlinear error dynamics for cycled data assimilation methods.  Inverse Problems, 29(2):025002, 2013.
  * [33] J. J. Moré.  The Levenberg-Marquardt algorithm: implementation and theory.  In Numerical analysis, pages 105–116. Springer, 1978.
  * [34] J. J. Moré and S. M. Wild.  Benchmarking derivative-free optimization algorithms.  SIAM Journal on Optimization, 20(1):172–191, 2009.
  * [35] N. K. Nichols.  Mathematical concepts of data assimilation.  In W. Lahoz, R. Swinbank, and B. Khattatov, editors, Data Assimilation: Making Sense of Observations, pages 13–39. Springer, 2010.
  * [36] J. Nocedal and S. J. Wright.  Numerical Optimization 2nd Edition.  Springer, 2006.
  * [37] F. Rabier.  Overview of global data assimilation developments in numerical weather-prediction centres.  Quarterly Journal of the Royal Meteorological Society, 131(613):3215–3233, 2005.
  * [38] F. Rabier, J.-N. Thépaut, and P. Courtier.  Extended assimilation and forecast experiments with a four-dimensional variational assimilation system.  Quarterly Journal of the Royal Meteorological Society, 124(550):1861–1887, 1998.
  * [39] F. Rawlins, S. Ballard, K. Bovis, A. Clayton, D. Li, G. Inverarity, A. Lorenc, and T. Payne.  The Met Office global four-dimensional variational data assimilation scheme.  Quarterly Journal of the Royal Meteorological Society, 133(623):347–362, 2007.
  * [40] D. F. Shanno and K.-H. Phua.  Remark on “Algorithm 500: Minimization of unconstrained multivariate functions [e4]”.  ACM Transactions on Mathematical Software (TOMS), 6(4):618–622, 1980.
  * [41] Z. Wang, K. Droegemeier, and L. White.  The adjoint Newton algorithm for large-scale unconstrained optimization in meteorology applications.  Computational Optimization and Applications, 10(3):283–320, 1998\.
  * [42] Z. Wang, I. Navon, X. Zou, and F. Le Dimet.  A truncated Newton optimization algorithm in meteorology applications with analytic Hessian/vector products.  Computational Optimization and Applications, 4(3):241–262, 1995\.
  * [43] P. Wolfe.  Convergence conditions for ascent methods.  SIAM Review, 11(2):226–235, 1969.
  * [44] X. Zou, I. M. Navon, M. Berger, K. H. Phua, T. Schlick, and F.-X. Le Dimet.  Numerical experience with limited-memory quasi-Newton and truncated Newton methods.  SIAM Journal on Optimization, 3(3):582–608, 1993.


# Overcoming Adversarial Attacks for Human-in-the-Loop Applications

Ryan McCoppin  Marla Kennedy  Platon Lukyanenko  Sean Kennedy

###### Abstract

Including human analysis has the potential to positively affect the robustness
of Deep Neural Networks and is relatively unexplored in the Adversarial
Machine Learning literature. Neural network visual explanation maps have been
shown to be prone to adversarial attacks. Further research is needed in order
to select robust visualizations of explanations for the image analyst to
evaluate a given model. These factors greatly impact Human-In-The-Loop (HITL)
evaluation tools due to their reliance on adversarial images, including
explanation maps and measurements of robustness. We believe models of human
visual attention may improve interpretability and robustness of human-machine
imagery analysis systems. Our challenge remains, how can HITL evaluation be
robust in this adversarial landscape?

Machine Learning, ICML, Adversarial Robustness, Human in the loop, active
learning, manipulated explanations, human vision, salience

##  1 Introduction

The adversarial advent has greatly impacted the machine learning landscape.
Adversarial images can degrade a model’s ability to perform its task. In
addition to degrading ML model performance, these perturbations are often
crafted to evade detection by image analysts. Additional data can be included
in image analysis applications to aid the analyst; including robustness
metrics and explanation maps. Even so, adversarial attacks have managed to
corrupt or evade many of these additional tools. HITL defenses are needed in
order to thwart attacks targeting the human element of computer vision. We
present the need for further research in HITL applications and human
observation in the adversarial domain.

##  2 Related Work

###  2.1 Attacking Active Learning

Humans are involved throughout the learning pipeline, including in curating
training data. Bespoke training data is unlabeled and is labeled at
substantial expense. Active learning involves more efficient labeling methods
and is largely separated from the field of adversarial machine learning. While
several studies have looked at misleading active learning querying (Miller et
al., 2017) there is much to explore in how active learning is affected by
adversarial attacks. Specifically, how can a model query and labeling be
enacted in an adversarial environment without poisoning the model.

###  2.2 Manipulating Model Explanations

Model utility hinges on user trust, which requires reasonable reliability
measures (e.g. confidence intervals) and understanding why it made a
particular decision. Past work has developed attacks that can modify
classifier explanations without changing the classifier prediction.
(Dombrowski et al., 2019). The relationship between model outputs and trust is
complicated, and a recent report found that user trust following attacks is
poorly studied (Xiong et al., 2022). Explanation maps which are expected to
reveal adversarial images (Ye et al., 2020) may be manipulated to disguise
adversarial attacks and model biases. This poses obstacles not only for model
trust but also for HITL evaluation (Figure 1). While solutions have been put
forward which provide robustness toward manipulation (Dombrowski et al.,
2022), the issue remains when analysts are looking at unfamiliar classes or
using traditional explanatory techniques.

Figure 1: Adversarial images can have manipulated explanations. Image A -
adversarial image; explanation reveals target class. Image B - adversarial
image; manipulated explanation hides target class. Manipulation based on
(Dombrowski et al., 2019)

##  3 Opportunity: Human-In-The-Loop Studies

###  3.1 Human Vision and Human Vision Models

Because humans are not fooled by adversarial images in the same way as deep
networks, humans and human vision models may contribute to adversarial
robustness. Human attention models can be used to predict where humans will
look when viewing a scene. Task-specific attention models can be built with
gaze tracking data or crowd-sourced using interfaces that direct users to
highlight (or deblur) areas of an image that are critical to their
classification decision (Linsley et al., 2017).

Disagreements between human attention models and model explanations may
indicate manipulated images, low salient classes or faulty models. Taking user
gaze into account could also improve model value by reducing user workload,
simplifying training, or improving user performance (Matzen et al., 2016).

###  3.2 Prototype Humans-in-the-Loop Tools

Active and Interactive machine learning methods allow users to label data
within an interface. This can be used to investigate attacks on
active/interactive labeling, user interfaces, and be used to train models from
user responses. One of our prototype HITL tools examines how users can
contribute to adversarial or poisoned image detection. Users assign ‘cards’
that contain images and metadata to ‘poisoned’ or ‘benign’ categories, as seen
in Figure 2. Visual explanations of a classifier’s decision, obtained via
Grad-CAM (Selvaraju et al., 2017), provide the user with regions of an image
the ML model focuses on during classification. As user annotated data is
collected, a poison detection ML model is updated via active learning. With
this tool, we aim to explore:

  * •

How do adversarial detection models compare to analyst detection capability?

  * •

What explanations do users find useful and convincing?

  * •

Do more robust explanations improve human performance?

  * •

Which attacks are really invisible? What tools can make them more visible?

Figure 2: A HITL labeling interface with explanation maps

##  4 Challenges

There are many challenges in HITL research concerning adversarial robustness
that have yet to find answers. These concerns include visualizing
perturbations, which explanation maps are robust and beneficial, the
reliability of active learning query methods and which additional metrics may
be beneficial to the analyst.

Generally, users are unable to detect adversarial images and thus rely on
”explanation maps” to provide visibility. These maps can also be adversarial
distorted (Dombrowski et al., 2019). If explanation maps are not always
reliable or truthful, what information can be provided to the user in order to
detect adversarial images and improve model robustness? In addition, is it
possible to gain a robustness to these distortions in a similar manner to
adversarial training? We believe human vision models may provide a key in
detecting and interpreting adversarial images. However, recent research has
shown that deep models of human attention are also vulnerable to adversarial
attack (Che et al., 2020). Can we improve the adversarial robustness of these
models and combine human and machine attention to identify adversarial images?
Research has already begun in adversarial robustness, but not in how it
relates to HITL evaluation and how analysts are able to use more robust
explanations.

## References

  * Che et al. (2020) Che, Z., Borji, A., Zhai, G., Ling, S., Li, J., and Le Callet, P.  A new ensemble adversarial attack powered by long-term gradient memories.  _AAAI Conference on Artificial Intelligence_ , 34:3405–3413, Apr. 2020.
  * Dombrowski et al. (2019) Dombrowski, A. K., Alber, M., Anders, C., Ackermann, M., Müller, K. R., and Kessel, P.  Explanations can be manipulated and geometry is to blame.  In _Advances in Neural Information Processing Systems 32_. Elsevier, 2019.
  * Dombrowski et al. (2022) Dombrowski, A. K., Anders, C. J., Müller, K. R., and Kessel, P.  Towards robust explanations for deep neural networks.  In _Pattern Recognition_. Elsevier, 2022.
  * Krizhevsky & Hinton (2009) Krizhevsky, A. and Hinton, G.  Learning multiple layers of features from tiny images.  Technical report, Computer Science Department, University of Toronto, Toronto, Canada, 2009.
  * Linsley et al. (2017) Linsley, D., Eberhardt, S., Sharma, T., Gupta, P., and Serre, T.  What are the visual features underlying human versus machine vision?  In _2017 IEEE International Conference on Computer Vision Workshops (ICCVW)_ , pp. 2706–2714, 2017.
  * Matzen et al. (2016) Matzen, L., Haass, M., Tran, J., and McNamara, L.  Using eye tracking metrics and visual saliency maps to assess image utility.  _Electronic Imaging_ , 2016:1–8, 02 2016.
  * Miller et al. (2017) Miller, D. J., Hu, X., Qiu, Z., and Kesidis, G.  Adversarial learning: A critical review and active learning study.  In _27th International Workshop on Machine Learning for Signal Processing_ , pp. 1–6. IEEE, 2017.
  * Russakovsky et al. (2015) Russakovsky, O., Deng, J., Su, H., Krause, J., Satheesh, S., Ma, S., Huang, Z., Karpathy, A., Khosla, A., Bernstein, M., Berg, A. C., and Fei-Fei, L.  Imagenet large scale visual recognition challenge.  In _International Journal of Computer Vision (IJCV)_. Springer, 2015\.
  * Selvaraju et al. (2017) Selvaraju, R. R., Cogswell, M., Das, A., Vedantam, R., Parikh, D., and Batra, D.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In _2017 IEEE International Conference on Computer Vision (ICCV)_ , pp. 618–626, 2017.
  * Xiong et al. (2022) Xiong, P., Buffett, S., Iqbal, S., Lamontagne, P., Mamun, M., and Molyneaux, H.  Towards a robust and trustworthy machine learning system development: An engineering perspective.  In _Journal of Information Security and Applications_ , pp. 65. Elsevier, 2022.
  * Ye et al. (2020) Ye, D., Chen, C., Liu, C., Wang, H., and Jiang, S.  Detection defense against adversarial attacks with saliency map.  In _International Journal of Intelligent Systems_. Wiley Online Library, 2020.