#  Convergence properties of a Gauss-Newton data-assimilation method

Nazanin Abedini<EMAIL_ADDRESS>Svetlana Dubinkina<EMAIL_ADDRESS>VU
Amsterdam, Department of Mathematics, De Boelelaan 1111, 1081 HV Amsterdam,
The Netherlands

###### Abstract

Four-dimensional weak-constraint variational data assimilation estimates a
state given partial noisy observations and dynamical model by minimizing a
cost function that takes into account both discrepancy between the state and
observations and model error over time. It can be formulated as a Gauss-Newton
iteration of an associated least-squares problem. In this paper, we introduce
a parameter in front of the observation mismatch and show analytically that
this parameter is crucial either for convergence to the true solution when
observations are noise-free or for boundness of the error when observations
are noisy with bounded observation noise. We also consider joint state-
parameter estimation. We illustrated theoretical results with numerical
experiments using the Lorenz 63 and Lorenz 96 models.

##  1 INTRODUCTION

Data assimilation (DA) estimates a state of a dynamical model given partial
noisy measurements (also called observations), e.g. [1]. A so-called
variational DA minimizes a cost function that is a difference between an
estimation and the observations provided that either the estimate is a
solution of the dynamical model (strong-constraint variational DA, see e.g.
[2, 3, 4, 5]) or the dynamical model equations are satisfied under model error
(weak-constraint variational DA, see e.g. [6, 7]).

Variational DA can be viewed as an approach of combining model and
observations based on Tikhonov-Philips regularization, which is used to solve
ill-posed inverse problems [8]. Inverse problems appear in diverse fields such
as geosciences and image reconstruction [9, 10, 11, 12]. Inverse problems are
concerned with seeking a (stationary) solution of a mathematical model given a
set of noisy and incomplete observations. Due to sparsity of observations, the
corresponding discrete inverse problem has a highly ill-conditioned
coefficient matrix. In order to obtain a stable solution to an ill-posed
inverse problem, regularization methods are required.

Motivated by the weak-constraint variational DA, we consider a DA method that
minimizes a cost function under assumption of model error. We are seeking a
solution over a time window at once akin to the four-dimensional variational
DA (WC4DVar). The main difference to WC4DVar is that the cost function we
consider has a parameter in front of the observation mismatch. As we will
show, this parameter plays a significant role in the convergence of the DA
method and has to be chosen carefully in order to achieve either convergence
or boundedness of the error. As in [13] we use the dynamical model to
regularize the least-squares minimization problem.

There are several research papers investigating error propagation of
variational data assimilation. The error is defined as a norm of the
difference between the DA estimate and the true trajectory, from which the
observations are generated. In [14], the authors considered the Lorenz 63
model and the Navier-Stocks equations and showed that the error converges to
zero asymptotically in time for noise-free observations. In [13, 15], it was
shown that for a nonlinear dynamical model the error is bounded under
assumption of contractivity in the high modes and bounded observation noise.
We do not assume contractivity but assume bounds on some operators. The other
difference in our analysis compared to the existing analysis of variational DA
is that we derive convergence and bounds on the error as iteration goes to
infinity and not as the time goes to infinity. The DA method considered in
this paper is based on a Gauss-Newton iteration of the least-squares
minimization problem, e.g. [16, 17], which is was also considered for
incremental four-dimensional DA [18] in [19, 20].

The paper is organised as follows. In Section 2, we describe the DA
minimization problem and introduce a Gauss-Newton DA iteration to solve it. In
Section 3, we derive an error convergence result for noise-free observations
and an error boundess result for noisy observations with bounded noise. In
Section 4, we extend the DA method to joint state-parameter estimation. In
Section 5, we illustrate the theoretical results with numercilac experiments
using the Lorenz 63 and Lorenz 96 models. Finally, we present our conclusions
in Section 6.

##  2 GAUSS-NEWTON DATA-ASSIMILATION METHOD

Let us consider the following nonlinear dynamical model

$\frac{du^{\dagger}}{dt}=f(u^{\dagger}),\quad
u^{\dagger}(t)\in\mathbb{R}^{n},\ t\in[0,T],$ (1)

where $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$. Since in many applications the
model is defined by the time-discretization, we consider data assimilation in
the context of a discrete deterministic model. Let
$0=t_{0}<t_{1}<\dots<t_{N-1}=T$ be an equidistant partition of $I=[0,T]$ with
$t_{j}=j\Delta t$, then the time-discretization of Equation (1) is

$u^{\dagger}_{j+1}=F_{j}(u^{\dagger}_{j}),\quad
u^{\dagger}_{j}\in\mathbb{R}^{n},\quad j=0,\dots,N-1,$ (2)

where $u^{\dagger}_{j+1}=u^{\dagger}(t_{j+1})$ and $F_{j}$ is twice
continuously differentiable $\forall j=0,\cdots,N$. Let us denote
$\mathbf{u}^{\dagger}=\\{u^{\dagger}_{0},\dots,u^{\dagger}_{N}\\}$ the true
solution of the model and assumed to be unknown. Suppose a sequence of noisy
observations $\mathbf{y}:=\\{y_{k_{0}},\dots,y_{k_{M}}\\}$, where $k_{0}\geq
0$ and $k_{M}\leq N$, related to $\mathbf{u}^{\dagger}$ is given as

$y_{j}=\mathcal{H}_{j}u^{\dagger}_{j}+\eta_{j},\quad
y_{j}\in\mathbb{R}^{b},\quad j=k_{0},\dots,k_{M},$ (3)

where $\mathcal{H}_{j}:\mathbb{R}^{n}\to\mathbb{R}^{b}$, $b\leq n$, is the
linear observation operator, and the observation noise $\eta_{j}$, either
determenistic or stochastic, is bounded. The goal of data assimilation is to
find a state $\mathbf{u}=\\{u_{0},\cdots,u_{N}\\}$ such that the distance
between the state $\mathbf{u}$ and the observation $\mathbf{y}$ is minimized.
The weak-constraint variational data assimilation (WC4DVar) [6, 3] minimizes
the distance to the observations $\mathbf{y}$ under the condition that the
estimate is a solution of the dynamical model (Equation (2)) under model
error. Namely, WC4DVar solves the following minimization problem

$\displaystyle\min_{\mathbf{u}\in\mathbb{R}^{nN}}\frac{1}{2}\\{\|G(\mathbf{u})\|^{2}+\|\mathbf{y}-H\mathbf{u}\|^{2}\\},$

where

$G(\mathbf{u})=\left(G_{0}(\mathbf{u})\ G_{1}(\mathbf{u})\ \cdots\
G_{N-1}(\mathbf{u})\right)^{T},\quad
G_{j}(\mathbf{u})=u_{j+1}-F_{j}(u_{j}),\quad j=0,\dots,N-1,$ (4)

and

$H=\left(\mathcal{H}_{k_{0}}\ \cdots\ \mathcal{H}_{k_{M}}\right)^{T}.$

In this paper, we abuse the notation of $L^{2}$-norm in $\mathbb{R}^{s}$ for
different values of dimension $s$. In WC4DVar the norms are typically weighted
by the error covariance matrices. We consider a similar minimization problem
with the sole difference of a parameter $\alpha$ in front of the observation
mismatch, namely

$\min_{\mathbf{u}\in\mathbb{R}^{nN}}\frac{1}{2}\\{\|G(\mathbf{u})\|^{2}+\alpha\|\mathbf{y}-H\mathbf{u}\|^{2}\\}.$
(5)

As we already mentioned in introduction, the parameter $\alpha$ is crucial for
the convergence and bound result on error.

We minimize the cost function (Equation (5)) as follows. We start with an
initial guess $\mathbf{u}^{(0)}=H^{T}\mathbf{y}+(I-H^{T}H)\mathbf{u}^{b}$,
where $\mathbf{y}$ is the observation and $\mathbf{u}^{b}$ is the background
solution, which is usually a forecast state from the previous analysis cycle.
Then the method proceeds by the iteration

$\mathbf{u}^{(k+1)}=\mathbf{u}^{(k)}-\left(G^{\prime}(\mathbf{u}^{(k)})^{T}G^{\prime}(\mathbf{u}^{(k)})+\alpha
H^{T}H\right)^{-1}\left(G^{\prime}(\mathbf{u}^{(k)})^{T}G(\mathbf{u}^{(k)})+\alpha
H^{T}(H\mathbf{u}^{(k)}-\mathbf{y})\right),$ (6)

where $k$ denotes the index of the Gauss-Newton’s iteration,
$G(\mathbf{u}^{(k)})$ defined in Equation (4), and
$G^{\prime}(\mathbf{u}^{(k)})$ is Jacobian of $G$ which has an $n(N-1)\times
nN$ block structure:

$G^{\prime}(\mathbf{u})=\begin{bmatrix}-F_{0}^{{}^{\prime}}(u_{0})&I&\\\
&-F_{1}^{{}^{\prime}}(u_{1})&I&\\\ &&\ddots&\ddots&\\\
&&&-F_{N-1}^{{}^{\prime}}(u_{N-1})&I\\\ \end{bmatrix}.$ (7)

##  3 ERROR ANALYSIS

We define the error between approximation $\mathbf{u}$ of the Gauss-Newton DA
method (Equation (6)) and the true solution $\mathbf{u}^{\dagger}$ by

$\mathbf{e}_{k}:=\mathbf{u}^{(k)}-\mathbf{u}^{\dagger},\quad\forall
k=0,1,\dots.$ (8)

We show that the norm of error (Equation (8)) either converges to zero or is
bounded, dependent whether the observations are noise-free or noisy.

###  3.1 Noise-free observations

We assume that the sequence of observations $\mathbf{y}$ is noise-free, thus
$\eta_{j}=0,\ j=0,\cdots,N$ in Equation (3). We show that the Gauss-Newton DA
method (Equation (6)) produces an accurate state estimation under the
vanishing noise assumption. In order to do this, we need the following
assumptions on $G^{\prime}$.

######  Assumption 3.1.

The Jacobian $G^{\prime}(\mathbf{u})$ defined in Equation (7) is globally
Lipschitz continuous with the Lipschitz constant denoted by $L_{1}>0$, namely

$\|G^{\prime}(\mathbf{u})-G^{\prime}(\mathbf{v})\|\leq
L_{1}\|\mathbf{u}-\mathbf{v}\|,\quad\forall\mathbf{u},\mathbf{v}\in\mathbb{R}^{nN},$
(9)

and there exists $\alpha>0$ such that

$\|(G^{\prime}(\mathbf{x})^{T}G^{\prime}(\mathbf{x})+\alpha
H^{T}H)^{-1}G^{\prime
T}(\mathbf{x})\|\leq(L_{1}c)^{-1},\quad\forall\mathbf{x}\in\mathbb{R}^{nN},$
(10)

with $L_{1}$ being the Lipschitz constant from Equation (9) and $c>0$ being an
upper bound on the norm of initial error $\|\mathbf{e}_{0}\|$.

We abuse the notation by denoting both vector and matrix norm by $\|\cdot\|$.

######  Theorem 1.

Let Assumption (3.1) hold and let observations (Equation (3)) to be noise-
free. Then norm of the error defined in Equation (8) converges to zero as
iteration goes to infinity.

###### Proof.

For the sake of simplicity, we use $G_{k}$ and $G^{{}^{\prime}}_{k}$ instead
of $G(\mathbf{u}^{(k)})$ and $G^{\prime}(\mathbf{u}^{(k)})$, respectively.
Substituting Equation (6) into Equation (8) we have

$\displaystyle\mathbf{e}_{k+1}$ $\displaystyle=$
$\displaystyle\mathbf{u}^{(k+1)}-\mathbf{u}^{\dagger}$ $\displaystyle=$
$\displaystyle\mathbf{u}^{(k)}-(G_{k}^{{}^{\prime}T}G^{{}^{\prime}}_{k}+\alpha
H^{T}H)^{-1}\left(G_{k}^{{}^{\prime}T}G_{k}-\alpha
H^{T}(\mathbf{y}-H\mathbf{u}^{(k)})\right)-\mathbf{u}^{\dagger}.$

By substituting Equation (3) with $\eta_{j}=0,\ \forall j=0,\cdots,N$ into the
expression above, we obtain the following

$\displaystyle\mathbf{e}_{k+1}$ $\displaystyle=$
$\displaystyle\mathbf{e}_{k}-(G_{k}^{{}^{\prime}T}G^{\prime}_{k}+\alpha
H^{T}H)^{-1}(G_{k}^{{}^{\prime}T}G_{k}+\alpha H^{T}H\mathbf{e}_{k}).$

Opening the brackets and using the following property

$\displaystyle I-(G_{k}^{{}^{\prime}T}G^{\prime}_{k}+\alpha H^{T}H)^{-1}\alpha
H^{T}H=(G_{k}^{{}^{\prime}T}G^{\prime}_{k}+\alpha
H^{T}H)^{-1}G_{k}^{{}^{\prime}T}G^{\prime}_{k},$

we get

$\displaystyle\mathbf{e}_{k+1}=(G_{k}^{{}^{\prime}T}G^{\prime}_{k}+\alpha
H^{T}H)^{-1}\left(G_{k}^{{}^{\prime}T}(G^{\prime}_{k}\mathbf{e}_{k}-G_{k})\right).$

Since $G(\mathbf{u}^{\dagger})=0$, we add it to the right-hand side and use
the mean value theorem to obtain the following:

$\displaystyle\mathbf{e}_{k+1}$
$\displaystyle=(G_{k}^{{}^{\prime}T}G^{\prime}_{k}+\alpha
H^{T}H)^{-1}(G_{k}^{{}^{\prime}T}G^{\prime}_{k}\mathbf{e}_{k}-G_{k}^{{}^{\prime}T}(G(\mathbf{u}^{(k)})-G(\mathbf{u}^{\dagger}))$
$\displaystyle=(G_{k}^{{}^{\prime}T}G^{\prime}_{k}+\alpha
H^{T}H)^{-1}G_{k}^{{}^{\prime}T}G^{\prime}_{k}\mathbf{e}_{k}$
$\displaystyle-(G_{k}^{{}^{\prime}T}G^{\prime}_{k}+\alpha
H^{T}H)^{-1}G_{k}^{{}^{\prime}T}\left(\int_{0}^{1}G^{\prime}(s\mathbf{u}^{(k)}+(1-s)\mathbf{u}^{\dagger})ds\right)\mathbf{e}_{k}$
$\displaystyle=(G_{k}^{{}^{\prime}T}G^{\prime}_{k}+\alpha
H^{T}H)^{-1}G_{k}^{{}^{\prime}T}\left(G^{\prime}_{k}-\int_{0}^{1}G^{\prime}(s\mathbf{u}^{(k)}+(1-s)\mathbf{u}^{\dagger})ds\right)\mathbf{e}_{k}$
$\displaystyle=(G_{k}^{{}^{\prime}T}G^{\prime}_{k}+\alpha
H^{T}H)^{-1}G_{k}^{{}^{\prime}T}\left(\int_{0}^{1}\left(G^{\prime}(\mathbf{u}^{(k)})-G^{\prime}\left(s\mathbf{u}^{(k)}+(1-s)\mathbf{u}^{\dagger}\right)\right)\mathbf{e}_{k}ds\right).$

Taking norm of both sides and using Lipschitz continuity (Equation (9)) on
$G^{\prime}$, we get

$\|\mathbf{e}_{k+1}\|\leq\frac{L_{1}}{2}\big{\lVert}(G_{k}^{{}^{\prime}T}G^{\prime}_{k}+\alpha
H^{T}H)^{-1}G_{k}^{{}^{\prime}T}\big{\rVert}\|\mathbf{e}_{k}\|^{2}.$

Using Equation (10), we conclude that

$\displaystyle\|\mathbf{e}_{k+1}\|\leq\frac{1}{2c}\|\mathbf{e}_{k}\|^{2}.$

Therefore, for $k=1,2,\dots$, we have

$\|\mathbf{e}_{k}\|\leq\left(\frac{1}{2c}\right)^{2^{k}-1}\|\mathbf{e}_{0}\|^{2^{k}}.$

Since $\|\mathbf{e}_{0}\|\leq c$, we get

$\displaystyle\|\mathbf{e}_{k}\|$ $\displaystyle\leq$
$\displaystyle\left(\frac{1}{2}\right)^{2^{k}-1}c,$

which leads to $\lim_{k\to\infty}\|\mathbf{e}_{k}\|=0$. ∎

###  3.2 Noisy observations

Next we consider noisy observations. We show that norm of the error defined in
Equation (8) is bounded. In order to prove this result, we require local
conditions on $G^{\prime}$ and bounded observation noise.

######  Assumption 3.2.

$G$ is continuously differentiable in the open convex set
$D\subset\mathbb{R}^{nN}$. The Jacobian $G^{\prime}(\mathbf{u})$ defined in
Equation (7) is locally Lipschitz continuous in $D$, with the Lipschitz
constant denoted by $L_{2}>0$, namely

$\|G^{\prime}(\mathbf{u})-G^{\prime}(\mathbf{v})\|\leq
L_{2}\|\mathbf{u}-\mathbf{v}\|,\quad\forall\mathbf{u},\mathbf{v}\in D,$ (11)

and there exists $0<\alpha<1$ such that

$\displaystyle\|(G^{\prime}(\mathbf{x})^{T}G^{\prime}(\mathbf{x})+\alpha
H^{T}H)^{-1}G^{\prime
T}(\mathbf{x})\|\leq(L_{2}c)^{-1},\quad\forall\mathbf{x}\in D,$ (12)
$\displaystyle\|H^{T}\eta\|\|(G^{\prime}(\mathbf{x})^{T}G^{\prime}(\mathbf{x})+\alpha
H^{T}H)^{-1}\|\leq c/2,\quad\forall\mathbf{x}\in D,$ (13)

with $L_{2}$ being the Lipschitz constant from Equation (11) and $c>0$ being
an upper bound on the norm of initial error $\|\mathbf{e}_{0}\|$.

######  Lemma 3.3.

(Dennis and Robert [21]) Let $G:\mathbb{R}^{l}\to\mathbb{R}^{m}$ be
continuously differentiable in the open convex set $D\subset\mathbb{R}^{l}$
and let the Jacobian of $G$ be Lipschitz continuous at $\mathbf{x}\in D$,
using a vector norm and the induced matrix operator norm and the Lipschitz
constant $L$. Then, for any $\mathbf{x+p}\in D$,

$\|G(\mathbf{x+p})-G(\mathbf{x})-G^{\prime}(\mathbf{x})\mathbf{p}\|\leq\frac{L}{2}\|\mathbf{p}\|^{2}.$

######  Theorem 2.

Let Assumption (3.2) hold and let observations (Equation (3)) to be noisy.
Then

$\limsup_{k\to\infty}\|\mathbf{e}_{k}\|\leq\frac{\alpha c}{1-\alpha},$ (14)

provided the iteration $\mathbf{u}^{(k)}$ does not leave the convex set $D$
for all $k$.

###### Proof.

For the sake of simplicity, we use $G_{k}$ and $G^{{}^{\prime}}_{k}$ instead
of $G(\mathbf{u}^{(k)})$ and $G^{\prime}(\mathbf{u}^{(k)})$, respectively.
Substituting Equation (6) into Equation (8) we have

$\displaystyle e_{k+1}$ $\displaystyle=$
$\displaystyle\mathbf{u}^{(k+1)}-\mathbf{u}^{\dagger}$ $\displaystyle=$
$\displaystyle\mathbf{u}^{(k)}-\mathbf{u}^{\dagger}-(G^{\prime
T}_{k}G^{\prime}_{k}+\alpha H^{T}H)^{-1}(G^{\prime T}_{k}G_{k}+\alpha
H^{T}(H\mathbf{u}^{(k)}-\mathbf{y})),$

Using Equation (3) in the above equation and adding and subtracting $G^{\prime
T}_{k}G^{\prime}_{k}\mathbf{e}_{k}$ leads to

$\displaystyle\mathbf{e}_{k+1}$ $\displaystyle=$
$\displaystyle\mathbf{e}_{k}-(G^{\prime T}_{k}G^{\prime}_{k}+\alpha
H^{T}H)^{-1}(G^{\prime T}_{k}G_{k}+\alpha H^{T}H\mathbf{e}_{k}-\alpha
H^{T}\eta))$ $\displaystyle=$ $\displaystyle\mathbf{e}_{k}-(G^{\prime
T}_{k}G^{\prime}_{k}+\alpha H^{T}H)^{-1}(G^{\prime
T}_{k}G^{\prime}_{k}\mathbf{e}_{k}+\alpha H^{T}H\mathbf{e}_{k}+G^{\prime
T}_{k}(G_{k}-G^{\prime}_{k}\mathbf{e}_{k})-\alpha H^{T}\eta))$
$\displaystyle=$ $\displaystyle-(G^{\prime T}_{k}G^{\prime}_{k}+\alpha
H^{T}H)^{-1}(G^{\prime
T}_{k}(G_{k}-G(\mathbf{u}^{\dagger})-G^{\prime}_{k}\mathbf{e}_{k})-\alpha
H^{T}\eta),$

since $G(\mathbf{u}^{\dagger})=0$. By taking norm of both sides of the above
equation and using Lemma 3.3, we get

$\|\mathbf{e}_{k+1}\|\leq\frac{L_{2}}{2}\|(G^{\prime
T}_{k}G^{\prime}_{k}+\alpha H^{T}H)^{-1}G^{\prime
T}_{k}\|\|\mathbf{e}_{k}\|^{2}+\alpha\|(G^{\prime T}_{k}G^{\prime}_{k}+\alpha
H^{T}H)^{-1}\|\|H^{T}\eta\|.$ (15)

Using Equations (12) and (13) in the above expression leads to

$\|\mathbf{e}_{k+1}\|\leq\frac{1}{2c}\|\mathbf{e}_{k}\|^{2}+\alpha\frac{c}{2},$

Using $(a+b)^{2}\leq 2a^{2}+2b^{2}$ in the inequality above gives rise to

$\|\mathbf{e}_{k}\|\leq
2^{-k}\left(\frac{1}{c}\right)^{2^{k}-1}\|\mathbf{e}_{0}\|^{2^{k}}+c\sum_{i=0}^{k-1}\left(\frac{\alpha^{2^{i}}}{2^{i+1}}\right),\quad\mbox{for}\quad
k=1,2,\dots.$

Since $\|\mathbf{e}_{0}\|\leq c$ and $2^{-i}<1$ for $\forall i=0,1,\dots$, we
get

$\|\mathbf{e}_{k}\|<2^{-k}c+c\sum_{i=0}^{k-1}\alpha^{2^{i}}.$

Since $0<\alpha<1$, we have
$\sum_{i=0}^{k-1}\alpha^{2^{i}}<\alpha\sum_{i=0}^{k}\alpha^{i}$ leading
consequently to

$\|\mathbf{e}_{k}\|<2^{-k}c+c\alpha\sum_{i=0}^{k}\alpha^{i},$

and in the limit of $k$ going to infinity

$\limsup_{k\to\infty}\|\mathbf{e}_{k}\|<\frac{\alpha c}{1-\alpha}.$

∎

######  Remark 3.4.

Note that in case of $\eta_{j}=0,\forall j=1,\cdots,N$, we have
$\|H^{T}\eta\|=0$ and Equation (13) is trivially satisfied, which in turn
implies assumption of $\alpha>0$ instead of $1>\alpha>0$. Furthermore,
Equations (11) and (12) are equivalent to Equations (9) and (10) but locally
and thus Theorem (2) is equivalent to Theorem (1) but locally.

##  4 JOINT STATE-PARAMETER ESTIMATION

Data assimilation can also be used if the dynamical model depends on uncertain
parameter. We extend the Gauss-Newton DA method Equation (6) to joint state-
parameter estimation. We consider

$G(\mathbf{u};\bm{\theta})=u_{n+1}-F_{n}(u_{n};{\theta}),$

where $\bm{\theta}=(\theta_{1},\theta_{2},\dots,\theta_{q})^{T}$ is an
uncertain parameter. Then the minimization problem becomes

$\min_{\\{\mathbf{u}\in\mathbb{R}^{nN},\
\bm{\theta}\in\mathbb{R}^{q}\\}}\\{\|G(\mathbf{u};\bm{\theta})\|^{2}+\alpha\|H\mathbf{u}-\mathbf{y}\|^{2}\\}.$

Starting with an initial guess for the state $\mathbf{u}^{(0)}$ and parameters
$\bm{\theta}^{(0)}$, the iteration proceeds as follows:

$\displaystyle\mathbf{u}^{(k+1)}$ $\displaystyle=$
$\displaystyle\mathbf{u}^{(k)}-\mathcal{L}(\mathbf{u}^{(k)};\bm{\theta}^{(k)}),$
(16) $\displaystyle\bm{\theta}^{(k+1)}$ $\displaystyle=$
$\displaystyle\bm{\theta}^{(k)}-\mathcal{S}(\mathbf{u}^{(k+1)};\bm{\theta}^{(k)}).$
(17)

Here $\mathcal{L}(\mathbf{u}^{(k)};\cdot)$ and
$\mathcal{S}(\cdot;\bm{\theta}^{(k)})$ are defined by

$\displaystyle\mathcal{L}(\mathbf{u}^{(k)};\cdot)$ $\displaystyle=$
$\displaystyle\left(G_{u}^{{}^{\prime}T}(\mathbf{u}^{(k)};\cdot)G_{u}^{{}^{\prime}}(\mathbf{u}^{(k)};\cdot)+\alpha
H^{T}H\right)^{-1}\bigg{(}G_{u}^{{}^{\prime}T}(\mathbf{u}^{(k)};\cdot)G(\mathbf{u}^{(k)};\cdot)+\alpha
H^{T}(H\mathbf{u}^{(k)}-\mathbf{y})\bigg{)},$
$\displaystyle\mathcal{S}(\cdot;\bm{\theta}^{(k)})$ $\displaystyle=$
$\displaystyle\left(G_{\theta}^{{}^{\prime}T}(\cdot;\bm{\theta}^{(k)})G_{\theta}^{{}^{\prime}}(\cdot;\bm{\theta}^{(k)})\right)^{-1}G_{\theta}^{{}^{\prime}T}(\cdot;\bm{\theta}^{(k)})G(\cdot;\bm{\theta}^{(k)}),$

respectively, $G^{\prime}_{u}$ and $G^{\prime}_{\theta}$ are derivatives of
$G$ with respect to $u$ and $\theta$, respectively.

######  Assumption 4.1.

We consider $G(\mathbf{u};\bm{\theta})=\mathcal{G}(\mathbf{u})+A\bm{\theta}$.
$\mathcal{G}$ is continuously differentiable in the open convex set
$D\subset\mathbb{R}^{nN}$ and Lipshitz in $D$, with the Lipschitz constant
denoted by $L_{0}>0$, namely

$\|\mathcal{G}(\mathbf{u})-\mathcal{G}(\mathbf{v})\|\leq
L_{0}\|\mathbf{u}-\mathbf{v}\|,\quad\forall\mathbf{u},\mathbf{v}\in D.$ (18)

The Jacobian $\mathcal{G}^{\prime}(\mathbf{u})$ is locally Lipschitz
continuous in $D$, with the Lipschitz constant denoted by $L_{3}>0$, namely

$\|\mathcal{G}^{\prime}(\mathbf{u})-\mathcal{G}^{\prime}(\mathbf{v})\|\leq
L_{3}\|\mathbf{u}-\mathbf{v}\|,\quad\forall\mathbf{u},\mathbf{v}\in D,$ (19)

and there exists $\alpha>0$ such that

$\|(\mathcal{G}^{\prime}(\mathbf{x})^{T}\mathcal{G}^{\prime}(\mathbf{x})+\alpha
H^{T}H)^{-1}\mathcal{G}^{\prime
T}(\mathbf{x})\|\leq(2L_{3}c)^{-1},\quad\forall\mathbf{x}\in D,$ (20)

with $L_{3}$ being the Lipschitz constant from Equation (19), and $c>0$ being
an upper bound on the norm of initial error $\|\mathbf{e}_{0}\|$. Furthermore,
$b/c<1$ where $b=\|A(A^{T}A)^{-1}A^{T}\|L_{0}/L_{3}$.

######  Theorem 3.

We consider $G(\mathbf{u};\bm{\theta})=\mathcal{G}(\mathbf{u})+A\bm{\theta}$
and noise-free observations of $\mathbf{u}$. We assume that Assumption 4.1
holds for $G(\mathbf{u};\bm{\theta})$. Then

$\limsup_{k\to\infty}\|\mathbf{e}_{k}\|<\frac{b/2}{1-b/c},$

and

$\limsup_{k\to\infty}\|\bm{\theta}^{(k)}-\bm{\theta}^{\dagger}\|<L_{0}\|(A^{T}A)^{-1}A^{T}\|\frac{b/2}{1-b/c}.$

Theorem (3) is proven in (B).

##  5 NUMERICAL EXPERIMENTS

In this section, we present numerical experiments to illustrate the
theoretical results of Section 3 and Section 4. First, we consider noise-free
observations. We illustrate the convergence Theorem (1) under different sizes
of observations. Second, we consider noisy observations and show the
theoretical bound (Equation (14)), computed numerically. Third, we consider a
dynamical model with uncertain parameters and estimate them using the Gauss-
Newton DA method (Equations (16) and (17)). For each set-up, we perform 100
numerical experiments with different realizations of truth
$\mathbf{u}^{\dagger}$, observations $\mathbf{y}$, and background solution
$\mathbf{u}^{b}$. For each solution $\mathbf{u}$, we compute cost function,
error with respect to the truth, error with respect to the truth of observed
variables, and error with respect to the truth of non-observed variables :

$\displaystyle C$
$\displaystyle=\|G(\mathbf{u})\|+\alpha\|\mathbf{y}-H\mathbf{u}\|,$ (21)
$\displaystyle\mathcal{E}$
$\displaystyle=\|\mathbf{u}^{\dagger}-\mathbf{u}\|,$ (22)
$\displaystyle\mathcal{E}^{O}$
$\displaystyle=\|H\mathbf{u}^{\dagger}-H\mathbf{u}\|,$ (23)
$\displaystyle\mathcal{E}^{N}$
$\displaystyle=\|(I-H^{T}H)(\mathbf{u}^{\dagger}-\mathbf{u})\|,$ (24)

respectively.

We compare the Gauss-Newton DA method to WC4DVar, which minimizes the
following cost function:

$J(u_{0};\\{y_{n}\\})=\frac{1}{2}\sum_{n=1}^{N}(y_{n}-Hu_{n})^{T}R^{-1}(y_{n}-Hu_{n})+\frac{1}{2}(u_{n}-F_{n}(u_{n-1})^{T}Q^{-1}(u_{n}-F_{n}(u_{n-1}),$
(25)

where $R$ is covariance matrix of observational error and $Q$ is covariance
matrix of model error (see, e.g., [2, 3, 4, 5, 18]). The main distinction of
the Gauss-Newton DA method from WC4DVar is in the tunable parameter $\alpha$,
which has a significant influence on error convergence or error boundedness.
The minimization of the WC4DVar cost function is done by a Matlab built-in
Levenberg-Marquardt algorithm.

We perform numerical experiments with the Lorenz 63 (L63) and Lorenz 96 (L96)
models. L63 is a chaotic model which is widely used as a toy model in data-
assimilation numerical experiments. It simulates atmospheric convection in a
simple way [22]. The model is described by the following ODEs

$\dot{x}_{1}=\sigma(x_{2}-x_{1}),\quad\dot{x}_{2}=x_{1}(\rho-
x_{3})-x_{2},\quad\dot{x}_{3}=x_{1}x_{2}-bx_{3}.$ (26)

We implement the L63 model with the standard parameters, $\sigma=10$,
$\rho=28$, and $b=8/3$. The differential equations are discretized with a
forward Euler scheme with time step $\Delta t=0.005$. The initial conditions
are random numbers, which are independently and identically distributed. We
generate observations by computing a solution of L63 on $t\in[0,100]$. The
observations are drawn at every tenth time step, which corresponds to
assimilation time window of 6 hours. For the Euler-discretized L63 model, the
Lipschitz condition is

$\|G^{\prime}(X)-G^{\prime}(Y)\|\leq\sqrt{2}\Delta t\|X-Y\|.$

(For derivation of the Lipschitz constant see A.1).

The L96 model [23] is a one-dimensional atmosphere model which is described by
the following ODEs

$\dot{x}_{l}=-x_{l-2}x_{l-1}+x_{l-1}x_{l+1}-x_{l}+\mathcal{F},\quad
l=1,\dots,d,$ (27)

where the dimension $d$ and forcing $\mathcal{F}$ are parameters. Cyclic
boundary conditions are imposed. We implement the L96 model with the standard
parameter choices $d=40$ and $\mathcal{F}=8$. We discretize the differential
equations with a forward Euler scheme with time step $\Delta t=0.0025$. The
initial conditions are random numbers which are independently and identically
distributed. We generate observations by computing a solution of L96 on
$t\in[0,100]$. The observations are drawn at every tenth time step, which
corresponds to assimilation time window of 3 hours. The Lipschitz condition of
the Euler-discretized L96 model is:

$\|G^{\prime}(X)-G^{\prime}(Y)\|\leq\sqrt{6}\Delta t\|X-Y\|.$

(For derivation of the Lipschitz constant see A.2).

###  5.1 State estimation given noise-free observations

We perform numerical experiments with noise-free observations, thus
$\eta_{j}=0,\ j=0,\cdots,N$ in Equation (3). In order to satisfy the
conditions of Theorem (1), we first need to compute the Lipschitz constant of
the Jacobian of the dynamical model. Next, we need to find a suitable
$\alpha$. We use Algorithm 1, where an upper bound on the initial error, $c$,
is chosen arbitrarily.

Algorithm 1 Finding parameter $\alpha$ in case of noise-free observations

Given the initial guess of the Gauss-Newton DA method, $\mathbf{u}^{(0)}$, an
upper bound on the initial error, $c$, and the Lipschitz constant $L_{1}$ that
satisfies Equation (9) in Assumption (3.1), we choose an arbitrary positive
$\alpha_{0}$, for example $\alpha_{0}=0.001$.

1:while ${\|(G^{\prime
T}(\mathbf{u}^{(0)})G^{\prime}(\mathbf{u}^{(0)})+\alpha_{0}H^{T}H)^{-1}G^{\prime
T}(\mathbf{u}^{(0)})\|>\left(L_{1}c\right)^{-1}}$ and $(\sim
isnan(\|(G^{\prime}(\mathbf{u}^{(0)})^{T}G^{\prime}(\mathbf{u}^{(0)})+\alpha_{0}H^{T}H)^{-1}G^{\prime
T}(\mathbf{u}^{(0)})\|))$ do

2: $\alpha_{0}\leftarrow 2*\alpha_{0}$

3:end while

4:if
$isnan(\|(G^{\prime}(\mathbf{u}^{(0)})^{T}G^{\prime}(\mathbf{u}^{(0)})+\alpha_{0}H^{T}H)^{-1}G^{\prime
T}(\mathbf{u}^{(0)})\|)$ then

5: Error: There is no $\alpha$.

6:else $\alpha\leftarrow\alpha_{0}$

7:end if

We use the same value of $\alpha$ throughout the iteration and check whether
Equation (10) is satisfied. If there exists an iteration $k$ such that
Equation (10) does not hold we terminate the iteration, otherwise the
iteration proceeds until a tolerance value $(10^{-14})$ is reached for the
change in the solution $\|\mathbf{u}^{(k+1)}-\mathbf{u}^{(k)}\|$.

We consider the L63 model and assume only the first component is observed,
i.e., $H=[1,0,0]$. In all of the experiments Equation (10) is satisfied
throughout the iteration. In Figure 1(a), we display error with respect to the
truth (Equation (22)) as a function of iteration. We plot median and plus and
minus one standard deviation over 100 simulations. We see that error is a
decreasing function, as we expect from Theorem (1). We also considered other
observation operators which depend only on second or third variable and
obtained equivalent results (not shown). Therefore, for the L63 model the
convergence Theorem (1) does not depend which variable is observed.

(a)

(b)

Figure 1: Application of the Gauss-Newton DA method to L63 with noise-free
observations, where only the first variable is observed. On the left: Error
(Equation (22)) as a function of iteration: median (dashed line), +/- one
standard deviation (shadowed area) over 100 simulations. On the right: Cost
function (Equation (21)) as a function of iteration for different values of
$\alpha$.

Furthermore, We perform numerical experiments with different values of
$\alpha$ that satisfy Equation (10). In Figure 1(b), we display cost function
(Equation (21)) as a function of iteration for different values of $\alpha$. A
faster convergence rate is achieved at $\alpha=0.004$, as to be expected since
$\alpha=0$ corresponds to the second order Newton method when observations are
complete.

Next, we perform numerical experiments with the L96 model and we assume every
second variable is observed. In Figure 2(a), we display error with respect to
the truth as a function of iteration. We plot median and plus and minus one
standard deviation over 100 simulations. We see that error is a decreasing
function, as we expect from Theorem (1). Comparing to Figure 1(a) we observe
that more iterations are needed to reach the desired tolerance in the L96
model due to the higher dimension.

(a)

(b)

Figure 2: Application of the Gauss-Newton DA method to L96 with noise-free
observations, where every second variable is observed. On the left: Error
(Equation (22)) as a function of iteration: median (dashed line), +/- one
standard deviation (shadowed area) over 100 simulations. On the right: Cost
function (Equation (21)) as a function of iteration for different values of
$\alpha$.

In Figure 2(b), we plot the cost function of the L96 model for different
values of $\alpha$ when observing every second variable. A faster convergence
rate is achieved at $\alpha=0.004$. Furthermore, we investigate the
convergence property under different sizes of observations. In Figure 3, we
plot the median value of error (Equation (22)) over 20 simulations as a
function of iteration in a semi-log scale. We see that when the size of
observations decreases more iterations are needed to obtain the desired
tolerance.

Figure 3: Application to L96 with noise-free observations. Error of the Gauss-
Newton DA method as a function of iterations for different sizes of
observations.

Furthermore, Equation (10) of Assumption (3.1) gives implicit requirement on
the size of observations. In numerical experiments with the L96 model, we see
that if the number of observations is less than eight then the matrix
$(G^{\prime}(\mathbf{x})^{T}G^{\prime}(\mathbf{x})+\alpha H^{T}H)^{-1}$ is
ill-posed and we are not able to find $\alpha$ that satisfies Equation (10).
Therefore, the conditions in Theorem (1) do not hold and error convergence is
not guaranteed.

###  5.2 State estimation given noisy observations

We perform numerical experiments with noisy observations, where
$\eta\sim\mathcal{N}(0,\gamma^{2}I)$ in Equation (3). We assume that
$\|H^{T}\eta\|$ is bounded and known, which is a strong assumption. In order
to satisfy the conditions of Theorem (2), we first need to compute the
Lipschitz constant of the Jacobian of the dynamical model, and then to find a
suitable $\alpha$. We use Algorithm 2 to find $\alpha$, where an upper bound
on the initial error, $c$, is chosen arbitrarily.

Algorithm 2 Finding parameter $\alpha$ in case of noisy observations

Given the initial guess of the Gauss-Newton DA method, $\mathbf{u}^{(0)}$, an
upper bound on the initial error, $c$, and the Lipschitz constant $L_{2}$ that
satisfies Equation (11) in Assumption (3.2), we choose an arbitrary positive
$\alpha_{0}$, for example $\alpha_{0}=0.001$.

1:while ${\|(G^{\prime
T}(\mathbf{u}^{(0)})G^{\prime}(\mathbf{u}^{(0)})+\alpha_{0}H^{T}H)^{-1}G^{\prime
T}(\mathbf{u}^{(0)})\|>(L_{2}c)^{-1}}$ and ${\|H^{T}\eta\|\|(G^{\prime
T}(\mathbf{u}^{(0)})G^{\prime}(\mathbf{u}^{(0)})+\alpha_{0}H^{T}H)^{-1}\|>c/2}$
and $(\sim
isnan(\|(G^{\prime}(\mathbf{u}^{(0)})^{T}G^{\prime}(\mathbf{u}^{(0)})+\alpha_{0}H^{T}H)^{-1}G^{\prime
T}(\mathbf{u}^{(0)})\|))$ do

2: $\alpha_{0}\leftarrow 2*\alpha_{0}$

3:end while

4:if
$isnan(\|(G^{\prime}(\mathbf{u}^{(0)})^{T}G^{\prime}(\mathbf{u}^{(0)})+\alpha_{0}H^{T}H)^{-1}G^{\prime
T}(\mathbf{u}^{(0)})\|)$ or $\alpha_{0}>1$ then

5: Error: There is no $\alpha$.

6:else $\alpha\leftarrow\alpha_{0}$

7:end if

We use the same value of $\alpha$ throughout the iteration and check whether
Equations (12) and (13) are satisfied. If there exists an iteration $k$ such
that either Equation (12) or Equation (13) does not hold, we terminate the
iteration, otherwise the iteration proceeds until it reaches the maximum value
of iterations, which is 20 for L63 and 70 for L96.

We consider the L63 model with noisy observations. We assume only the first
component is observed, i.e., $H=[1,0,0]$ and $\gamma=0.01$ in Equation (3). In
all of the experiments Equations (12) and (13) are satisfied throughout the
iteration. In Figure 4, we display error with respect to the truth (Equation
(22)) as a function of iteration. We plot median (dashed line) and plus and
minus one standard deviation (shadowed area) over 100 simulations. In Figure
4, we also display theoretical bound (Equation (14)), we plot median (green
line) over 100 simulations which depends on different $\alpha$s. We see that
error is below the theoretical bound, as we expect from Theorem (2). On the
right of Figure 4, we display error of observed components (Equation (23)) as
a function of iteration. We plot median and plus and minus one standard
deviation over 100 simulations. We see that error of observed components
(dashed line) is less than observation error
$\|\mathbf{y}-H\mathbf{u}^{\dagger}\|$ (red line). Furthermore, comparing the
right panel to the left panel of Figure 4 we see that the error of observed
components is less than the total error.

Next, we perform numerical experiments with the L96 model with noisy
observations. We assume every second variable is observed and $\gamma=0.01$ in
Equation (3). On the left of Figure 5, we display error with respect to the
truth as a function of iteration. We plot median and plus and minus one
standard deviation over 100 simulations. We see that error is a decreasing
function for the first five iterations and after that it remains bounded.
Comparing to Figure 4 we observe that due to the higher dimensions, more
iterations are needed to reach the desired tolerance in the L96 model.
Moreover, as we expect from Theorem (2), the median of error is below the
median of theoretical bound (green line) over 100 simulations. On the right of
Figure 5, we display error of observed components (Equation (23)) as a
function of iteration. We see that error of observed components is less than
the total error (left panel of Figure 5), as it is expected. Furthermore, we
see from both Figures 4 and 5 that the theoretical bound is not tight as to be
expected.

Figure 4: Application of the Gauss-Newton DA method to L63 with noisy
observations, where only the first variable is observed and $\gamma=0.01$ in
Equation (3). On the left: Error (Equation (22)) as a function of iteration:
median (dashed line), +/- one standard deviation (shadowed area) over 100
simulations and a theoretical bound (Equation (14)) in green. On the right:
Error of observed components (Equation (23)) as a function of iteration and
observation error in red. Figure 5: Application of the Gauss-Newton DA method
to L96 with noisy observations, where every second variable is observed and
$\gamma=0.01$ in Equation (3). On the left: Error (Equation (22)) as a
function of iteration: median (dashed line), +/- one standard deviation
(shadowed area) over 100 simulations and a theoretical bound Equation (14) in
green. On the right: Error of observed components (Equation (23)) and
observation error in red.

Finally, we investigate error (Equation (22)) as a function of observation
noise level $\gamma$. From Theorem (2) follows that the convergence result is
achieved as $\|H^{T}\eta\|$ goes to zero. In Figure 6, we display error
(Equation (22)) for different values of $\gamma$ for the L63 model on the left
and for the L96 model on the right. As we expected from Remark (3.4), we see
that error decreases as observation noise level decreases and the convergence
result achieves when $\gamma$ is close to zero.

Figure 6: Error with respect to the truth for different observation noise
level. On the left: Application to L63 with noisy observations, we observe
only the first variable. On the right: Application to L96 with noisy
observations, we observe every second variable.

###  5.3 Comparison to weak-constraint 4DVar

Weak-constraint four-dimensional variational method is one of the well-known
data-assimilation methods for estimating initial condition in weather
forecasting applications. It minimizes the cost function (Equation (25)) under
assumption of imperfect model dynamics which is also the goal of the Gauss-
Newton DA method. We compare Gauss-Newton DA method to WC4DVar method. We
perform numerical experiments using the L63 and L96 models with the same
parameters as in the previous section. In these experiments, we use identical
data, models, and windows for both methods.

In Figure 7 and Figure 8 we plot errors for the L63 and L96 models,
respectively. On the left of the figures, we plot errors with respect to the
truth of observed variables (Equation (23)). On the right of the figures, we
plot errors with respect to the truth of non-observed variables (Equation
(24)). We see that error of Gauss-Newton DA method is significantly less than
the error of WC4DVar method for both observed variables and non-observed
variables. We see that the error of both methods is below the observation
error.

Figure 7: Application to L63. Error as a function of time: median (dashed
line) +/- one standard deviation over 20 simulations with length of
assimilation window, 2.5. On the left: error with respect to the truth of
observed variables. On the right: error with respect to the truth of non-
observed variables. The Gauss-Newton DA method is in grey, WC4DVar method is
in blue, and the observational error is in red. Figure 8: Application to L96.
Error as a function of time: median (dashed line) +/- one standard deviation
over 20 simulations with length of assimilation window, 1.25. On the left:
error with respect to the truth of observed variables. On the right: error
with respect to the truth of non-observed variables. The Gauss-Newton DA
method is in grey, WC4DVar method is in blue, and the observational error is
in red.

###  5.4 Parameter estimation

As described in Section 4, the Gauss-Newton DA method can be applied to the
problem of joint state-parameter estimation. In these experiments, we use the
Gauss-Newton DA method (Equations (16) and (17)) to estimate a parameter
$\sigma$ of the L63 model which we assume to be uncertain. Different values of
initial $\sigma$ were chosen, 5, 15, and 20, with the true $\sigma$ being 10.
We alternate between Equations (16) and (17) until a termination condition is
satisfied. The termination condition is satisfied when either the number of
iterations reaches the maximum (500) or the distance between two successive
approximations of the uncertain parameter is less than desired tolerance
$(10^{-3})$. We consider both noisy and noise-free observations and assume the
first and second components of the state are observed. For each numerical
experiment and fixed initial $\sigma$, we perform 100 realizations. In Table I
and II, we display the error of state estimation, which is the median of
Equation (22) over 100 simulations, and estimated $\sigma$ for different
initial parameters. Table I is for noise-free observations and Table II is for
noisy observations with $\gamma=\sqrt{3}$.

In Table I, we consider noise-free observations. We assume the first and
second components of the state are observed. The results of $\sigma$
estimation are shown in Table I.

In Table II, we consider noisy observations. The observations of state are
obtained from the true state by adding i.i.d zero mean Gaussian noise as in
Equation (3) with covariance matrix $\Gamma=\gamma^{2}I$, where $\gamma^{2}$
is the variance of the noise procedure. In this numerical experiment, we
consider noisy observations with $\gamma=\sqrt{3}$ in Equation (3). It should
be noted that similar results can be obtained for $\rho$ or $\beta$ estimation
of the L63 model (not shown).

Table I: Estimation of $\sigma$ by the Gauss-Newton DA method for noise-free observations, we observe the first and second components of the state. Property |  Initial guess for $\sigma$
---|---
Initial $\sigma$ | 5 | 15 | 20
Error between estimated and true solution | $0.4025$ | $0.4256$ | $0.4267$
Estimated $\sigma$ | 9.7465 | 10.2746 | 10.2785
Table II: Estimation of $\sigma$ by Gauss-Newton DA method for noisy observations. Property |  Initial guess for $\sigma$
---|---
Initial $\sigma$ | 5 | 15 | 20
Observation error | 3
Error between estimated and true solution | $2.7287$ | $6.7806$ | $3.8887$
Estimated $\sigma$ | 9.7380 | 10.4505 | 10.8239

##  6 CONCLUSION

We considered a cost function of a data-assimilation problem consisting of
observation mismatch, $\|H\mathbf{u}-\mathbf{y}\|$ weighted by a parameter
$\alpha$ and model mismatch $\|G(\mathbf{u})\|$. By solving this minimization
problem approximately using an iterative algorithm, Gauss-Newton DA method, we
found an estimate for a considered system. One of the advantage of this cost
function is using the information of available observations in each iteration
which can help to find a more accurate approximated solution for a physical
system. We then established error bounds with respect to the truth under some
conditions in two cases, noise-free and noisy observations. In case of noise-
free observations, we proved the method is convergent and we obtained an upper
bound in case of noisy observations. In numerical experiments, we applied
Gauss-Newton DA method for two models, Lorenz 63 and Lorenz 96. We observed
the numerical evidence confirmed the theoretical results and we demonstrated
that the results compare favorably to those of WC4DVar. Moreover, we showed
that how the size of observations could affect on the convergence results. We
observed that if the size of observations is less than eight in the Lorenz 96
model, we do not have a convergence result. Therefore future directions
include finding a method for deriving an approximation based on low-order of
observations.

##  Appendix A Lipschitz constants for the Lorenz 63 and the Lorenz 96 models

Here we compute the Lipschitz constant for the Lorenz 63 and the Lorenz 96
models discretized with forward Euler.

###  A.1 Lipschitz constant for the L63 model discretized with forward Euler

Consider the discretized form of Equation (26)

$\displaystyle X_{1}(t_{k+1})$ $\displaystyle=$ $\displaystyle
X_{1}(t_{k})+\Delta t\sigma(X_{2}(t_{k})-X_{1}(t_{k})),$ (28) $\displaystyle
X_{2}(t_{k+1})$ $\displaystyle=$ $\displaystyle X_{2}(t_{k})+\Delta
t\left(X_{1}(t_{k})(\rho-X_{3}(t_{k}))-X_{2}(t_{k})\right),$ $\displaystyle
X_{3}(t_{k+1})$ $\displaystyle=$ $\displaystyle X_{3}(t_{k})+\Delta
t\left(X_{1}(t_{k})X_{2}(t_{k})-bX_{3}(t_{k})\right),$

with the following matrix notation

$\mathbf{X}_{k+1}=F(\mathbf{X}_{k}),\quad
G(\mathbf{X}):=\mathbf{X}_{k+1}-F(\mathbf{X}_{k})$

where $\mathbf{X}_{k}=\left(X_{1}(t_{k}),X_{2}(t_{k}),X_{3}(t_{k})\right)^{T}$
and $F(\mathbf{X}_{k})$ is the right-hand side of (28). From the definition of
Jacobian of $G$ we have

$G^{\prime}(\mathbf{X})-G^{\prime}(\mathbf{Y})=\begin{bmatrix}-F^{\prime}(\mathbf{X}_{0})+F^{\prime}(\mathbf{Y}_{0})&&\\\
&-F^{\prime}(\mathbf{X}_{1})+F^{\prime}(\mathbf{Y}_{1})&&\\\ &&\ddots&&\\\
&&&-F^{\prime}(\mathbf{X}_{N-1})+F^{\prime}(\mathbf{Y}_{N-1})&\\\
\end{bmatrix}.$

with

$F^{\prime}(\mathbf{X}_{k})=\begin{pmatrix}1-\sigma\Delta t&\sigma\Delta
t&0\\\ \Delta t(\rho-X_{3}(t_{k}))&1-\Delta t&-\Delta tX_{1}(t_{k})\\\ \Delta
tX_{2}(t_{k})&\Delta tX_{1}(t_{k})&1-b\Delta t\\\ \end{pmatrix}.$ (29)

From (29) we get

$\displaystyle\|F^{\prime}(\mathbf{X}_{k})-F^{\prime}(\mathbf{Y}_{k})\|^{2}_{F}=\sum_{j=1}^{n}\sum_{i=1}^{m}|F_{ij}|^{2}\
$ $\displaystyle=\Delta
t^{2}\bigg{(}|Y_{3}(t_{k})-X_{3}(t_{k})|^{2}+|Y_{1}(t_{k})-X_{1}(t_{k})|^{2}$
$\displaystyle+|Y_{2}(t_{k})-X_{2}(t_{k})|^{2}+|Y_{1}(t_{k})-X_{1}(t_{k})|^{2}\bigg{)}$
$\displaystyle\leq 2\Delta
t^{2}\|\mathbf{X}_{k}-\mathbf{Y}_{k}\|^{2}_{L^{2}}.$ (30)

Therefore, we have

$\|G^{\prime}(\mathbf{X})-G^{\prime}(\mathbf{Y})\|_{F}\leq\sqrt{2}\Delta
t\|\mathbf{X}-\mathbf{Y}\|_{L^{2}}.$

Since $\|A\|_{L^{2}}\leq\|A\|_{F}$, by substituting it into the expression
above we obtain the following

$\|G^{\prime}(\mathbf{X})-G^{\prime}(\mathbf{Y})\|_{L^{2}}\leq\sqrt{2}\Delta
t\|\mathbf{X}-\mathbf{Y}\|_{L^{2}}.$

###  A.2 Lipschitz constant for the L96 model discretized with forward Euler

Consider the discretized form of Equation (27)

$X_{l}(t_{k+1})=X_{l}(t_{k})+\left(-X_{l-2}(t_{k})X_{l-1}(t_{k})+X_{l-1}(t_{k})X_{l+1}(t_{k})-X_{l}(t_{k})+\mathcal{F}\right)\Delta
t,\ l=1,\dots,40.$ (31)

with the following matrix notation

$\mathbf{X}_{k+1}=F(\mathbf{X}_{k}),\qquad
G(\mathbf{X}):=\mathbf{X}_{k+1}-F(\mathbf{X}_{k})$

where $\mathbf{X}_{k}=\left(X_{1}(t_{k}),\dots,X_{40}(t_{k})\right)^{T}$ and
$F(\mathbf{X}_{k})$ is the right-hand side of (31). From the definition of
Jacobian of $G$ we have

$G^{\prime}(\mathbf{X})-G^{\prime}(\mathbf{Y})=\begin{bmatrix}-F^{\prime}(\mathbf{X}_{0})+F^{\prime}(\mathbf{Y}_{0})&&\\\
&-F^{\prime}(\mathbf{X}_{1})+F^{\prime}(\mathbf{Y}_{1})&&\\\ &&\ddots&&\\\
&&&-F^{\prime}(\mathbf{X}_{N-1})+F^{\prime}(\mathbf{Y}_{N-1})&\\\
\end{bmatrix}.$

with

$F^{\prime}(\mathbf{X}_{k})=\begin{pmatrix}1-\Delta t&\Delta
tX_{40}&0&\dots&0&-\Delta tX_{40}&\Delta t(X_{2}-X_{39})\\\ \Delta
t(X_{3}-X_{40})&1-\Delta t&\Delta tX_{1}&0&\dots&0&-\Delta tX_{1}\\\
\vdots&\dots&&\dots&&&\vdots\\\ 0&\dots&0&-\Delta X_{38}&\Delta
t(X_{40}-X_{37})&1-\Delta t&\Delta tX_{38}\\\ \Delta tX_{39}&0&\dots&0&-\Delta
tX_{39}&\Delta t(X_{1}-X_{38})&1-\Delta t\\\ \end{pmatrix}.$

From the expression above we get

$\displaystyle\|F^{\prime}(\mathbf{X}_{k})-F^{\prime}(\mathbf{Y}_{k})\|_{F}^{2}$
$\displaystyle=\Delta t^{2}|X_{40}-Y_{40}|^{2}+\Delta
t^{2}|X_{40}-Y_{40}|^{2}$ $\displaystyle+\Delta
t^{2}|(X_{2}-Y_{2})-(X_{39}-Y_{39})|^{2}$ $\displaystyle+\Delta
t^{2}|(X_{3}-Y_{3})-(X_{40}-Y_{40})|^{2}$ $\displaystyle+\Delta
t^{2}|X_{1}-Y_{1}|^{2}+\Delta t^{2}|X_{1}-Y_{1}|^{2}$ $\displaystyle+...$
$\displaystyle+\Delta t^{2}|X_{38}-Y_{38}|^{2}+\Delta
t^{2}|(X_{40}-Y_{40})-(X_{37}-Y_{37})|^{2}$ $\displaystyle+\Delta
t^{2}|X_{38}-Y_{38}|^{2}+\Delta t^{2}|X_{39}-Y_{39}|^{2}$
$\displaystyle+\Delta t^{2}|X_{39}-Y_{39}|^{2}+\Delta
t^{2}|(X_{1}-Y_{1})-(X_{38}-Y_{38})|^{2}$ $\displaystyle\leq 2\Delta
t^{2}\|\mathbf{X}_{k}-\mathbf{Y}_{k}\|^{2}_{L^{2}}$ $\displaystyle+\\{2\Delta
t^{2}|X_{2}-Y_{2}|^{2}+2\Delta t^{2}|X_{39}-Y_{39}|^{2}+2\Delta
t^{2}|X_{3}-Y_{3}|^{2}$ $\displaystyle+2\Delta
t^{2}|X_{40}-Y_{40}|^{2}+...+2\Delta t^{2}|X_{40}-Y_{40}|^{2}$
$\displaystyle+2\Delta t^{2}|X_{37}-Y_{37}|^{2}+2\Delta
t^{2}|X_{1}-Y_{1}|^{2}+2\Delta t^{2}|X_{38}-Y_{38}|^{2}\\}$ $\displaystyle\leq
2\Delta t^{2}\|\mathbf{X}_{k}-\mathbf{Y}_{k}\|^{2}_{L^{2}}+2\Delta
t^{2}(2\|\mathbf{X}_{k}-\mathbf{Y}_{k}\|^{2}_{L^{2}})$ $\displaystyle=6\Delta
t^{2}\|\mathbf{X}_{k}-\mathbf{Y}_{k}\|^{2}_{L^{2}}.$

On the other hand,

$\|G^{\prime}(\mathbf{X}_{k})-G^{\prime}(\mathbf{Y}_{k})\|^{2}_{L^{2}}\leq\|G^{\prime}(\mathbf{X}_{k})-G^{\prime}(\mathbf{Y}_{k})\|^{2}_{F}=\sum_{j=1}^{n}\sum_{i=1}^{m}|F_{ij}|^{2}.$

Combining the two above expressions leads to

$\|G^{\prime}(\mathbf{X})-G^{\prime}(\mathbf{Y})\|_{L^{2}}\leq\sqrt{6}\Delta
t\|\mathbf{X}-\mathbf{Y}\|_{L^{2}}.$

##  Appendix B Bound on the parameter estimation error

Here we prove Theorem (3).

###### Proof.

By setting
$G(\mathbf{u}^{(k)};\bm{\theta}^{(k-1)})=\mathcal{G}(\mathbf{u}^{(k)})+A\bm{\theta}^{(k-1)}$
in Equation (17), we get

$\bm{\theta}^{(k)}=\bm{\theta}^{(k-1)}-(A^{T}A)^{-1}A^{T}\left(\mathcal{G}(\mathbf{u}^{(k)})+A\bm{\theta}^{(k-1)}\right)$

Let us define $\mathbf{E}_{k}=\bm{\theta}^{(k)}-\bm{\theta}^{\dagger}$, then
we have

$\displaystyle\mathbf{E}_{k}$
$\displaystyle=\mathbf{E}_{k-1}-(A^{T}A)^{-1}A^{T}\left(\mathcal{G}(\mathbf{u}^{(k)})+A\bm{\theta}^{(k-1)}\right)$
$\displaystyle=\mathbf{E}_{k-1}-(A^{T}A)^{-1}A^{T}\left(\mathcal{G}(\mathbf{u}^{(k)})+A\bm{\theta}^{(k-1)}-\mathcal{G}(\mathbf{u}^{\dagger})-A\bm{\theta}^{\dagger}\right),$

where we used
$G(\mathbf{u}^{\dagger};\bm{\theta}^{\dagger})=\mathcal{G}(\mathbf{u}^{\dagger})+A\bm{\theta}^{\dagger}=0$.
Therefore, we conclude that

$\mathbf{E}_{k}=-(A^{T}A)^{-1}A^{T}\left(\mathcal{G}(\mathbf{u}^{(k)})-\mathcal{G}(\mathbf{u}^{\dagger})\right).$
(32)

On the other hand, using Equation (16) in the definition of the state error
$\mathbf{e}_{k}=\mathbf{u}^{k}-\mathbf{u}^{\dagger}$ leads to

$\displaystyle\mathbf{e}_{k+1}=\mathbf{e}_{k}-\left(\mathcal{G}^{\prime
T}\mathcal{G}^{\prime}+\alpha H^{T}H\right)^{-1}\left(\mathcal{G}^{\prime
T}\left(\mathcal{G}(\mathbf{u}^{(k)})+A\bm{\theta}^{(k)}\right)+\alpha
H^{T}(H\mathbf{u}^{(k)}-\mathbf{y})\right),$

where $\mathcal{G}^{\prime}$ and $\mathcal{G}^{\prime T}$ are computed at
$\mathbf{u}^{(k)}$. By substituting Equation (3) with $\eta_{j}=0,\ \forall
j=0,\cdots,N$ into the expression above, we obtain the following

$\displaystyle\mathbf{e}_{k+1}$ $\displaystyle=$
$\displaystyle\mathbf{e}_{k}-\left(\mathcal{G}^{\prime
T}\mathcal{G}^{\prime}+\alpha H^{T}H\right)^{-1}\left(\mathcal{G}^{\prime
T}\left(\mathcal{G}(\mathbf{u}^{(k)})+A\bm{\theta}^{(k)}\right)+\alpha
H^{T}H\mathbf{e}_{k}\right)$ $\displaystyle=$
$\displaystyle\mathbf{e}_{k}-\left(\mathcal{G}^{\prime
T}\mathcal{G}^{\prime}+\alpha H^{T}H\right)^{-1}\left(\mathcal{G}^{\prime
T}\left(\mathcal{G}(\mathbf{u}^{(k)})+A\bm{\theta}^{(k)}-\mathcal{G}(\mathbf{u}^{\dagger})-A\bm{\theta}^{\dagger}\right)+\alpha
H^{T}H\mathbf{e}_{k}\right),$

where we used $\mathcal{G}(\mathbf{u}^{\dagger})+A\bm{\theta}^{\dagger}=0$.
Opening the brackets and using the following property

$\displaystyle I-(\mathcal{G}^{{}^{\prime}T}\mathcal{G}^{\prime}+\alpha
H^{T}H)^{-1}\alpha
H^{T}H=(\mathcal{G}^{{}^{\prime}T}\mathcal{G}^{\prime}+\alpha
H^{T}H)^{-1}\mathcal{G}^{{}^{\prime}T}\mathcal{G}^{\prime},$

we get

$\mathbf{e}_{k+1}=\left(\mathcal{G}^{\prime T}\mathcal{G}^{\prime}+\alpha
H^{T}H\right)^{-1}\mathcal{G}^{\prime
T}\left(\mathcal{G}^{\prime}\mathbf{e}_{k}-\mathcal{G}(\mathbf{u}^{(k)})+\mathcal{G}(\mathbf{u}^{\dagger})-AE_{k}\right).$

Substituting Equation (32) into the expression above leads to

$\mathbf{e}_{k+1}=(\mathcal{G}^{{}^{\prime}T}\mathcal{G}^{\prime}+\alpha
H^{T}H)^{-1}\mathcal{G}^{{}^{\prime}T}\left[\mathcal{G}^{\prime}\mathbf{e}_{k}-\mathcal{G}(\mathbf{u}^{(k)})+\mathcal{G}(\mathbf{u}^{\dagger})+A(A^{T}A)^{-1}A^{T}\left(\mathcal{G}(\mathbf{u}^{(k)})-\mathcal{G}(\mathbf{u}^{\dagger})\right)\right].$

Taking norm of both sides, using assumptions (18)–(20) and Lemma 3.3, we get

$\|\mathbf{e}_{k+1}\|\leq\frac{1}{2L_{3}c}\left(\frac{L_{3}}{2}\|\mathbf{e}_{k}\|^{2}+L_{0}\|A(A^{T}A)^{-1}A^{T}\|\|\mathbf{e}_{k}\|\right),$

which we rewrite into the following expression

$\|\mathbf{e}_{k+1}\|\leq\frac{1}{4c}\left(\|\mathbf{e}_{k}\|+b\right)^{2}\leq\frac{1}{2c}\|\mathbf{e}_{k}\|^{2}+\frac{1}{2c}b^{2},$

where $b=\|A(A^{T}A)^{-1}A^{T}\|L_{0}/L_{3}$. Next, we recursively obtain the
following inequality

$\|\mathbf{e}_{k}\|\leq
2^{-k}\left(\frac{1}{c}\right)^{2^{k}-1}\|\mathbf{e}_{0}\|^{2^{k}}+\frac{c}{2}\sum_{i=1}^{k}\left(\frac{b}{c}\right)^{2^{i}},\quad\mbox{for}\quad
k=1,2,\dots$

Since by assumption $b/c<1$ and $\|e_{0}\|<c$, we have that

$\|\mathbf{e}_{k}\|<2^{-k}c+\frac{b}{2}\sum_{i=0}^{k}\left(\frac{b}{c}\right)^{i},$

and in the limit of $k$ goes to infinity

$\limsup_{k\to\infty}\|\mathbf{e}_{k}\|<\frac{b/2}{1-b/c}.$

Furthermore, taking norm of Equation (32) leads to

$\|\mathbf{E}_{k}\|\leq L_{0}\|(A^{T}A)^{-1}A^{T}\|\|\mathbf{e}_{k}\|.$

By taking limit of both sides of the above inequality, as $k\to\infty$, we
conclude that

$\displaystyle\limsup_{k\to\infty}\|\mathbf{E}_{k}\|<L_{0}\|(A^{T}A)^{-1}A^{T}\|\frac{b/2}{1-b/c}.$

∎

## References

  * [1] A. H. Jazwinski, Stochastic processes and filtering theory, Mathematics in science and engineering, Academic press, New York, 1970, uKM.
  * [2] J. M. Lewis, J. C. Derber, The use of adjoint equations to solve a variational adjustment problem with advective constraints, Tellus A (1985).
  * [3] O. Talagrand, Assimilation of observations, an introduction (gtspecial issueltdata assimilation in meteology and oceanography: Theory and practice), Journal of the Meteorological Society of Japan. Ser. II (1997).
  * [4] O. Talagrand, P. Courtier, Variational assimilation of meteorological observations with the adjoint vorticity equation. i: Theory, Quarterly Journal of the Royal Meteorological Society (1987).
  * [5] Y. Sasaki, Some basic formalisms in numerical variational analysis, Monthly Weather Review (1970).
  * [6] Y. Tr’emolet, Accounting for an imperfect model in 4d-var, Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography 132 (621) (2006) 2483–2504.
  * [7] Y. Trémolet, Model-error estimation in 4d-var, Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography 133 (626) (2007) 1267–1280.
  * [8] C. Johnson, N. Nichols, B. Hoskins, Very large inverse problems in atmosphere and ocean modelling, International journal for numerical methods in fluids 47 (8-9) (2005) 759–771.
  * [9] A. Lorenc, S. Ballard, R. Bell, N. Ingleby, P. Andrews, D. Barker, J. Bray, A. Clayton, T. Dalby, D. Li, et al., The met. office global three-dimensional variational data assimilation scheme, Quarterly Journal of the Royal Meteorological Society 126 (570) (2000) 2991–3012.
  * [10] N. Nichols, Data assimilation: Aims and basic concepts, in: Data Assimilation for the Earth System, Springer, 2003, pp. 9–20.
  * [11] H. W. Engl, Inverse problems and their regularization, Computational mathematics driven by industrial problems (2000).
  * [12] J. Kaipio, E. Somersalo, Statistical and computational inverse problems, Springer Science & Business Media, 2006.
  * [13] C. E. Brett, K. F. Lam, K. Law, D. McCormick, M. R. Scott, A. Stuart, Accuracy and stability of filters for dissipative pdes, Physica D: Nonlinear Phenomena 245 (1) (2013) 34–45.
  * [14] K. Hayden, E. Olson, E. S. Titi, Discrete data assimilation in the lorenz and 2d navier–stokes equations, Physica D: Nonlinear Phenomena 240 (18) (2011) 1416–1425.
  * [15] A. Moodey, A. Lawless, R. Potthast, P. Van Leeuwen, Nonlinear error dynamics for cycled data assimilation methods, Inverse Problems 29 (2) (2013) 25002.  doi:https://dx.doi.org/10.1088/0266-5611/29/2/025002.
  * [16] S. Gratton, A. S. Lawless, N. K. Nichols, Approximate gauss–newton methods for nonlinear least squares problems, SIAM Journal on Optimization 18 (1) (2007) 106–132.
  * [17] C. Cartis, M. H. Kaouri, A. S. Lawless, N. K. Nichols, Convergent least-squares optimisation methods for variational data assimilation, arXiv preprint arXiv:2107.12361 (2021).
  * [18] P. Courtier, J.-N. Thépaut, A. Hollingsworth, A strategy for operational implementation of 4d-var, using an incremental approach, Quarterly Journal of the Royal Meteorological Society 120 (519) (1994) 1367–1387.
  * [19] A. Lawless, S. Gratton, N. Nichols, An investigation of incremental 4d-var using non-tangent linear models, Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography 131 (606) (2005) 459–476.
  * [20] A. Lawless, S. Gratton, N. Nichols, Approximate iterative methods for variational data assimilation, International journal for numerical methods in fluids 47 (10-11) (2005) 1129–1135.
  * [21] J. E. Dennis Jr, R. B. Schnabel, Numerical methods for unconstrained optimization and nonlinear equations, SIAM, 1996.
  * [22] E. N. Lorenz, Deterministic Nonperiodic Flow, Journal of Atmospheric Sciences 20 (1963) 130–148.  doi:10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2.
  * [23] E. N. Lorenz, Predictability - a problem partly solved, in: T. Palmer, R. Hagedorn (Eds.), Proceedings of seminar on Predictability, Vol. 1, ECMWF, Cambridge University Press, Reading, UK, 1996, pp. 1–18.