# A Tutorial on Bayesian Data Assimilation

Colin Grudzien  Center for Western Weather and Water Extremes (CW3E), Scripps
Institution of Oceanography, University of California San Diego, San Diego,
CA, USA  Department of Mathematics and Statistics, University of Nevada, Reno,
Reno, NV, USA  Marc Bocquet  CEREA, École des Ponts and EDF R&D, Île-de-
France, France

## Abstract

This tutorial provides a broad introduction to Bayesian data assimilation that
will be useful to practitioners, in interpreting algorithms and results, and
for theoretical studies developing novel schemes with an understanding of the
rich history of geophysical data assimilation and its current directions. The
simple case of data assimilation in a ‘perfect’ model is primarily discussed
for pedagogical purposes. Some mathematical results are derived at a high-
level in order to illustrate key ideas about different estimators. However,
the focus of this work is on the intuition behind these methods, where more
formal and detailed treatments of the data assimilation problem can be found
in the various references. In surveying a variety of widely used data
assimilation schemes, the key message of this tutorial is how the Bayesian
analysis provides a consistent framework for the estimation problem and how
this allows one to formulate its solution in a variety of ways to exploit the
operational challenges in the geosciences.

##  3.1 Introduction

In applications such as short- to medium-range weather prediction, _data
assimilation_ (DA) provides a means to sequentially and recursively update
forecasts of a time-varying physical process with newly incoming information
[26, 51, 4], typically Earth observations. The Bayesian approach to DA is
widely adopted [59] because it provides a unified treatment of tools from
statistical estimation, nonlinear optimisation and machine learning for
handling such a problem. This tutorial illustrates how this approach can be
utilised to develop and interpret a variety of widely used DA algorithms, both
classical and those at the current state-of-the-art.

Suppose that the time-dependent physical states to be modelled can be written
as a vector, $\boldsymbol{x}_{k}\in\mathbb{R}^{N_{x}}$, where $k$ labels some
time $t_{k}$. Formally, the time-evolution of these states is represented with
the nonlinear map $\mathcal{M}$,

$\displaystyle\boldsymbol{x}_{k}=\mathcal{M}_{k}\left(\boldsymbol{x}_{k-1},\boldsymbol{\lambda}\right)+\boldsymbol{\eta}_{k}$
(3.1)

where: (i) $\boldsymbol{x}_{k-1}$ is the state variable vector at an earlier
time $t_{k-1}$; (ii) $\boldsymbol{\lambda}$ is a vector of uncertain static
physical parameters but on which the time evolution depends; and (iii)
$\boldsymbol{\eta}_{k}$ is an additive (for simplicity but other choices are
possible), stochastic noise term, representing errors in our model for the
physical process. Define $\Delta_{t}:=t_{k}-t_{k-1}$ to be a fixed-length
forecast horizon in this tutorial, though none of the following results
require this to be fixed in practice.

The basic goal of sequential DA is to estimate the random state vector
$\boldsymbol{x}_{k}$, given a prior distribution on
$\left(\boldsymbol{x}_{k-1},\boldsymbol{\lambda}\right)$, and given knowledge
of $\mathcal{M}_{k}$ and knowledge of how $\boldsymbol{\eta}_{k}$ is
statistically distributed. At time $t_{k-1}$, a forecast is made for the
distribution of $\boldsymbol{x}_{k}$ utilising the prior knowledge, which
includes the physics-based model. For simplicity, most of this tutorial is
restricted to the case where $\boldsymbol{\lambda}$ is a known constant, and
the forecast model is _perfect_ , i.e.,

$\displaystyle\boldsymbol{x}_{k}=\mathcal{M}_{k}\left(\boldsymbol{x}_{k-1}\right).$
(3.2)

However, a strength of Bayesian analysis is how it easily extends to include a
general treatment of model errors and the estimation of model parameters –
see, e.g., [4] for a more general introduction.

While a forecast is made, one acquires a collection of observations of the
real-world process. This is written as the observation vector
$\boldsymbol{y}_{k}\in\mathbb{R}^{N_{y}}$, which is related to the state
vector by

$\displaystyle\boldsymbol{y}_{k}=\mathcal{H}_{k}\left(\boldsymbol{x}_{k}\right)+\boldsymbol{\epsilon}_{k}.$
(3.3)

The (possibly nonlinear) map
$\mathcal{H}_{k}:\mathbb{R}^{N_{x}}\rightarrow\mathbb{R}^{N_{y}}$ relates the
physical states being modelled, $\boldsymbol{x}_{k}$, to the values that are
actually observed, $\boldsymbol{y}_{k}$. Typically, in geophysical
applications, observations are not $1:1$ with the state variables; while the
data dimension can be extremely large, $N_{y}\ll N_{x}$ so this information is
sparse relative to the model state dimension. The term
$\boldsymbol{\epsilon}_{k}$ in Eq. (3.3) is an additive, stochastic noise term
representing errors in the measurements, or a mismatch between the state
variable representation and the observation [[, the representation
error,]]janjic2018.

Therefore, some time after the real-life physical system has reached time
$t_{k}$, one has a forecast distribution for the states $\boldsymbol{x}_{k}$,
generated by the prior on $\boldsymbol{x}_{k-1}$ and the physics-based model
$\mathcal{M}$, and the observations $\boldsymbol{y}_{k}$ with some associated
uncertainty. The goal of Bayesian DA is to estimate the posterior distribution
for $\boldsymbol{x}_{k}$ conditioned on $\boldsymbol{y}_{k}$, or some
statistics of this distribution.

##  3.2 Hidden Markov Models and Bayesian Analysis

###  3.2.1 The Observation-Analysis-Forecast Cycle

Recursive estimation of the distribution for $\boldsymbol{x}_{k}$ conditional
on $\boldsymbol{y}_{k}$ (i.e., assuming $\boldsymbol{y}_{k}$ is known) can be
described as an observation-analysis-forecast cycle [89]. Given the forecast-
prior for the model state, and the likelihood function for the observation,
Bayes’ law updates the prior for the modelled state to the posterior
conditioned on the observation. Bayes’ law is a simple re-arrangement of
_conditional probability_ , defined by [53] as

$\displaystyle\mathcal{P}(A|B):=\frac{\mathcal{P}(A,B)}{\mathcal{P}(B)},$
(3.4)

for two events $A$ and $B$ where the probability of $B$ is non-zero.
Intuitively, this says that the probability of an event $A$, given knowledge
that an event $B$ occurs, is equal to the probability of $A$ occurring
relative to a sample space restricted to the event $B$. Using the symmetry in
the joint event $\mathcal{P}(A,B)=\mathcal{P}(B,A)$, Bayes’ law is written

$\displaystyle\mathcal{P}(A|B)=\frac{\mathcal{P}(B|A)\mathcal{P}(A)}{\mathcal{P}(B)}.$
(3.5)

In the observation-analysis-forecast cycle, $A$ is identified with the state
vector (seen as a random vector) taking its value in a neighbourhood of
$\boldsymbol{x}_{k}$, $B$ is identified with the observation vector (seen as a
random vector) taking its value in a neighbourhood of $\boldsymbol{y}_{k}$.
The power of this statement is in how it describes an ‘inverse’ probability –
while the _posterior_ , $\mathcal{P}(A|B)$, on the left-hand-side may not be
directly accessible, often the _likelihood_ $\mathcal{P}(B|A)$ and the _prior_
$\mathcal{P}(A)$ are easy to compute, and this is sufficient to develop a
variety of probabilistic DA techniques.

Figure 3.1: Conceptual diagram of the observation-analysis-forecast cycle. The
y-axis represents a variable of the state and observation vectors; the x-axis
represents time. Ellipses represent the spread of the observation / posterior
/ forecast errors. Original figure, adapted from [18].

A conceptual diagram of this process is pictured in Fig. 3.1. Given the
initial first prior (represented as the first ‘posterior’ at time $t_{0}$), a
forecast-prior for the model state at $t_{1}$, $\mathcal{P}(A)$, is produced
with the numerical model. At the update time, there is an (possibly indirect
and noisy) observation of the physical state with an associated likelihood
$\mathcal{P}(B|A)$. The posterior for the model state conditioned on this
observation $\mathcal{P}(A|B)$, commonly denoted the _analysis_ in geophysical
DA, is used to initialise the subsequent numerical forecast. Recursive
estimates of the current modelled state can be performed in this fashion, but
a related question regards the past states. The newly received observation
gives information about the model states at past times and this allows one to
produce a retrospective posterior estimate for the past states. Recursive
estimation of the present state using these incoming observations is commonly
known as _filtering_. Conditional estimation of a past state given a time
series including future observations is commonly known as _smoothing_.

It is important to recognise that the filtering _probability density function_
(pdf) for the current time
$p(\boldsymbol{x}_{k}|\boldsymbol{y}_{k},\boldsymbol{y}_{k-1},\boldsymbol{y}_{k-2},\ldots)$
is actually just a marginal of the joint posterior pdf over all states in the
current _data assimilation window_ (DAW), i.e., the window of lagged past and
current states being estimated. In Fig. 3.1, the DAW is the time window
$\\{t_{1},t_{2},t_{3}\\}$. The conditional pdf for the model state at time
$t_{3}$, given observations in the DAW is written in terms of the joint
posterior as

$\displaystyle
p(\boldsymbol{x}_{3}|\boldsymbol{y}_{3},\boldsymbol{y}_{2},\boldsymbol{y}_{1})=\iiint\\!p(\boldsymbol{x}_{3},\boldsymbol{x}_{2},\boldsymbol{x}_{1},\boldsymbol{x}_{0}|\boldsymbol{y}_{3},\boldsymbol{y}_{2},\boldsymbol{y}_{1})\,\mathrm{d}\boldsymbol{x}_{2}\mathrm{d}\boldsymbol{x}_{1}\mathrm{d}\boldsymbol{x}_{0},$
(3.6)

by averaging out the past history of the model state from the joint posterior
in the integrand. A smoothing estimate may be produced in a variety of ways,
exploiting different formulations of the Bayesian problem. One may estimate
only a _marginal pdf_ as on the left-hand-side of Eq. (3.6), or the entire
_joint posterior pdf_ as in the integrand above [3, 22, 24]. The approach
chosen may strongly depend on whether the DAW is static or is advanced in
time. Particularly, if one produces a smoothing estimate for the state at
times $t_{1}$ through $t_{3}$, one may subsequently shift the DAW so that in
the new cycle, the posterior for times $t_{2}$ through $t_{4}$ is estimated –
this type of analysis is known as _fixed-lag smoothing_. This tutorial
considers how one can utilise a Bayesian _maximum a posteriori_ (MAP)
formalism to efficiently solve the filtering and smoothing problems, using the
various tools of statistical estimation, nonlinear optimisation and machine
learning.

###  3.2.2 A Generic Hidden Markov Model

Recall the perfect physical process model, and the noisy observation model,

$\displaystyle\boldsymbol{x}_{k}$
$\displaystyle=\mathcal{M}_{k}\left(\boldsymbol{x}_{k-1}\right),$ (3.7a)
$\displaystyle\boldsymbol{y}_{k}$
$\displaystyle=\mathcal{H}_{k}\left(\boldsymbol{x}_{k}\right)+\boldsymbol{\epsilon}_{k}.$
(3.7b)

Denote the sequence of the process model states and observation model states
between time $t_{k}$ and time $t_{l}$, for $k<l$, as

$\displaystyle\boldsymbol{x}_{l:k}:=\left\\{\boldsymbol{x}_{l},\boldsymbol{x}_{l-1},\cdots,\boldsymbol{x}_{k}\right\\},$
$\displaystyle\boldsymbol{y}_{l:k}:=\left\\{\boldsymbol{y}_{l},\boldsymbol{y}_{l-1},\cdots,\boldsymbol{y}_{k}\right\\}.$
(3.8)

For arbitrary $l\in\mathbb{N}$, assume that the sequence of observation error

$\displaystyle\\{\boldsymbol{\epsilon}_{l},\boldsymbol{\epsilon}_{l-1},\cdots,1\\}$
(3.9)

is independent-in-time (i.e., a white process).

The above formulation is a type of _hidden Markov model_ , where the dynamic
state variables $\boldsymbol{x}_{k}$ are known as the hidden variables because
they are not directly observed. A Markov model is a type of ‘memoryless’
process, described this way because of how the conditional probability for the
state is represented between different times [75, see Chapter 4].
Particularly, if $\boldsymbol{x}_{k:1}$ is a Markov process, the Markov
property is defined as

$\displaystyle p\left(\boldsymbol{x}_{k}|\boldsymbol{x}_{k-1:0}\right)$
$\displaystyle=p\left(\boldsymbol{x}_{k}|\boldsymbol{x}_{k-1}\right).$ (3.10)

The above statement says that, given knowledge of the state
$\boldsymbol{x}_{k-1}$, the conditional probability for $\boldsymbol{x}_{k}$
is independent of the past history of the state before time $t_{k-1}$,
representing the probabilistic analogue of an initial value problem.

Applying the Markov property recursively with the definition of the
conditional pdf yields

$\displaystyle
p\left(\boldsymbol{x}_{L:0}\right)=p(\boldsymbol{x}_{0})\prod_{k=1}^{L}p(\boldsymbol{x}_{k}|\boldsymbol{x}_{k-1}).$
(3.11)

Therefore, the joint pdf for a forecast of the model state can be written as
the product of the first prior at time $t_{0}$, representing the uncertainty
of the data used to initialise the model forecast, and the product of the
_Markov transition pdfs_ , describing the evolution of the state between
discrete times.

With the perfect state model, as in Eq. (3.7a) above, the transition
probability for some subset
$\mathrm{d}\boldsymbol{x}\subset\mathbb{R}^{N_{x}}$ is written

$\displaystyle\mathcal{P}\left(\boldsymbol{x}_{k}\in\mathrm{d}\boldsymbol{x}|\boldsymbol{x}_{k-1}\right)=\delta_{\mathcal{M}_{k}\left(\boldsymbol{x}_{k-1}\right)}(\mathrm{d}\boldsymbol{x}),$
(3.12)

with $\delta_{\boldsymbol{v}}$ referring to the _Dirac measure_ at
$\boldsymbol{v}\in\mathbb{R}^{N_{x}}$. The Dirac measure satisfies

$\displaystyle\int\\!f(\boldsymbol{x})\delta_{\boldsymbol{v}}\left(\mathrm{d}\boldsymbol{x}\right)=f\left(\boldsymbol{v}\right),$
(3.13)

where this is a singular measure, to be understood by the integral equation.
Accordingly, the transition pdf is often written proportional as

$\displaystyle
p(\boldsymbol{x}_{k}|\boldsymbol{x}_{k-1})\propto\delta\left\\{\boldsymbol{x}_{k}-\mathcal{M}_{k}\left(\boldsymbol{x}_{k-1}\right)\right\\}$
(3.14)

where $\delta$ represents the _Dirac distribution_. Heuristically, this is
known as the ‘function’ defined by the property

$\displaystyle\delta_{\epsilon}(\boldsymbol{x})=\begin{cases}\frac{1}{\epsilon}&\boldsymbol{x}\in\left[-\epsilon,+\epsilon\right]\\\
0&\text{else}\end{cases},$
$\displaystyle\delta(\boldsymbol{x})=\lim_{\epsilon\rightarrow
0^{+}}\delta_{\epsilon}(\boldsymbol{x}).$ (3.15)

However, this is just a convenient abuse of notations, as the Dirac measure
does not have a pdf with respect to the standard Lebesgue measure. Rather, the
Dirac distribution is understood through the generalised function theory of
distributions [85, see Section 3.4] as a type of integral kernel satisfying

$\displaystyle\int\\!f(\boldsymbol{x}_{k})\delta\left\\{\boldsymbol{x}_{k}-\mathcal{M}_{k}\left(\boldsymbol{x}_{k-1}\right)\right\\}\mathrm{d}\boldsymbol{x}_{k}=f\left(\mathcal{M}_{k}\left(\boldsymbol{x}_{k-1}\right)\right).$
(3.16)

Equation (3.16) represents the existence and uniqueness of the solution to an
initial value problem in deterministic systems of ordinary and partial
differential equations, where the knowledge of the state
$\boldsymbol{x}_{k-1}$ completely determines the state $\boldsymbol{x}_{k}$
via the perfect forecast model $\mathcal{M}_{k}$. Particularly, there is
probability $1$ of the state following the unique solution to the time
evolution, and probability $0$ of all other outcomes. Each Markov transition
pdf represents the evolution of an initial condition with respect to the
dynamical model, conditional ultimately on an uncertain outcome of the initial
data from the first prior. While the perfect model assumption is used for
simplicity, the decomposition of the forecast pdf above can be derived for
erroneous models under the additional assumption that the model errors are
independent-in-time. Like the decomposition of the forecast pdf, for the
observation likelihood we can write

$\displaystyle
p\left(\boldsymbol{y}_{k}|\boldsymbol{x}_{k},\boldsymbol{y}_{k-1:1}\right)=p\left(\boldsymbol{y}_{k}|\boldsymbol{x}_{k}\right)$
(3.17)

due to the independence assumption on the observation errors, and the
relationship between $\boldsymbol{x}_{k}$ and $\boldsymbol{y}_{k}$. This says
that the knowledge of the physical state $\boldsymbol{x}_{k}$ completely
determines the likelihood of the observation $\boldsymbol{y}_{k}$. Indeed,

$\displaystyle\boldsymbol{\epsilon}_{k}=\boldsymbol{y}_{k}-\mathcal{H}_{k}(\boldsymbol{x}_{k}),$
(3.18)

which follows a known distribution and is independent of the other observation
error outcomes by assumption.

Consider thus how to estimate the filtering pdf
$p\left(\boldsymbol{x}_{k}|\boldsymbol{y}_{k:1}\right)$. Using the definition
of the conditional pdf, one has

$\displaystyle p\left(\boldsymbol{x}_{k}|\boldsymbol{y}_{k:1}\right)$
$\displaystyle=\frac{p\left(\boldsymbol{y}_{k:1},\boldsymbol{x}_{k}\right)}{p\left(\boldsymbol{y}_{k:1}\right)}.$
(3.19)

Rewriting these pdfs as conditional pdfs, and by using the independence
assumption,

$\displaystyle p\left(\boldsymbol{x}_{k}|\boldsymbol{y}_{k:1}\right)$
$\displaystyle=\frac{p\big{(}\boldsymbol{y}_{k},(\boldsymbol{x}_{k},\boldsymbol{y}_{k-1:1})\big{)}}{p\left(\boldsymbol{y}_{k:1}\right)}$
(3.20a)
$\displaystyle=\frac{p\left(\boldsymbol{y}_{k}|\boldsymbol{x}_{k},\boldsymbol{y}_{k-1:1}\right)p\left(\boldsymbol{x}_{k},\boldsymbol{y}_{k-1:1}\right)}{p\left(\boldsymbol{y}_{k:1}\right)}=\frac{p\left(\boldsymbol{y}_{k}|\boldsymbol{x}_{k}\right)p\left(\boldsymbol{x}_{k},\boldsymbol{y}_{k-1:1}\right)}{p\left(\boldsymbol{y}_{k:1}\right)}.$
(3.20b)

Writing the joint pdfs again in terms of conditional pdfs

$\displaystyle p\left(\boldsymbol{x}_{k}|\boldsymbol{y}_{k:1}\right)$
$\displaystyle=\frac{p\left(\boldsymbol{y}_{k}|\boldsymbol{x}_{k}\right)p\left(\boldsymbol{x}_{k}|\boldsymbol{y}_{k-1:1}\right)p\left(\boldsymbol{y}_{k-1:1}\right)}{p\left(\boldsymbol{y}_{k}|\boldsymbol{y}_{k-1:1}\right)p\left(\boldsymbol{y}_{k-1:1}\right)}$
(3.21a)
$\displaystyle=\frac{p\left(\boldsymbol{y}_{k}|\boldsymbol{x}_{k}\right)p\left(\boldsymbol{x}_{k}|\boldsymbol{y}_{k-1:1}\right)}{p\left(\boldsymbol{y}_{k}|\boldsymbol{y}_{k-1:1}\right)}.$
(3.21b)

Now, suppose that the posterior pdf
$p(\boldsymbol{x}_{k-1}|\boldsymbol{y}_{k-1:1})$ at the last observation time
$t_{k-1}$ is already computed – then the model forecast of this pdf is given
by averaging over the state at time $t_{k-1}$ with respect to the Markov
transition pdf,

$\displaystyle p(\boldsymbol{x}_{k}|\boldsymbol{y}_{k-1:1})$
$\displaystyle=\int\\!p(\boldsymbol{x}_{k}|\boldsymbol{x}_{k-1})p(\boldsymbol{x}_{k-1}|\boldsymbol{y}_{k-1:1})\,\mathrm{d}\boldsymbol{x}_{k-1},$
(3.22)

yielding the forecast-prior. The filtering pdf, on the left-hand-side of Eq.
(3.21) is written in terms of: (i) the likelihood of the observed data given
the model forecast, $p\left(\boldsymbol{y}_{k}|\boldsymbol{x}_{k}\right)$;
(ii) the forecast-prior given the last best estimate of the state,
$p\left(\boldsymbol{x}_{k}|\boldsymbol{y}_{k-1:1}\right)$; and (iii) the
marginal of the joint pdf
$p(\boldsymbol{y}_{k},\boldsymbol{x}_{k}|\boldsymbol{y}_{k-1:1})$, integrating
out the hidden variables,

$\displaystyle
p\left(\boldsymbol{y}_{k}|\boldsymbol{y}_{k-1:1}\right)=\int\\!p(\boldsymbol{y}_{k}|\boldsymbol{x}_{k})p(\boldsymbol{x}_{k}|\boldsymbol{y}_{k-1:1})\,\mathrm{d}\boldsymbol{x}_{k}.$
(3.23)

This type of pdf, only depending on the observations, is called an _evidence_
[[, e.g.,]]carrassi2017.

Typically, the pdf in the denominator of Eq. (3.21) is mathematically
intractable. However, the denominator is independent of the hidden variable
$\boldsymbol{x}_{k}$ by construction – the free argument in the pdf on the
left-hand-side is the model state $\boldsymbol{x}_{k}$ and the purpose of the
denominator on the right-hand-side is only to normalise the integral of the
posterior pdf to $1$. Instead, as a proportionality statement,

$\displaystyle p\left(\boldsymbol{x}_{k}|\boldsymbol{y}_{k:1}\right)$
$\displaystyle\propto
p\left(\boldsymbol{y}_{k}|\boldsymbol{x}_{k}\right)p\left(\boldsymbol{x}_{k}|\boldsymbol{y}_{k-1:1}\right)$
(3.24)

one can devise the Bayesian MAP estimate as the choice of
$\overline{\boldsymbol{x}}_{k}$ that maximises the posterior pdf, but written
in terms of the two right-hand-side components in Eq. (3.24). For the purpose
of maximising the posterior pdf, the denominator leads to insignificant
constants that can be discarded. Thus, in order to compute the MAP
sequentially and recursively in time, one can develop a recursion in
proportionality as above. However, note that the evidence can be estimated and
used for other significant purposes still within a Bayesian framework [20].

###  3.2.3 Linear-Gaussian Models

Generally, the filtering pdf
$p\left(\boldsymbol{x}_{k}|\boldsymbol{y}_{k:1}\right)$ has no analytical
solution, i.e., no explicit expression. However, when the models are linear,
i.e. both the state and observation models are written as matrix actions

$\displaystyle\boldsymbol{x}_{k}$
$\displaystyle=\mathbf{M}_{k}\boldsymbol{x}_{k-1},$ (3.25a)
$\displaystyle\boldsymbol{y}_{k}$
$\displaystyle=\mathbf{H}_{k}\boldsymbol{x}_{k}+\boldsymbol{\epsilon}_{k},$
(3.25b)

and the error pdfs are _Gaussian_ :

$\displaystyle p(\boldsymbol{x}_{0})$
$\displaystyle=n\left(\boldsymbol{x}_{0}|\overline{\boldsymbol{x}}_{0},\mathbf{B}_{0}\right),\quad
p\left(\boldsymbol{y}_{k}|\boldsymbol{x}_{k}\right)=n\left(\boldsymbol{y}_{k}|\mathbf{H}_{k}\boldsymbol{x}_{k},\mathbf{R}_{k}\right),$
(3.26) $\displaystyle n(\boldsymbol{z}|\overline{\boldsymbol{z}},\mathbf{B})$
$\displaystyle:=\frac{1}{\sqrt{\left(2\pi\right)^{N_{z}}\det\left(\mathbf{B}\right)}}\exp\left\\{-\frac{1}{2}\left(\overline{\boldsymbol{z}}-\boldsymbol{z}\right)^{\top}\mathbf{B}^{-1}\left(\overline{\boldsymbol{z}}-\boldsymbol{z}\right)\right\\},$
(3.27)

then the forecast and posterior pdfs are Gaussian at all times, and are
parametrised in terms of their mean and covariance. In particular, Gaussian
distributions are closed under affine transformations, i.e., maps of the form

$\displaystyle\boldsymbol{f}(\boldsymbol{x})=\mathbf{A}\boldsymbol{x}+\boldsymbol{b},$
(3.28)

corresponding to a linear transformation when $\boldsymbol{b}$ is a vector of
zeros, i.e., $\boldsymbol{b}=\boldsymbol{0}$ [87, see Theorem 3.3.3]. If
$\boldsymbol{x}$ is distributed with pdf
$p(\boldsymbol{x})=n\left(\boldsymbol{x}|\overline{\boldsymbol{x}},\mathbf{B}\right)$,
then the random vector
$\boldsymbol{y}:=\mathbf{A}\boldsymbol{x}+\boldsymbol{b}$ is distributed with
pdf

$\displaystyle
p(\boldsymbol{y})=n\left(\boldsymbol{y}|\mathbf{A}\overline{\boldsymbol{x}}+\boldsymbol{b},\mathbf{A}\mathbf{B}\mathbf{A}^{\top}\right).$
(3.29)

Suppose that the last analysis pdf is given as

$\displaystyle
p(\boldsymbol{x}_{k-1}|\boldsymbol{y}_{k-1:1})=n\left(\boldsymbol{x}_{k-1}|\overline{\boldsymbol{x}}_{k-1}^{\mathrm{a}},\mathbf{B}_{k-1}^{\mathrm{a}}\right),$
(3.30)

parametrised in terms of the analysis mean
$\overline{\boldsymbol{x}}^{\mathrm{a}}_{k-1}$ and analysis error covariance
$\mathbf{B}^{\mathrm{a}}_{k-1}$. Then, the forecast-prior pdf is written as

$\displaystyle
p(\boldsymbol{x}_{k}|\boldsymbol{y}_{k-1:1})=n\left(\boldsymbol{x}_{k}|\overline{\boldsymbol{x}}^{\mathrm{f}}_{k-1},\mathbf{B}^{\mathrm{f}}_{k}\right),$
(3.31)

where the forecast mean and forecast error covariance are defined by

$\displaystyle\overline{\boldsymbol{x}}_{k}^{\mathrm{f}}:=\mathbf{M}_{k}\overline{\boldsymbol{x}}^{\mathrm{a}}_{k-1},$
(3.32a)
$\displaystyle\mathbf{B}_{k}^{\mathrm{f}}:=\mathbf{M}_{k}\mathbf{B}^{\mathrm{a}}_{k-1}\mathbf{M}_{k}^{\top},$
(3.32b)

respectively.

Similarly, the conditional and marginal distributions of a Gaussian random
vector are also Gaussian and their pdf has an analytical form. Suppose that
$\boldsymbol{z}\in\mathbb{R}^{N_{z}}$ is an arbitrary Gaussian random vector,
partitioned as

$\displaystyle\boldsymbol{z}:=\begin{pmatrix}\boldsymbol{x}\\\
\boldsymbol{y}\end{pmatrix}$ $\displaystyle
p(\boldsymbol{z}):=n\left(\boldsymbol{z}\Big{|}\begin{pmatrix}\overline{\boldsymbol{x}}\\\
\overline{\boldsymbol{y}}\end{pmatrix},\begin{pmatrix}\boldsymbol{\Sigma}_{xx}&\boldsymbol{\Sigma}_{xy}\\\
\boldsymbol{\Sigma}_{yx}&\boldsymbol{\Sigma}_{yy}\end{pmatrix}\right),$ (3.33)

with the dimensions given as $N_{z}=N_{x}+N_{y}$ and

$\displaystyle\boldsymbol{x},\overline{\boldsymbol{x}}\in\mathbb{R}^{N_{x}},$
$\displaystyle\boldsymbol{y},\overline{\boldsymbol{y}}\in\mathbb{R}^{N_{y}},$
$\displaystyle\boldsymbol{\Sigma}_{xx}\in\mathbb{R}^{N_{x}\times N_{x}},$
(3.34a)
$\displaystyle\boldsymbol{\Sigma}_{xy}=\boldsymbol{\Sigma}_{yx}^{\top}\in\mathbb{R}^{N_{x}\times
N_{y}},$ $\displaystyle\boldsymbol{\Sigma}_{yy}\in\mathbb{R}^{N_{y}\times
N_{y}}.$ (3.34b)

Then, the general form of the pdf for $\boldsymbol{x}$ conditioned on the
outcome of $\boldsymbol{y}$ is given by the Gaussian

$\displaystyle
p(\boldsymbol{x}|\boldsymbol{y})=n\left(\boldsymbol{x}|\overline{\boldsymbol{x}}+\boldsymbol{\Sigma}_{xy}\boldsymbol{\Sigma}_{yy}^{-1}\left(\boldsymbol{y}-\overline{\boldsymbol{y}}\right),\boldsymbol{\Sigma}_{xx}-\boldsymbol{\Sigma}_{xy}\boldsymbol{\Sigma}_{yy}^{-1}\boldsymbol{\Sigma}_{yx}\right),$
(3.35)

where the covariance matrix
$\boldsymbol{\Sigma}_{xx}-\boldsymbol{\Sigma}_{xy}\boldsymbol{\Sigma}_{yy}^{-1}\boldsymbol{\Sigma}_{yx}$
is called the _Schur complement_ [[, see, e.g., Theorem 3.3.4
of]]tong2012multivariate. Noting the form of the linear observation model in
Eq. (3.25), and the independence of the observation errors, it is easy to see
that the vector composed of
$\begin{pmatrix}\boldsymbol{x}_{k}^{\top}&\boldsymbol{y}_{k}^{\top}\end{pmatrix}^{\top}$
is jointly Gaussian. Relying on the identifications

$\displaystyle\boldsymbol{\Sigma}_{xy}=\mathbf{B}_{k}^{\mathrm{f}}\mathbf{H}_{k}^{\top},$
$\displaystyle\boldsymbol{\Sigma}_{yy}=\mathbf{H}_{k}\mathbf{B}_{k}^{\mathrm{f}}\mathbf{H}_{k}^{\top}+\mathbf{R}_{k},$
(3.36)

the posterior pdf at time $t_{k}$ is derived as the Gaussian with analysis
mean and analysis error covariance given by

$\displaystyle\overline{\boldsymbol{x}}_{k}^{\mathrm{a}}$
$\displaystyle:=\overline{\boldsymbol{x}}_{k}^{\mathrm{f}}+\mathbf{B}_{k}^{\mathrm{f}}\mathbf{H}_{k}^{\top}\left(\mathbf{H}_{k}\mathbf{B}_{k}^{\mathrm{f}}\mathbf{H}_{k}^{\top}+\mathbf{R}_{k}\right)^{-1}\left(\boldsymbol{y}_{k}-\mathbf{H}_{k}\overline{\boldsymbol{x}}_{k}^{\mathrm{f}}\right),$
(3.37a) $\displaystyle\mathbf{B}_{k}^{\mathrm{a}}$
$\displaystyle:=\mathbf{B}_{k}^{\mathrm{f}}-\mathbf{B}_{k}^{\mathrm{f}}\mathbf{H}_{k}^{\top}\left(\mathbf{H}_{k}\mathbf{B}_{k}^{\mathrm{f}}\mathbf{H}_{k}^{\top}+\mathbf{R}_{k}\right)^{-1}\mathbf{H}_{k}\mathbf{B}_{k}^{\mathrm{f}}.$
(3.37b)

Defining
$\mathbf{K}_{k}:=\mathbf{B}_{k}^{\mathrm{f}}\mathbf{H}_{k}^{\top}\left(\mathbf{H}_{k}\mathbf{B}_{k}^{\mathrm{f}}\mathbf{H}_{k}^{\top}+\mathbf{R}_{k}\right)^{-1}$
as the Kalman gain, Eqs. (3.37) yield the classical _Kalman filter_ (KF)
update, and this derivation inductively defines the forecast and posterior
distribution for $\boldsymbol{x}_{k}$ at all times. The KF is recognised thus
as the parametric representation of a linear-Gaussian hidden Markov model,
providing a recursion on the first two moments for the forecast and posterior.

This analysis extends, as with the classical KF, easily to incorporate
additive model errors as in Eq. (3.1), see e.g., [3]. However, there are many
ways to formulate the KF and there are some drawbacks of the above approach.
Even when the model forecast equations themselves are linear, if they
functionally depend on an uncertain parameter vector, defined
$\mathbf{M}_{k}(\boldsymbol{\lambda})$, the joint estimation problem of
$(\boldsymbol{x}_{k},\boldsymbol{\lambda})$ can become highly nonlinear and an
iterative approach to the joint estimation may be favourable. Therefore, this
tutorial develops the subsequent extensions to the KF with the MAP approach,
which coincides with the development in least-squares and nonlinear
optimisation.

##  3.3 Least-Squares and Nonlinear Optimisation

###  3.3.1 The 3D Cost Function from Gaussian Statistics

As seen in the previous section, the linear-Gaussian analysis admits an
analytical solution for the forecast and posterior pdf at all times. However,
an alternative approach for deriving the KF is formed using the MAP estimation
in proportionality. Consider that the _natural logarithm_ ($\log$) is
monotonic, so that an increase in the input of the argument corresponds
identically to an increase in the output of the function. Therefore,
maximising the posterior pdf as in Eq. (3.24) is equivalent to maximising

$\displaystyle\log\left(p\left(\boldsymbol{y}_{k}|\boldsymbol{x}_{k}\right)p\left(\boldsymbol{x}_{k}|\boldsymbol{y}_{k-1:1}\right)\right)$
$\displaystyle=\log\left(p\left(\boldsymbol{y}_{k}|\boldsymbol{x}_{k}\right)\right)+\log\left(p\left(\boldsymbol{x}_{k}|\boldsymbol{y}_{k-1:1}\right)\right).$
(3.38)

Given the form of the the multivariate Gaussian pdf in Eq. (3.27), maximising
the above log-posterior is equivalent to minimising the following least-
squares cost function, derived in proportionality to the minus-log-posterior:

$\displaystyle\mathcal{J}_{\mathrm{KF}}(\boldsymbol{x}_{k})$
$\displaystyle=\frac{1}{2}\parallel\overline{\boldsymbol{x}}_{k}^{\mathrm{f}}-\boldsymbol{x}_{k}\parallel_{\mathbf{B}_{k}^{\mathrm{f}}}^{2}+\frac{1}{2}\parallel\boldsymbol{y}_{k}-\mathbf{H}_{k}\boldsymbol{x}_{k}\parallel_{\mathbf{R}_{k}}^{2}.$
(3.39)

For an arbitrary positive definite matrix $\mathbf{A}$ the _Mahalanobis
distance_ [61] with respect to $\mathbf{A}$ is defined as

$\displaystyle\parallel\boldsymbol{v}\parallel_{\mathbf{A}}:=\sqrt{\boldsymbol{v}^{\top}\mathbf{A}^{-1}\boldsymbol{v}}.$
(3.40)

The above distances are weighted Euclidean norms with: (i)
$\parallel\circ\parallel_{\mathbf{B}_{k}^{\mathrm{f}}}$ weighting relative to
the forecast spread; and (ii) $\parallel\circ\parallel_{\mathbf{R}_{k}}$
weighting relative to the observation imprecision. The MAP state thus
interpolates the forecast mean and the observation relative to the uncertainty
in each piece of data. Due to the unimodality of the Gaussian, and its
symmetry about its mean, it is clear that the conditional mean
$\overline{\boldsymbol{x}}^{\mathrm{a}}_{k}$ is also the MAP state.

While this cost function analysis provides the solution to finding the first
moment of the Gaussian posterior, this does not yet address how to find the
posterior error covariance. In this linear-Gaussian setting, it is easily
shown that the analysis error covariance is actually given by

$\displaystyle\begin{split}\mathbf{B}_{k}^{\mathrm{a}}:=&\left[\left(\mathbf{B}_{k}^{\mathrm{f}}\right)^{-1}+\mathbf{H}_{k}^{\top}\mathbf{R}_{k}^{-1}\mathbf{H}_{k}\right]^{-1}=\boldsymbol{\Xi}_{\mathcal{J}_{\mathrm{KF}}}^{-1},\end{split}$
(3.41)

where $\boldsymbol{\Xi}_{\mathcal{J}_{\mathrm{KF}}}$ refers to the _Hessian_
of the least-squares cost function in Eq. (3.39), i.e., the matrix of its
mixed second partial derivatives in the state vector variables. This is a
fundamental result that links the recursive analysis in time to the MAP
estimation performed with linear least-squares, see, e.g., Section 6.2 of [73]
for further details.

A more general result from maximum likelihood estimation extends this analysis
as an approximation to nonlinear state and observation models, and to non-
Gaussian error distributions. Define $\hat{\boldsymbol{\theta}}$ to be the
_maximum likelihood estimator_ (MLE) for some unknown parameter vector
$\boldsymbol{\theta}\in\mathbb{R}^{N_{\theta}}$, where
$\hat{\boldsymbol{\theta}}$ depends on the realisation of some arbitrarily
distributed random sample of observed data
$\\{\boldsymbol{z}_{i}\\}_{i=1}^{N}$. Then, under fairly general regularity
conditions, the MLE satisfies

$\displaystyle\mathcal{I}\left(\boldsymbol{\theta}\right)^{-\frac{1}{2}}\left(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}\right)\rightarrow_{d}\mathcal{N}\left(\boldsymbol{0},\mathbf{I}_{N_{\theta}}\right)$
(3.42)

where: (i) Eq. (3.42) refers to convergence in distribution of the random
vector
$\mathcal{I}\left(\boldsymbol{\theta}\right)^{-\frac{1}{2}}\left(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}\right)$
as the sample size $N\rightarrow\infty$; (ii) $\mathbf{I}_{N_{\theta}}$ is the
identity matrix in the dimension of $\boldsymbol{\theta}$; (iii)
$\mathcal{N}(\boldsymbol{0},\mathbf{I}_{N_{\theta}})$ is the multivariate
Gaussian distribution with mean zero and covariance equal to the identity
matrix, i.e., with pdf
$n(\boldsymbol{x}|\boldsymbol{0},\mathbf{I}_{N_{\theta}})$; and (iv)
$\mathcal{I}(\boldsymbol{\theta})$ is the _Fisher information matrix_. The
Fisher information matrix is defined as the expected value of the Hessian of
the minus-log-likelihood, taken over the realisations of the observed data
$\\{\boldsymbol{z}_{i}\\}_{i=1}^{N}$, and with respect to the true parameter
vector $\boldsymbol{\theta}$. It is common to approximate the distribution of
the MLE as

$\displaystyle\hat{\boldsymbol{\theta}}\sim\mathcal{N}\left(\boldsymbol{\theta},I\left(\boldsymbol{\hat{\theta}}\right)\right)$
(3.43)

where $I\left(\hat{\boldsymbol{\theta}}\right)$ is the observed Fisher
information, i.e., the realisation of the Hessian of the minus-log-likelihood
given the observed sample data. This approximation improves in large sample
sizes like the central limit theorem, see, e.g, Chapter 9 of [67] for the
details of the above discussion. In nonlinear optimisation, the approximate
distribution of the MLE therefore relates the geometry in the neighbourhood of
of a local minimiser to the variation in the optimal estimate. Particularly,
the curvature of the cost function level contours around the local minimum is
described by the Hessian while the spread of the distribution of the MLE is
described by the analysis error covariance.

###  3.3.2 3D-VAR and the Extended Kalman Filter

One of the benefits of the above cost function approach is that this
immediately extends to handle a nonlinear observation operator
$\mathcal{H}_{k}$ as a problem of nonlinear least-squares. When the
observations are related nonlinearly to the state vector, the Bayesian
posterior is no longer generally Gaussian, but the interpretation of the
analysis state $\overline{\boldsymbol{x}}^{\mathrm{a}}_{k}$ interpolating
between a background proposal $\overline{\boldsymbol{x}}_{k}^{\mathrm{f}}$ and
the observed data $\boldsymbol{y}_{k}$ relative to their respective
uncertainties remains valid. The use of the nonlinear least-squares cost
function can be considered as making a Gaussian approximation, similar to the
large sample theory in maximum likelihood estimation. However, the forecast
and analysis background error covariances, $\mathbf{B}^{\mathrm{f/a}}_{k}$,
take on different interpretations depending on the DA scheme.

In full-scale geophysical models, the computation and storage of the
background error covariance $\mathbf{B}_{k}^{\mathrm{f/a}}$ is rarely feasible
due to its large size – in practice, this is usually treated abstractly in its
action as an operator by preconditioning the optimisation [82]. One
traditional approach to handle this reduced representation is to use a static-
in-time background for the cost function in Eq. (3.39), rendering the _three-
dimensional, variational_ (3D-VAR) cost function

$\displaystyle\mathcal{J}_{\mathrm{3D}\text{-}\mathrm{VAR}}(\boldsymbol{x}):=\frac{1}{2}\parallel\overline{\boldsymbol{x}}_{k}^{\mathrm{f}}-\boldsymbol{x}\parallel_{\mathbf{B}_{\mathrm{3D}\text{-}\mathrm{VAR}}}^{2}+\frac{1}{2}\parallel\boldsymbol{y}_{k}-\mathcal{H}_{k}\left(\boldsymbol{x}\right)\parallel^{2}_{\mathbf{R}_{k}}.$
(3.44)

The 3D-VAR background error covariance is typically defined as a posterior,
‘climatological’ covariance, taken with respect to a long-time average over
the modelled state. One way that this is roughly estimated in meteorology is
by averaging over a _reanalysis data set_ , which is equivalent to averaging
with respect to a joint smoothing posterior over a long history for the system
[51, and references therein]. If one assumes that the DA cycle represents a
stationary, ergodic process in its posterior statistics, and that this admits
an invariant, ergodic measure, $\pi^{\ast}$, one would define the 3D-VAR
background error covariance matrix with the expected value with respect to
this measure as

$\displaystyle\overline{\boldsymbol{x}}^{\ast}$
$\displaystyle:=\mathbb{E}_{\pi^{\ast}}\left[\boldsymbol{x}\right],$ (3.45a)
$\displaystyle\mathbf{B}_{\mathrm{3D}\text{-}\mathrm{VAR}}$
$\displaystyle:=\mathbb{E}_{\pi^{\ast}}\left[\left(\boldsymbol{x}-\overline{\boldsymbol{x}}^{\ast}\right)\left(\boldsymbol{x}-\overline{\boldsymbol{x}}^{\ast}\right)^{\top}\right].$
(3.45b)

The 3D-VAR cycle is defined where $\overline{\boldsymbol{x}}^{\mathrm{a}}_{k}$
is the _minimising argument (argmin)_ of Eq. (3.44) and subsequently
$\overline{\boldsymbol{x}}_{k+1}^{\mathrm{f}}:=\mathcal{M}_{k+1}\left(\overline{\boldsymbol{x}}^{\mathrm{a}}_{k}\right)$,
i.e., the forecast mean is estimated with the background control trajectory
propagated through the fully nonlinear model. Therefore, at every iteration of
the algorithm, 3D-VAR can be understood to treat the proposal
$\overline{\boldsymbol{x}}^{\mathrm{f}}_{k}$ as a random draw from the
invariant, climatological-posterior measure of the system; the optimised
solution $\overline{\boldsymbol{x}}_{k}^{\mathrm{a}}$ is thus the state that
interpolates the forecast mean and the currently observed data as if the
forecast mean was drawn randomly from the climatological-posterior.

For illustration, consider a cost-effective, explicit formulation of 3D-VAR
where the background error covariance matrix is represented as a Cholesky
factorisation. Here, one performs the optimisation by writing the modelled
state as a perturbation of the forecast mean,

$\displaystyle\boldsymbol{x}_{k}:=\overline{\boldsymbol{x}}_{k}^{\mathrm{f}}+\boldsymbol{\Sigma}\boldsymbol{w},$
(3.46)

where
$\mathbf{B}_{\mathrm{3D}\text{-}\mathrm{VAR}}:=\boldsymbol{\Sigma}\boldsymbol{\Sigma}^{\top}$.
The weight vector $\boldsymbol{w}$ gives the linear combination of the columns
of the matrix factor that describe $\boldsymbol{x}_{k}$ as a perturbation. One
can iteratively optimise the 3D-VAR cost function as such with a locally
quadratic approximation. Let $\overline{\boldsymbol{x}}^{i}_{k}$ be the $i$-th
iterate for the proposal value at time $t_{k}$, defining
$\overline{\boldsymbol{x}}^{0}_{k}:=\overline{\boldsymbol{x}}_{k}^{\mathrm{f}}$.
The $i+1$-st iteration is defined (up to a subtlety pointed out by [62]) as

$\displaystyle\overline{\boldsymbol{x}}^{i+1}_{k}=\overline{\boldsymbol{x}}^{i}_{k}+\boldsymbol{\Sigma}\overline{\boldsymbol{w}}^{i+1},$
(3.47)

where $\overline{\boldsymbol{w}}^{i+1}$ is the argmin of the following
_incremental cost function_

$\displaystyle\mathcal{J}_{\mathrm{3D}\text{-}\mathrm{VARI}}(\boldsymbol{w}):=\frac{1}{2}\parallel\boldsymbol{w}\parallel^{2}+\frac{1}{2}\parallel\boldsymbol{y}_{k}-\mathcal{H}_{k}\left(\overline{\boldsymbol{x}}^{i}_{k}\right)-\mathbf{H}_{k}\boldsymbol{\Sigma}\boldsymbol{w}\parallel^{2}_{\mathbf{R}_{k}}.$
(3.48)

This approximate cost function follows from truncating the Taylor expansion

$\displaystyle\mathcal{H}_{k}(\boldsymbol{x})=\mathcal{H}_{k}\left(\overline{\boldsymbol{x}}^{i}_{k}\right)+\mathbf{H}_{k}\boldsymbol{\Sigma}\boldsymbol{w}+\mathcal{O}\left(\parallel\boldsymbol{w}\parallel^{2}\right)$
(3.49)

and the definition of the Mahalanobis distance. Notice that
$\overline{\boldsymbol{w}}^{i+1}$ is defined uniquely as Eq. (3.48) is
quadratic with respect to the weights, allowing one to use fast and adequate
methods such as the conjugate gradient [66, see Chapter 5]. The iterations are
set to terminate when $\parallel\overline{\boldsymbol{w}}^{i}\parallel$ is
sufficiently small, representing a negligible change from the last proposal.
This approach to nonlinear optimisation is known as the incremental approach,
and corresponds mathematically to the Gauss-Newton optimisation method, an
efficient simplification of Newton’s descent method for nonlinear least-
squares [66, see Chapter 10]. As described above, it consists of an iterative
linearisation of the nonlinear optimisation problem, which is a concept
utilised throughout this tutorial.

The 3D-VAR approach described above played an important role in early DA
methodology, and aspects of this technique are still widely used in ensemble-
based DA when using the regularisation technique known as covariance
hybridisation [44, 60, 68]. However, optimising the state with the
climatological-posterior background neglects important information in the
time-dependent features of the forecast spread, described sometimes as the
‘errors of the day’ [23]. The _extended Kalman filter_ (EKF) can be
interpreted to generalise the 3D-VAR cost function by including a time-
dependent background, like the KF, by approximating its time-evolution at
first order with the tangent linear model as follows.

Suppose that the equations of motion are generated by a nonlinear function,
independent of time for simplicity,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{x}:=\boldsymbol{f}(\boldsymbol{x}),$
$\displaystyle\boldsymbol{x}_{k}=\mathcal{M}_{k}(\boldsymbol{x}_{k-1}):=\int_{t_{k-1}}^{t_{k}}\\!\boldsymbol{f}(\boldsymbol{x})\,\mathrm{d}t+\boldsymbol{x}_{k-1}.$
(3.50)

One can extend the linear-Gaussian approximation for the forecast pdf by
modelling the state as a perturbation of the analysis mean,

$\displaystyle\boldsymbol{x}_{k-1}:=\overline{\boldsymbol{x}}^{\mathrm{a}}_{k-1}+\boldsymbol{\delta}_{k-1}\sim\mathcal{N}\left(\overline{\boldsymbol{x}}^{\mathrm{a}}_{k-1},\mathbf{B}^{\mathrm{a}}_{k-1}\right),$
(3.51a) $\displaystyle\Leftrightarrow\,\,\,$
$\displaystyle\boldsymbol{\delta}_{k-1}\sim\mathcal{N}\left(\boldsymbol{0},\mathbf{B}^{\mathrm{a}}_{k-1}\right).$
(3.51b)

The evolution of the perturbation $\boldsymbol{\delta}_{k-1}$ is written via
Taylor’s theorem as

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\delta}_{k-1}$
$\displaystyle:=\frac{\mathrm{d}}{\mathrm{d}t}\left(\boldsymbol{x}_{k-1}-\overline{\boldsymbol{x}}_{k-1}^{\mathrm{a}}\right)$
(3.52a)
$\displaystyle=\boldsymbol{f}\left(\boldsymbol{x}_{k-1}\right)-\boldsymbol{f}\left(\overline{\boldsymbol{x}}_{k-1}^{\mathrm{a}}\right)$
(3.52b)
$\displaystyle=\nabla_{\boldsymbol{x}}\boldsymbol{f}\left(\overline{\boldsymbol{x}}_{k-1}^{\mathrm{a}}\right)\boldsymbol{\delta}_{k-1}+\mathcal{O}\left(\parallel\boldsymbol{\delta}_{k-1}\parallel^{2}\right),$
(3.52c)

where
$\nabla_{\boldsymbol{x}}\boldsymbol{f}(\overline{\boldsymbol{x}}_{k-1}^{\mathrm{a}})$
is the Jacobian equation with dependence on the underlying analysis mean
state. The linear evolution defined by the truncated Taylor expansion about
the underlying reference solution $\overline{\boldsymbol{x}}$,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\delta}:=\nabla_{\boldsymbol{x}}\boldsymbol{f}(\overline{\boldsymbol{x}})\cdot\boldsymbol{\delta},$
(3.53)

is known as the _tangent linear model_.

Making the approximation of the tangent linear model for the first order
evolution of the modelled state,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{x}\approx\boldsymbol{f}(\overline{\boldsymbol{x}})+\nabla_{\boldsymbol{x}}\boldsymbol{f}(\overline{\boldsymbol{x}})\cdot\boldsymbol{\delta}$
(3.54a) $\displaystyle\Rightarrow$
$\displaystyle\int_{t_{k-1}}^{t_{k}}\\!\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{x}\,\mathrm{d}t\approx\int_{t_{k-1}}^{t_{k}}\\!\boldsymbol{f}(\overline{\boldsymbol{x}})\,\mathrm{d}t+\int_{t_{k-1}}^{t_{k}}\\!\nabla_{\boldsymbol{x}}\boldsymbol{f}(\overline{\boldsymbol{x}})\cdot\boldsymbol{\delta}\,\mathrm{d}t$
(3.54b) $\displaystyle\Rightarrow$
$\displaystyle\boldsymbol{x}_{k}\approx\mathcal{M}_{k}\left(\overline{\boldsymbol{x}}_{k-1}\right)+\mathbf{M}_{k}\boldsymbol{\delta}_{k-1}$
(3.54c)

where $\mathbf{M}_{k}$ is the resolvent of the tangent linear model. Given
that Gaussians are closed under affine transformations, the EKF approximation
for the (perfect) evolution of the state vector is defined as

$\displaystyle p(\boldsymbol{x}_{k}|\boldsymbol{y}_{k-1:1})\approx
n\left(\boldsymbol{x}_{k}|\mathcal{M}_{k}\left(\overline{\boldsymbol{x}}^{\mathrm{a}}_{k-1}\right),\mathbf{M}_{k}\mathbf{B}^{\mathrm{a}}_{k-1}\mathbf{M}^{\top}_{k}\right).$
(3.55)

Respectively, define the EKF analysis cost function as

$\displaystyle\mathcal{J}_{\mathrm{EKF}}(\boldsymbol{x}):=\frac{1}{2}\parallel\overline{\boldsymbol{x}}_{k}^{\mathrm{f}}-\boldsymbol{x}\parallel_{\mathbf{B}^{\mathrm{f}}_{k}}^{2}+\frac{1}{2}\parallel\boldsymbol{y}_{k}-\mathcal{H}_{k}\left(\boldsymbol{x}\right)\parallel^{2}_{\mathbf{R}_{k}},$
(3.56)

where

$\displaystyle\overline{\boldsymbol{x}}_{k}^{\mathrm{f}}$
$\displaystyle:=\mathcal{M}_{k}\left(\overline{\boldsymbol{x}}^{\mathrm{a}}_{k-1}\right),$
(3.57a) $\displaystyle\mathbf{B}_{k}^{\mathrm{f}}$
$\displaystyle:=\mathbf{M}_{k}\mathbf{B}^{\mathrm{a}}_{k-1}\mathbf{M}^{\top}_{k}.$
(3.57b)

A similar, locally quadratic, weight-space optimisation can be performed
versus the nonlinear observation operator using a matrix decomposition for the
time-varying background forecast error covariance, defining the state as a
perturbation using this matrix factor similar to Eq. (3.46). This gives the
square root EKF, which was used historically to improve the stability over the
direct approach [86].

The accuracy and stability of the above EKF depends strongly on the length of
the forecast horizon $\Delta_{t}$ [64]. For short-range forecasting, the
perturbation dynamics of the tangent linear model can be an adequate
approximation, and is an underlying approximation for most operational DA
[18]. However, the explicit linear approximation can degrade quickly, and
especially when the mean state is not accurately known. The perturbation
dynamics in the linearisation about the approximate mean can differ
substantially from the nonlinear dynamics of the true system as the
approximate mean state diverges from the modelled system [38, see Chapters 7 -
8]. These stability issues, and the computational cost of explicitly
representing the evolution of the background covariance in the tangent linear
model, limits the use of the EKF for geophysical DA. Nonetheless, this
technique provides important intuition that is developed later in this
tutorial.

###  3.3.3 The Gaussian 4D Cost Function

Consider again the perfect, linear-Gaussian model represented by Eq. (3.25).
Using the Markov assumption, Eq. (3.10), and the independence of observation
errors, Eq. (3.17), recursively for the hidden Markov model, the joint
posterior decomposes as

$\displaystyle p(\boldsymbol{x}_{L:0}|\boldsymbol{y}_{L:1})$
$\displaystyle\propto\underbrace{p(\boldsymbol{x}_{0})}_{(i)}\underbrace{\left[\prod_{k=1}^{L}p(\boldsymbol{x}_{k}|\boldsymbol{x}_{k-1})\right]}_{(ii)}\underbrace{\left[\prod_{k=1}^{L}p(\boldsymbol{y}_{k}|\boldsymbol{x}_{k})\right]}_{(iii)}$
(3.58)

where (i) is the prior for the initial state $\boldsymbol{x}_{0}$; (ii) is the
free-forecast with the perfect model $\mathbf{M}_{k}$ depending on some
outcome drawn from the prior; and (iii) is the joint likelihood of the
observations in the DAW, given the background forecast.

Define the composition of the linear model forecast from time $t_{k-1}$ to
$t_{l}$ as

$\displaystyle\mathbf{M}_{l:k}:=\mathbf{M}_{l}\cdots\mathbf{M}_{k},$
$\displaystyle\mathbf{M}_{k:k}:=\mathbf{I}_{N_{x}}.$ (3.59)

Using the perfect, linear model hypothesis, note that

$\displaystyle\boldsymbol{x}_{k}:=\mathbf{M}_{k:1}\boldsymbol{x}_{0}$ (3.60)

for every $k$. Therefore, the transition pdfs in Eq. (3.58) are reduced to a
trivial condition by re-writing

$\displaystyle p(\boldsymbol{x}_{L:0}|\boldsymbol{y}_{L:1})$
$\displaystyle\propto
p(\boldsymbol{x}_{0})\left[\prod_{k=1}^{L}p(\boldsymbol{x}_{k}|\mathbf{M}_{k-1:1}\boldsymbol{x}_{0})\right]\left[\prod_{k=1}^{L}p(\boldsymbol{y}_{k}|\boldsymbol{x}_{k})\right]$
(3.61a) $\displaystyle\propto
p(\boldsymbol{x}_{0})\left[\prod_{k=1}^{L}p(\boldsymbol{y}_{k}|\mathbf{M}_{k:1}\boldsymbol{x}_{0})\right]$
(3.61b)

as this pdf evaluates to zero whenever
$\boldsymbol{x}_{k}\neq\mathbf{M}_{k:1}\boldsymbol{x}_{0}$. Given a Gaussian
prior, and Gaussian observation error distributions, the minus-log-posterior
_four-dimensional_ (4D) cost function is derived as

$\displaystyle\mathcal{J}_{\mathrm{4D}}(\boldsymbol{x}_{0}):=\frac{1}{2}\parallel\overline{\boldsymbol{x}}_{0}-\boldsymbol{x}_{0}\parallel^{2}_{\mathbf{B}_{0}}+\frac{1}{2}\sum_{k=1}^{L}\parallel\boldsymbol{y}_{k}-\mathbf{H}_{k}\mathbf{M}_{k:1}\boldsymbol{x}_{0}\parallel_{\mathbf{R}_{k}}^{2}.$
(3.62)

Notice that the 4D cost function in Eq. (3.62) is actually quadratic with
respect to the initial condition $\boldsymbol{x}_{0}$. Therefore, the
smoothing problem in the perfect, linear-Gaussian model has a unique, optimal
initial condition that minimises the sum-of-square deviations from the prior
mean, with distance weighted with respect to the prior covariance, and the
observations, with distance weighted with respect to the observation error
covariances. The optimal, smoothed initial condition using observation
information up to time $t_{L}$ is denoted
$\overline{\boldsymbol{x}}^{\mathrm{s}}_{0|L}$; this gives the smoothed model
states at subsequent times by the perfect model evolution,

$\displaystyle\overline{\boldsymbol{x}}^{\mathrm{s}}_{k|L}:=\mathbf{M}_{k:1}\overline{\boldsymbol{x}}^{\mathrm{s}}_{0|L}.$
(3.63)

The derivation above is formulated for a smoothing problem with a static-in-
time DAW, in which an entire time series of observations
$\boldsymbol{y}_{L:1}$ is available, and for which the estimation may be
performed ‘offline’. However, a simple extension of this analysis accommodates
DAWs that are sequentially shifted in time, allowing an ‘online’ smoothing
estimate to be formed analogously to sequential filtering.

Fixed-lag smoothing sets the length of the DAW, $L$, to be fixed for all time,
while the underlying states are cycled through this window as it shifts
forwards in time. Given the _lag_ of length $L$, suppose that a _shift_ $S$ is
defined for which $1\leq S\leq L$. It is convenient to consider an
algorithmically stationary DAW, referring to the time indices
$\\{t_{1},\cdots,t_{L}\\}$. In a given cycle, the joint posterior
$p(\boldsymbol{x}_{L:1}|\boldsymbol{y}_{L:1})$ is estimated. After the
estimate is produced, the DAW is subsequently shifted in time by
$S\times\Delta_{t}$, and all states are re-indexed by $t_{k}\longleftarrow
t_{k+S}$ to begin the next cycle. For a lag of $L$ and a shift of $S$, the
observation vectors at times $\\{t_{L-S+1},\cdots,t_{L}\\}$ correspond to
observations newly entering the DAW for the analysis performed up to time
$t_{L}$. When $S=L$, the DAWs are disconnected and adjacent in time, whereas
for $S<L$ there is an overlap between the estimated states in sequential DAWs.
Figure 3.2 provides a schematic of how the DAW is shifted for a lag of $L=5$
and shift $S=2$. Following the common convention in DA that there is no
observation at time zero, in addition to the DAW, $\\{t_{1},\cdots,t_{L}\\}$,
states at time $t_{0}$ may be estimated or utilised in order to connect
estimates between adjacent / overlapping DAWs.

Figure 3.2: Three cycles of a shift $S=2$, lag $L=5$ fixed-lag smoother, cycle
number is increasing top to bottom. Time indices on the left-hand margin
indicate the current time for the associated cycle of the algorithm. New
observations entering the current DAW are shaded black. Source: [41], adapted
from [4].

Consider the algorithmically stationary DAW, and suppose that the current
analysis time is $t_{L}$, where the joint posterior pdf
$p(\boldsymbol{x}_{L-S:1-S}|\boldsymbol{y}_{L-S:1-S})$ is available from the
last fixed-lag smoothing cycle at analysis time $t_{L-S}$. Using the
independence of observation errors and the Markov assumption recursively,

$\displaystyle\begin{split}p(\boldsymbol{x}_{L:1}&|\boldsymbol{y}_{L:1-S})\propto\\\
&\int\\!\underbrace{\mathrm{d}\boldsymbol{x}_{0}\,p\left(\boldsymbol{x}_{0}|\boldsymbol{y}_{L-S:1-S}\right)}_{(i)}\underbrace{\left[\prod_{k=1}^{L}p(\boldsymbol{x}_{k}|\boldsymbol{x}_{k-1})\right]}_{(ii)}\underbrace{\left[\prod_{k=L-S+1}^{L}p(\boldsymbol{y}_{k}|\boldsymbol{x}_{k})\right]}_{(iii)},\end{split}$
(3.64)

where: (i) now represents averaging out the initial condition at time $t_{0}$
with respect to the marginal smoothing pdf for
$\boldsymbol{x}^{\mathrm{s}}_{0|L-S}$ over the last DAW; (ii) represents the
free forecast of the smoothed estimate for
$\boldsymbol{x}^{\mathrm{s}}_{0|L-S}$; and term (iii) represents the joint
likelihood of the newly incoming observations to the DAW given the forecasted
model state. Noting that $p(\boldsymbol{x}_{L:1}|\boldsymbol{y}_{L:1})\propto
p(\boldsymbol{x}_{L:1}|\boldsymbol{y}_{L:1-S})$, this provides a recursive
form of the 4D cost function, shifted sequentially in time,

$\displaystyle\mathcal{J}_{\mathrm{4D}\text{-}\mathrm{seq}}(\boldsymbol{x}_{0}):=\frac{1}{2}\parallel\overline{\boldsymbol{x}}_{0|L-S}^{\mathrm{s}}-\boldsymbol{x}_{0}\parallel^{2}_{\mathbf{B}_{0|L-S}^{\mathrm{s}}}+\frac{1}{2}\sum_{k=L-S+1}^{L}\parallel\boldsymbol{y}_{k}-\mathbf{H}_{k}\mathbf{M}_{k:1}\boldsymbol{x}_{0}\parallel_{\mathbf{R}_{k}}^{2},$
(3.65)

where $\mathbf{B}_{0|L-S}^{\mathrm{s}}$ refers to the smoothed error
covariance from the last DA cycle. As with the KF cost function, the posterior
error covariance, conditioning on observations up to time $t_{L}$, is
identified with

$\displaystyle\mathbf{B}^{\mathrm{s}}_{0|L}:=\boldsymbol{\Xi}_{\mathcal{J}_{\mathrm{4D}\text{-}\mathrm{seq}}}^{-1}.$
(3.66)

Given the perfect, linear-Gaussian model assumption, the mean and covariance
are propagated to time $t_{S}$ via

$\displaystyle\overline{\boldsymbol{x}}^{\mathrm{s}}_{S|L}$
$\displaystyle:=\mathbf{M}_{S:1}\overline{\boldsymbol{x}}_{0|L}^{\mathrm{s}},$
(3.67a) $\displaystyle\mathbf{B}_{S|L}^{\mathrm{s}}$
$\displaystyle:=\mathbf{M}_{S:1}\mathbf{B}_{0|L}^{\mathrm{s}}\mathbf{M}_{S:1}^{\top},$
(3.67b)

and states are re-indexed as $t_{k+S}\longleftarrow t_{k}$ to initialise the
next cycle. This provides a 4D derivation of the _Kalman smoother_ (KS)
assuming a perfect model, though not all developments take the above approach,
using a global analysis over all observations in the current DAW at once.
Other formulations use an alternating: (i) forwards sequential filtering pass
to update the current forecast; and (ii) a backwards-in-time sequential
filtering pass over the DAW to condition lagged states on the new information.
Examples of smoothers that follow this alternating analysis include the
_ensemble Kalman smoother_ (EnKS) [32] and the ensemble _Rauch-Tung-Striebel_
(RTS) smoother [70], with the EnKS to be discussed in Section 3.4.3.

###  3.3.4 Incremental 4D-VAR

The method of incremental _four-dimensional, variational_ (4D-VAR) data
assimilation [57, 83, 25] is a classical and widely used DA technique that
extends the linear-Gaussian 4D analysis to nonlinear settings, both for fixed
and sequential DAWs. Modern formulations of the 4D-VAR analysis furthermore
handle model errors, as in weak-constraint 4D-VAR [88, 29, 55], and may
include time-varying background error covariances by combining an ensemble of
4D-VAR [15], which is known as _ensemble of data assimilation_ (EDA).

Assuming now that the state and observation models are nonlinear, as in Eq.
(3.7), denote the composition of the nonlinear model forecast from time
$t_{k-1}$ to time $t_{l}$ as

$\displaystyle\mathcal{M}_{l:k}:=\mathcal{M}_{l}\circ\cdots\circ\mathcal{M}_{k},$
$\displaystyle\mathcal{M}_{k:k}:=\mathbf{I}_{N_{x}}.$ (3.68)

Then, the linear-Gaussian 4D cost function is formally extended as a Gaussian
approximation in 4D-VAR with

$\displaystyle\mathcal{J}_{\mathrm{4D}\text{-}\mathrm{VAR}}(\boldsymbol{x}_{0}):=\frac{1}{2}\parallel\overline{\boldsymbol{x}}_{0}-\boldsymbol{x}_{0}\parallel^{2}_{\mathbf{B}_{\mathrm{4D}\text{-}\mathrm{VAR}}}+\frac{1}{2}\sum_{k=1}^{L}\parallel\boldsymbol{y}_{k}-\mathcal{H}_{k}\circ\mathcal{M}_{k:1}(\boldsymbol{x}_{0})\parallel_{\mathbf{R}_{k}}^{2},$
(3.69)

where
$\mathbf{B}_{\mathrm{4D}\text{-}\mathrm{VAR}}\longleftarrow\mathbf{B}_{\mathrm{3D}\text{-}\mathrm{VAR}}$.
A similar indexing to equation to Eq. (3.65) gives the sequential form over
newly incoming observations.

The incremental linearisation that was performed in 3D-VAR in Eq. (3.48) forms
the basis for the classical technique of incremental 4D-VAR. Suppose again
that an explicit Cholesky factor for the background covariance is given
$\mathbf{B}_{\mathrm{4D}\text{-}\mathrm{VAR}}:=\boldsymbol{\Sigma}\boldsymbol{\Sigma}^{\top}$,
where the state is written as a perturbation of the $i$-th iterate as

$\displaystyle\boldsymbol{x}_{0}:=\overline{\boldsymbol{x}}^{i}_{0}+\boldsymbol{\Sigma}\boldsymbol{w}.$
(3.70)

Taking a Taylor expansion of the composition of the nonlinear observation
model and the nonlinear state model at the $i$-th iterate,

$\displaystyle\mathcal{H}_{k}\circ\mathcal{M}_{k:1}\left(\boldsymbol{x}_{0}\right)=\mathcal{H}\circ\mathcal{M}_{k:1}\left(\overline{\boldsymbol{x}}^{i}_{0}\right)+\mathbf{H}_{k}\mathbf{M}_{k:1}\boldsymbol{\Sigma}\boldsymbol{w}+\mathcal{O}\left(\parallel\boldsymbol{w}\parallel^{2}\right),$
(3.71)

the incremental cost function is rendered

$\displaystyle\mathcal{J}_{\mathrm{4D}\text{-}\mathrm{VARI}}(\boldsymbol{w})=\frac{1}{2}\parallel\boldsymbol{w}\parallel^{2}+\sum_{k=1}^{L}\parallel\boldsymbol{y}_{k}-\mathcal{H}_{k}\circ\mathcal{M}_{k:1}\left(\overline{\boldsymbol{x}}^{i}_{0}\right)-\mathbf{H}_{k}\mathbf{M}_{k:1}\boldsymbol{\Sigma}\boldsymbol{w}\parallel_{\mathbf{R}_{k}}^{2}.$
(3.72)

The argmin of Eq. (3.72) is defined as $\overline{\boldsymbol{w}}^{i}$, where
$\overline{\boldsymbol{x}}^{i+1}$ is defined as in Eq. (3.47).

It is important to remember that in Eq. (3.72), the terms $\mathbf{M}_{k:1}$
involve the calculation of the tangent linear model with reference to the
underlying nonlinear solution
$\mathcal{M}_{L:1}\left(\overline{\boldsymbol{x}}^{i}_{0}\right)$, simulated
freely over the entire DAW. A key aspect to the efficiency of the incremental
4D-VAR approach is in the adjoint model computation of the gradient to this
cost function, commonly known as _backpropagation_ in statistical and machine
learning [74, see Chapter 7]. The _adjoint model_ with respect to the proposal
is defined as

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\tilde{\boldsymbol{\delta}}=\left(\nabla_{\boldsymbol{x}}\boldsymbol{f}(\overline{\boldsymbol{x}}^{i})\right)^{\top}\tilde{\boldsymbol{\delta}}.$
(3.73)

The solution to the linear adjoint model above defines the adjoint resolvent
matrix $\mathbf{M}_{k}^{\top}$, i.e., the transpose of the tangent linear
resolvent. Notice that for the composition of the tangent linear model
forecasts with $k\leq l$, the adjoint is given as

$\displaystyle\mathbf{M}_{l:k}^{\top}:=\left(\mathbf{M}_{l}\mathbf{M}_{l-1}\cdots\mathbf{M}_{k+1}\mathbf{M}_{k}\right)^{\top}=\mathbf{M}_{k}^{\top}\mathbf{M}_{k+1}^{\top}\cdots\mathbf{M}_{l-1}^{\top}\mathbf{M}_{l}^{\top}.$
(3.74)

This means that the adjoint model variables $\tilde{\boldsymbol{\delta}}_{l}$
are propagated by the linear resolvents of the adjoint model, but applied in
reverse chronological order from the tangent linear model from last time
$t_{l}$ to the initial time $t_{k-1}$,

$\displaystyle\tilde{\boldsymbol{\delta}}_{k-1}:=\mathbf{M}^{\top}_{l:k}\tilde{\boldsymbol{\delta}}_{l}$
(3.75)

transmitting the sensitivity from a future time back to a perturbation of
$\overline{\boldsymbol{x}}_{k-1}^{i}$ [51, see Section 6.3].

Define the _innovation_ vector of the $i$-th iterate and the $k$-th
observation vector as

$\displaystyle\overline{\boldsymbol{\delta}}^{i}_{k}:=\boldsymbol{y}_{k}-\mathcal{H}_{k}\circ\mathcal{M}_{k:1}\left(\overline{\boldsymbol{x}}^{i}_{0}\right).$
(3.76)

Then the gradient of Eq. (3.72) with respect to the weight vector is written

$\displaystyle\nabla_{\boldsymbol{w}}\mathcal{J}_{\mathrm{4D}\text{-}\mathrm{VARI}}=\boldsymbol{w}-\sum_{k=1}^{L}\boldsymbol{\Sigma}^{\top}\mathbf{M}_{k:1}^{\top}\mathbf{H}_{k}^{\top}\mathbf{R}_{k}^{-1}\left[\overline{\boldsymbol{\delta}}^{i}_{k}-\mathbf{H}_{k}\mathbf{M}_{k:1}\boldsymbol{\Sigma}\boldsymbol{w}\right].$
(3.77)

Making the substitution

$\displaystyle\boldsymbol{\Delta}_{k}:=\mathbf{R}_{k}^{-1}\left[\overline{\boldsymbol{\delta}}^{i}_{k}-\mathbf{H}_{k}\mathbf{M}_{k:1}\boldsymbol{\Sigma}\boldsymbol{w}\right],$
(3.78)

notice that the gradient is written in Horner factorisation as

$\displaystyle\nabla_{\boldsymbol{w}}\mathcal{J}_{\mathrm{4D}\text{-}\mathrm{VARI}}=\boldsymbol{w}-\boldsymbol{\Sigma}^{\top}\mathbf{M}_{1}^{\top}\left[\mathbf{H}^{\top}_{1}\boldsymbol{\Delta}_{1}+\mathbf{M}^{\top}_{2}\left[\mathbf{H}^{\top}_{2}\boldsymbol{\Delta}_{2}+\cdots+\left[\mathbf{M}_{L}^{\top}\mathbf{H}^{\top}_{L}\boldsymbol{\Delta}_{L}\right]\right]\right].$
(3.79)

Defining the adjoint variable as
$\tilde{\boldsymbol{\delta}}_{L}:=\mathbf{H}_{L}^{\top}\boldsymbol{\Delta}_{L}$,
with recursion in $t_{k}$,

$\displaystyle\tilde{\boldsymbol{\delta}}_{k}$
$\displaystyle:=\mathbf{H}_{k}^{\top}\boldsymbol{\Delta}_{k}+\mathbf{M}_{k+1}^{\top}\tilde{\boldsymbol{\delta}}_{k+1},$
(3.80)

the gradient of the cost function is written as

$\displaystyle\nabla_{\boldsymbol{w}}\mathcal{J}_{\mathrm{4D}\text{-}\mathrm{VARI}}:=\boldsymbol{w}-\boldsymbol{\Sigma}\tilde{\boldsymbol{\delta}}_{0}.$
(3.81)

Using the above definitions, the gradient of the incremental cost function is
computed from the following steps: (i) a forwards pass of the free evolution
of $\boldsymbol{x}_{0}^{i}$ under the nonlinear forecast model, computing the
innovations as in Eq. (3.76); (ii) the propagation of the perturbation
$\boldsymbol{\Sigma}\boldsymbol{w}$ in the tangent linear model and the
linearised observation operator with respect to the proposal,
$\overline{\boldsymbol{x}}^{i}_{k}$, in order to compute the terms in Eq.
(3.78); and (iii) the back propagation of the sensitivities in the adjoint
model, Eq. (3.73), to obtain the adjoint variables recursively as in Eq.
(3.80) and the gradient from Eq. (3.81).

The benefit of the above approach is that this gives an extremely efficient
calculation of the gradient, provided that the tangent linear and adjoint
models are available, which, in turn, is key to very efficient numerical
optimisations of cost functions. However, the trade off is that the tangent
linear and adjoint models of a full scale geophysical state model, and of the
observation model, require considerable development time and expertise.
Increasingly, these models can be computed abstractly by differentiating a
computer program alone, in what is known as _automatic differentiation_ of
code [39, 40, 46, 6]. When tangent linear and adjoint models are not
available, one alternative is to use ensemble sampling techniques in the fully
nonlinear model $\mathcal{M}$ alone. The ensemble-based analysis provides a
complementary approach to the explicit tangent linear and adjoint model
analysis – this approach can be developed independently, or hybridised with
the use of the tangent linear and adjoint models as in various flavours of
hybrid _ensemble-variational_ (EnVAR) techniques [4, 5]. This tutorial focuses
on the independent development of EnVAR estimators in the following section.

##  3.4 Bayesian Ensemble-Variational Estimators

###  3.4.1 The Ensemble Transform Kalman Filter

Consider once again the perfect, linear-Gaussian model in Eq. (3.25). Rather
than explicitly computing the evolution of the background mean and error
covariance in the linear model as in Eq. (3.32) and the KF, one can
alternatively estimate the state mean and the error covariances using a
statistical sampling approach. Let
$\left\\{\boldsymbol{x}_{k,i}^{\mathrm{f/a}}\right\\}_{i=1}^{N_{e}}$ be
replicates of the model state, _independently and identically distributed_
(iid) according to the distribution

$\displaystyle\boldsymbol{x}_{k,i}^{\mathrm{f/a}}\sim\mathcal{N}\left(\overline{\boldsymbol{x}}^{\mathrm{f/a}}_{k},\mathbf{B}_{k}^{\mathrm{f/a}}\right).$
(3.82)

Given the iid assumption above, the ensemble-based mean and the ensemble-based
covariance

$\displaystyle\hat{\boldsymbol{x}}^{\mathrm{f/a}}_{k}$
$\displaystyle:=\frac{1}{N_{e}}\sum_{i=1}^{N_{e}}\boldsymbol{x}_{k,i}^{\mathrm{f/a}},$
(3.83a) $\displaystyle\mathbf{P}_{k}^{\mathrm{f/a}}$
$\displaystyle:=\frac{1}{N_{e}-1}\sum_{i=1}^{N_{e}}\left(\boldsymbol{x}_{k,i}^{\mathrm{f/a}}-\hat{\boldsymbol{x}}^{\mathrm{f/a}}_{k}\right)\left(\boldsymbol{x}_{k,i}^{\mathrm{f/a}}-\hat{\boldsymbol{x}}^{\mathrm{f/a}}_{k}\right)^{\top},$
(3.83b)

are (asymptotically) consistent estimators of the background, i.e.,

$\displaystyle\mathbb{E}\left[\hat{\boldsymbol{x}}^{\mathrm{f/a}}_{k}\right]=\overline{\boldsymbol{x}}^{\mathrm{f/a}}_{k},$
$\displaystyle\lim_{N_{e}\rightarrow\infty}\mathbb{E}\left[\mathbf{P}_{k}^{\mathrm{f/a}}\right]=\mathbf{B}_{k}^{\mathrm{f/a}},$
(3.84)

where the expectation is over the possible realisations of the random sample.
Particularly, the multivariate central limit theorem gives that

$\displaystyle\left(\mathbf{P}_{k}^{\mathrm{f/a}}\right)^{-\frac{1}{2}}\left(\boldsymbol{x}^{\mathrm{f/a}}_{k}-\hat{\boldsymbol{x}}^{\mathrm{f/a}}_{k}\right)\rightarrow_{d}\mathcal{N}\left(\boldsymbol{0},\mathbf{I}_{N_{x}}\right),$
(3.85)

referring to convergence in distribution as $N_{e}\rightarrow\infty$ [45, see
Section 6.2].

Using the relationships in Eq. (3.83), these estimators are efficiently
encoded as linear operations on the ensemble matrix. Define the _ensemble
matrix_ and the _perturbation matrix_ as

$\displaystyle\mathbf{E}^{\mathrm{f/a}}_{k}$
$\displaystyle:=\begin{pmatrix}\boldsymbol{x}_{k,1}^{\mathrm{f/a}}&\cdots&\boldsymbol{x}_{k,N_{e}}^{\mathrm{f/a}}\end{pmatrix}\in\mathbb{R}^{N_{x}\times
N_{e}},$ (3.86) $\displaystyle\mathbf{X}^{\mathrm{f/a}}_{k}$
$\displaystyle:=\begin{pmatrix}\boldsymbol{x}_{k,1}^{\mathrm{f/a}}-\hat{\boldsymbol{x}}_{k}^{\mathrm{f/a}}&\cdots&\boldsymbol{x}_{k,N_{e}}^{\mathrm{f/a}}-\hat{\boldsymbol{x}}_{k}^{\mathrm{f/a}}\end{pmatrix}\in\mathbb{R}^{N_{x}\times
N_{e}},$ (3.87)

i.e., as the arrays with the columns given as the ordered replicates of the
model state and their deviations from the ensemble mean respectively. For
$\boldsymbol{1}$ defined as the vector composed entirely of ones, define the
following linear operations conformally in their dimensions

$\displaystyle\hat{\boldsymbol{x}}_{k}^{\mathrm{f/a}}$
$\displaystyle=\mathbf{E}^{\mathrm{f/a}}_{k}\boldsymbol{1}/N_{e},$ (3.88a)
$\displaystyle\mathbf{X}_{k}^{\mathrm{f/a}}$
$\displaystyle=\mathbf{E}_{k}^{\mathrm{f/a}}\left(\mathbf{I}_{N_{e}}-\boldsymbol{1}\boldsymbol{1}^{\top}/N_{e}\right),$
(3.88b) $\displaystyle\mathbf{P}_{k}^{\mathrm{f/a}}$
$\displaystyle=\left(\mathbf{X}_{k}^{\mathrm{f/a}}\right)\left(\mathbf{X}_{k}^{\mathrm{f/a}}\right)^{\top}/(N_{e}-1).$
(3.88c)

The operator $\mathbf{I}_{N_{e}}-\boldsymbol{1}\boldsymbol{1}^{\top}/N_{e}$ is
the orthogonal complementary projection operator to the span of the vector of
ones, known as the _centring operator_ in statistics [45, see Section 6.1];
this has the effect of transforming the ensemble to mean zero.

Notice, when $N_{e}\leq N_{x}$ the ensemble-based error covariance has rank of
at most $N_{e}-1$ irrespective of the rank of the background error covariance,
corresponding to the one degree of freedom lost in computing the ensemble
mean. Therefore, to utilise the ensemble error covariance in the least-squares
optimisation as with the KF in Eq. (3.39), a new construction is necessary.
For a generic matrix $\mathbf{A}\in\mathbb{R}^{N\times M}$ with full column
rank $M$, define the (left Moore-Penrose) _pseudo-inverse_ [63, see page 423]

$\displaystyle\mathbf{A}^{\dagger}:=\left(\mathbf{A}^{\top}\mathbf{A}\right)^{-1}\mathbf{A}^{\top}.$
(3.89)

In particular, $\mathbf{A}^{\dagger}\mathbf{A}=\mathbf{I}_{M}$ and the
orthogonal projector into the column span of $\mathbf{A}$ is defined by
$\mathbf{A}\mathbf{A}^{\dagger}$. When $\mathbf{A}$ has full column rank as
above, define the Mahalanobis ‘distance’ with respect to
$\mathbf{G}=\mathbf{A}\mathbf{A}^{\top}$ as

$\displaystyle\parallel\boldsymbol{v}\parallel_{\mathbf{G}}:=\sqrt{\left(\mathbf{A}^{\dagger}\boldsymbol{v}\right)^{\top}\left(\mathbf{A}^{\dagger}\boldsymbol{v}\right)}.$
(3.90)

Note that in the case that $\mathbf{G}$ does not have full column rank, i.e.,
$N>M$, this is not a true norm on $\mathbb{R}^{N}$ as it is degenerate in the
null space of $\mathbf{A}^{\dagger}$. This instead represents a lift of a non-
degenerate norm in the column span of $\mathbf{A}$ to $\mathbf{R}^{N}$. In the
case that $\boldsymbol{v}$ is in the column span of $\mathbf{A}$, one can
equivalently write

$\displaystyle\boldsymbol{v}=\mathbf{A}\boldsymbol{w},$ (3.91a)
$\displaystyle\parallel\boldsymbol{v}\parallel_{\mathbf{G}}=\parallel\boldsymbol{w}\parallel,$
(3.91b)

for a vector of weights $\boldsymbol{w}\in\mathbf{R}^{M}$.

The _ensemble Kalman filter_ (EnKF) cost function for the linear-Gaussian
model is defined

$\displaystyle\mathcal{J}_{\mathrm{EnKF}}(\boldsymbol{x}_{k}):=\frac{1}{2}\parallel\hat{\boldsymbol{x}}_{k}^{\mathrm{f}}-\boldsymbol{x}_{k}\parallel_{\mathbf{P}_{k}^{\mathrm{f}}}^{2}+\frac{1}{2}\parallel\boldsymbol{y}_{k}-\mathbf{H}_{k}\boldsymbol{x}_{k}\parallel_{\mathbf{R}_{k}}^{2}$
(3.92a) $\displaystyle\Leftrightarrow\,\,$
$\displaystyle\mathcal{J}_{\mathrm{EnKF}}(\boldsymbol{w}):=\frac{1}{2}(N_{e}-1)\parallel\boldsymbol{w}\parallel^{2}+\frac{1}{2}\parallel\boldsymbol{y}_{k}-\mathbf{H}_{k}\hat{\boldsymbol{x}}^{\mathrm{f}}_{k}-\mathbf{H}_{k}\mathbf{X}_{k}^{\mathrm{f}}\boldsymbol{w}\parallel_{\mathbf{R}_{k}}^{2},$
(3.92b)

where the model state is written as a perturbation of the ensemble mean

$\displaystyle\boldsymbol{x}_{k}=\hat{\boldsymbol{x}}_{k}^{\mathrm{f}}+\mathbf{X}_{k}^{\mathrm{f}}\boldsymbol{w}.$
(3.93)

Notice, $\boldsymbol{w}\in\mathbb{R}^{N_{e}}$, giving the linear combination
of the ensemble perturbations used to represent the model state.

Define $\hat{\boldsymbol{w}}$ to be the argmin of the cost function in Eq.
(3.92b). [48] and [8] demonstrate that, up to a gauge transformation,
$\hat{\boldsymbol{w}}$ yields the argmin of the state-space cost function, Eq.
(3.92a), when the estimate is restricted to the ensemble span. Equation
(3.92b) is quadratic in $\boldsymbol{w}$ and can be solved to render

$\displaystyle\hat{\boldsymbol{w}}$
$\displaystyle:=\boldsymbol{0}-\boldsymbol{\Xi}^{-1}_{\mathcal{J}_{\mathrm{EnKF}}}\nabla_{\boldsymbol{w}}\mathcal{J}_{\mathrm{EnKF}}|_{\boldsymbol{w}=\boldsymbol{0}},$
(3.94a) $\displaystyle\mathbf{T}$
$\displaystyle:=\boldsymbol{\Xi}^{-\frac{1}{2}}_{\mathcal{J}_{\mathrm{EnKF}}},$
(3.94b) $\displaystyle\mathbf{P}_{k}^{\mathrm{a}}$
$\displaystyle=\left(\mathbf{X}_{k}^{\mathrm{f}}\mathbf{T}\right)\left(\mathbf{X}_{k}^{\mathrm{f}}\mathbf{T}\right)^{\top}/(N_{e}-1),$
(3.94c)

corresponding to a single iteration of Newton’s descent algorithm in Eq.
(3.94a), initialised with the ensemble mean, to find the optimal weights.

The linear _ensemble transform Kalman filter_ (ETKF) equations are then given
by [7, 48]

$\displaystyle\mathbf{E}_{k}^{\mathrm{f}}$
$\displaystyle=\mathbf{M}_{k}\mathbf{E}_{k-1}^{\mathrm{a}},$ (3.95a)
$\displaystyle\mathbf{E}^{\mathrm{a}}_{k}$
$\displaystyle=\hat{\boldsymbol{x}}_{k}^{\mathrm{f}}\boldsymbol{1}^{\top}+\mathbf{X}_{k}^{\mathrm{f}}\left(\hat{\boldsymbol{w}}\boldsymbol{1}^{\top}+\sqrt{N_{e}-1}\mathbf{T}\mathbf{U}\right),$
(3.95b)

where $\mathbf{U}\in\mathbb{R}^{N_{e}\times N_{e}}$ can be any mean-
preserving, orthogonal transformation, i.e.,
$\mathbf{U}\boldsymbol{1}=\boldsymbol{1}$. The simple choice of
$\mathbf{U}:=\mathbf{I}_{N_{e}}$ is sufficient, but it has been demonstrated
that choosing a random, mean-preserving orthogonal transformation at each
analysis as above can improve the accuracy and robustness of the ETKF,
smoothing out higher-order artefacts in the empirical covariance estimate
[79].

Notice that Eq. (3.95b) is written equivalently as a single right ensemble
transformation:

$\displaystyle\mathbf{E}_{k}^{\mathrm{a}}$
$\displaystyle=\mathbf{E}_{k}^{\mathrm{f}}\boldsymbol{\Psi}_{k},$ (3.96a)
$\displaystyle\boldsymbol{\Psi}_{k}$
$\displaystyle:=\boldsymbol{1}\boldsymbol{1}^{\top}/N_{e}+\left(\mathbf{I}_{N_{e}}-\boldsymbol{1}\boldsymbol{1}^{\top}/N_{e}\right)\left(\hat{\boldsymbol{w}}\boldsymbol{1}^{\top}+\sqrt{N_{e}-1}\mathbf{T}\mathbf{U}\right),$
(3.96b)

where the columns are approximately distributed as in Eq. (3.82), with the
(asymptotic) consistency as with the central limit theorem, Eq. (3.85).
However, in the small feasible sample sizes for realistic geophysical models,
this approximation can lead to a systematic underestimation of the uncertainty
of the analysis state, where the ensemble-based error covariance can become
overly confident in its own estimate with artificially small variances.
Covariance inflation is a technique that is widely used to regularise the
ensemble-based error covariance by increasing the empirical variances of the
estimate. This can be used to handle the inaccuracy of the estimator due to
the finite sample size approximation of the background mean and error
covariance as in Eq. (3.85), as well as inaccuracies due to a variety of other
sources of error [18, 71, 84].

A Bayesian hierarchical approach can model the inaccuracy in the approximation
error due to the finite sample size by including a prior additionally on the
background mean and error covariance
$p\left(\overline{\boldsymbol{x}}_{k}^{\mathrm{f}},\mathbf{B}_{k}^{\mathrm{f}}\right)$,
as in the finite-size ensemble Kalman filter formalism of [8], [12] and [11].
Mathematical results demonstrate that covariance inflation can ameliorate the
systematic underestimation of the variances due to model error in the presence
of a low rank ensemble [42]. In the presence of significant model error, the
finite-size analysis is extended by the variant developed by [71].

The linear transform of the ensemble matrix in Eq. (3.96) is key to the
efficiency of the (deterministic) EnKF presented above. The ensemble-based
cost function Hessian
$\boldsymbol{\Xi}_{\mathcal{J}_{\mathrm{EnKF}}}\in\mathbb{R}^{N_{e}\times
N_{e}}$, where $N_{e}\ll N_{x}$ for typical geophysical models. Particularly,
the cost of computing the optimal weights as in Eq. (3.94a) and computing the
transform matrix $\mathbf{T}$ in Eq. (3.94b) are both subordinate to the cost
of the eigenvalue decomposition of the Hessian at
$\mathcal{O}\left(N_{e}^{3}\right)$ floating point operations (flops), or to a
randomised singular value decomposition [34]. However, the extremely low
ensemble size means that the correction to the forecast is restricted to the
low-dimensional ensemble span, which may fail to correct directions of rapidly
growing errors due to the rank deficiency. While chaotic, dissipative dynamics
implies that the background covariance $\mathbf{B}_{k}^{\mathrm{f/a}}$ has
spectrum concentrated on a reduced rank subspace [19, and referenced therein],
covariance hybridisation [68] or localisation [76] are used in practice to
regularise the estimator’s extreme rank deficiency, and the spurious empirical
correlations that occur as a result of the degenerate sample size.

When extended to nonlinear dynamics, as in Eq. (3.7a), the EnKF can be seen to
make an ensemble-based approximation to the EKF cost function, where the
forecast ensemble is defined by

$\displaystyle\boldsymbol{x}_{k,i}^{\mathrm{f}}:=\mathcal{M}_{k-1}\left(\boldsymbol{x}_{k-1,i}^{\mathrm{a}}\right),$
$\displaystyle\mathbf{E}_{k}^{\mathrm{f}}:=\begin{pmatrix}\boldsymbol{x}_{k,1}^{\mathrm{f}}&\cdots&\boldsymbol{x}_{k,N_{e}}^{\mathrm{f}}\end{pmatrix}.$
(3.97)

The accuracy of the linear-Gaussian approximation to the dynamics of the
nonlinear evolution of the ensemble, like the approximation of the EKF,
depends strongly on the length of the forecast horizon $\Delta_{t}$. When the
ensemble mean is a sufficiently accurate approximation of the mean state, and
if the ensemble spread is of the same order as the error in the mean estimate
[90], a similar approximation can be made for the ensemble evolution at first
order as with Eq. (3.53) but linearised about the ensemble mean as discussed
later in Eq. (3.112). Despite the similarity to the EKF, in the moderately
nonlinear dynamics present in medium- to longer-range forecast horizons, the
EnKF does not suffer from the same inaccuracy as the EKF in truncating the
time-evolution at first order [31]. The forecast ensemble members themselves
evolve fully nonlinearly, tracking the higher-order dynamics, but the analysis
update based on the Gaussian approximation becomes increasingly biased and
fails to discriminate features like multi-modality of the posterior, even
though this is occasionally an asset with higher-order statistical artefacts
[[, see the discussion in]]lawson2004.

###  3.4.2 The Maximum Likelihood Ensemble Filter

The EnKF filter analysis automatically accommodates weak nonlinearity in the
state model without using the tangent linear model, as the filtering cost
function has no underlying dependence on the state model. However, the EnKF
analysis must be adjusted to account for nonlinearity in the observation
model. When a nonlinear observation operator is introduced, as in Eq. (3.7b),
the EnKF cost function can be re-written in the incremental analysis as

$\displaystyle\mathcal{J}_{\mathrm{EnKFI}}(\boldsymbol{w}):=\frac{1}{2}(N_{e}-1)\parallel\boldsymbol{w}\parallel^{2}+\frac{1}{2}\parallel\boldsymbol{y}_{k}-\mathcal{H}_{k}\left(\hat{\boldsymbol{x}}_{k}^{\mathrm{i,f}}\right)-\mathbf{H}_{k}\mathbf{X}_{k}^{\mathrm{f}}\boldsymbol{w}\parallel_{\mathbf{R}_{k}}^{2},$
(3.98)

where $\hat{\boldsymbol{x}}_{k}^{i,\mathrm{f}}$ refers to the $i$-th iteration
for the forecast mean. Note that this does not refer to an estimate derived by
an iterated ensemble forecast through the nonlinear state model – rather, this
is an iteration only with respect to the estimate of the optimal weights for
the forecast perturbations. If $\hat{\boldsymbol{w}}^{i}$ is defined as the
argmin of the cost function in Eq. (3.98), then the iterations of the ensemble
mean are given as

$\displaystyle\hat{\boldsymbol{x}}_{k}^{i+1,\mathrm{f}}:=\hat{\boldsymbol{x}}^{i,\mathrm{f}}_{k}+\mathbf{X}_{k}^{\mathrm{f}}\hat{\boldsymbol{w}}^{i}.$
(3.99)

When $\parallel\hat{\boldsymbol{w}}^{i}\parallel$ is sufficiently small, the
optimisation terminates and the transform and the ensemble update can be
performed as with the ETKF as in Eqs. (3.94b) and (3.96). However, the direct,
incremental approach above used the computation of the Jacobian of the
observation operator $\mathbf{H}_{k}$.

The _maximum likelihood ensemble filter_ (MLEF) of [92] and [93] is an
estimator designed to perform the above incremental analysis in the ETKF
formalism, but without taking an explicit Taylor expansion for the observation
operator above. The method is designed to approximate the directional
derivative of the nonlinear observation operator with respect to the ensemble
perturbations,

$\displaystyle\mathbf{H}_{k}\mathbf{X}_{k}^{\mathrm{f}}:=\nabla|_{\hat{\boldsymbol{x}}^{i,\mathrm{f}}_{k}}\left[\mathcal{H}_{k}\right]\mathbf{X}^{\mathrm{f}}_{k}$
(3.100)

equivalent to computing the ensemble sensitivities of the map linearised about
the ensemble-based mean. This is often performed with an explicit finite-
differences approximation between the ensemble members and the ensemble mean,
mapped to the observation space. Particularly, one may write

$\displaystyle\mathcal{H}_{k}\left(\boldsymbol{x}_{k}^{i,\mathrm{f}}\right)\approx\hat{\boldsymbol{y}}^{i}_{k}$
$\displaystyle:=\mathcal{H}_{k}\left(\hat{\boldsymbol{x}}_{k}^{i,\mathrm{f}}\boldsymbol{1}^{\top}+\epsilon\mathbf{X}_{k}^{\mathrm{f}}\right)\boldsymbol{1}/N_{e}$
(3.101a)
$\displaystyle\mathbf{H}_{k}\mathbf{X}_{k}^{\mathrm{f}}\approx\widetilde{\mathbf{Y}}_{k}$
$\displaystyle:=\frac{1}{\epsilon}\mathcal{H}_{k}\left(\hat{\boldsymbol{x}}_{k}^{i,\mathrm{f}}\boldsymbol{1}^{\top}+\epsilon\mathbf{X}_{k}^{\mathrm{f}}\right)\left(\mathbf{I}_{N_{e}}-\boldsymbol{1}\boldsymbol{1}^{\top}/N_{e}\right),$
(3.101b)

where $\epsilon$ is a small constant that re-scales the ensemble perturbations
to approximate infinitesimals about the mean, and rescales the finite
differences about the ensemble mean in the observation space. This technique
is used, e.g., in a modern form of the MLEF algorithm, based on the analysis
of the iterative ensemble Kalman filter and smoother [4, see Section 6.7.2.1].
The approximation above easily generalises to a 4D analysis, taking the
directional derivative with respect to the state and observation models
simultaneously, as in Eq. (3.111). However, in the smoothing problem, one can
also extend the above analysis in terms of an alternating forwards filtering
pass and backwards filtering pass to estimate the joint posterior. This
tutorial returns to the 4D analysis with the iterative ensemble Kalman filter
and smoother in Section 3.4.4 after an interlude on the retrospective analysis
of the EnKS in the following section.

###  3.4.3 The Ensemble Transform Kalman Smoother

The EnKS extends the filter analysis in the ETKF over the smoothing DAW by
sequentially re-analysing past states with future observations with an
additional filtering pass over the DAW backwards-in-time. This analysis is
performed retrospectively in the sense that the filter cycle of the ETKF is
left unchanged, while an additional inner-loop of the DA cycle performs an
update on the estimated lagged state ensembles within the DAW, stored in
memory. This can be formulated both for a fixed DAW and for fixed-lag
smoothing, where only minor modifications are needed. Consider here the
algorithmically stationary DAW $\\{t_{1},\cdots,t_{L}\\}$ of fixed-lag
smoothing, with a shift $S$ and lag $L$, and where it is assumed that $S=1\leq
L$. The fixed-lag smoothing cycle of the EnKS begins by estimating the joint
posterior pdf $p\left(\boldsymbol{x}_{L:1}|\boldsymbol{y}_{L:1}\right)$
recursively, given the joint posterior estimate over the last DAW
$p\left(\boldsymbol{x}_{L-1:0}|\boldsymbol{y}_{L-1:0}\right)$.

Given $p(\boldsymbol{x}_{L-1:0},\boldsymbol{y}_{L-1:0})$, one can write the
filtering pdf up to proportionality:

$\displaystyle p(\boldsymbol{x}_{L}|\boldsymbol{y}_{L:0})$
$\displaystyle\propto
p(\boldsymbol{y}_{L}|\boldsymbol{x}_{L},\boldsymbol{y}_{L-1:0})p(\boldsymbol{x}_{L},\boldsymbol{y}_{L-1:0})$
(3.102a)
$\displaystyle\propto\underbrace{p(\boldsymbol{y}_{L}|\boldsymbol{x}_{L})}_{(i)}\underbrace{\int
p(\boldsymbol{x}_{L}|\boldsymbol{x}_{L-1})p(\boldsymbol{x}_{L-1:0}|\boldsymbol{y}_{L-1:0})\mathrm{d}\boldsymbol{x}_{L-1:0}}_{(ii)},$
(3.102b)

as the product of (i) the likelihood of the observation $\boldsymbol{y}_{L}$
given $\boldsymbol{x}_{L}$; and (ii) the forecast for $\boldsymbol{x}_{L}$
using the transition pdf on the last joint posterior estimate, marginalising
out the past history of the model state $\boldsymbol{x}_{L-1:0}$. Recalling
that $p(\boldsymbol{x}_{L}|\boldsymbol{y}_{L:1})\propto
p(\boldsymbol{x}_{L}|\boldsymbol{y}_{L:0})$, this provides a means to estimate
the filter marginal of the joint posterior. An alternating filtering pass,
backwards-in-time, completes the smoothing cycle by estimating the joint
posterior pdf $p(\boldsymbol{x}_{L:1},\boldsymbol{y}_{L:1})$.

Consider that the marginal smoother pdf is proportional to

$\displaystyle p(\boldsymbol{x}_{L-1}|\boldsymbol{y}_{L:0})$
$\displaystyle\propto
p(\boldsymbol{y}_{L}|\boldsymbol{x}_{L-1},\boldsymbol{y}_{L-1:0})p(\boldsymbol{x}_{L-1},\boldsymbol{y}_{L-1:0})$
(3.103a)
$\displaystyle\propto\underbrace{p(\boldsymbol{y}_{L}|\boldsymbol{x}_{L-1})}_{(i)}\underbrace{p(\boldsymbol{x}_{L-1}|\boldsymbol{y}_{L-1:0})}_{(ii)},$
(3.103b)

where: (i) is the likelihood of the observation $\boldsymbol{y}_{L}$ given the
past state $\boldsymbol{x}_{L-1}$; (ii) is the marginal pdf for
$\boldsymbol{x}_{L-1}$ from the last joint posterior. The corresponding
linear-Gaussian Bayesian MAP cost function is given for the retrospective
analysis of the KS as

$\displaystyle\mathcal{J}_{\mathrm{KS}}(\boldsymbol{x}_{L-1})=\frac{1}{2}\parallel\boldsymbol{x}_{L-1}-\overline{\boldsymbol{x}}_{L-1|L-1}^{\mathrm{s}}\parallel^{2}_{\mathbf{B}_{L-1|L-1}^{\mathrm{s}}}+\frac{1}{2}\parallel\boldsymbol{y}_{L}-\mathbf{H}_{L}\mathbf{M}_{L}\boldsymbol{x}_{L-1}\parallel_{\mathbf{R}_{L}}^{2},$
(3.104)

where $\overline{\boldsymbol{x}}_{L-1|L-1}^{\mathrm{s}}$ and
$\mathbf{B}^{\mathrm{s}}_{L-1|L-1}$ are the mean and covariance of the
marginal smoother pdf $p(\boldsymbol{x}_{L-1}|\boldsymbol{y}_{L-1:0})$. Define
the matrix decomposition with the factorisation (e.g., a Cholesky
decomposition):

$\displaystyle\mathbf{B}_{L-1|L-1}^{\mathrm{s}}=\boldsymbol{\Sigma}_{L-1|L-1}^{\mathrm{s}}\left(\boldsymbol{\Sigma}_{L-1|L-1}^{\mathrm{s}}\right)^{\top},$
(3.105)

and write
$\boldsymbol{x}_{L-1}=\overline{\boldsymbol{x}}_{L-1|L-1}^{\mathrm{s}}+\boldsymbol{\Sigma}_{L-1|L-1}^{\mathrm{s}}\boldsymbol{w}$,
rendering the cost function as

$\displaystyle\mathcal{J}_{\mathrm{KS}}(\boldsymbol{w})$
$\displaystyle=\frac{1}{2}\parallel\boldsymbol{w}\parallel^{2}+\frac{1}{2}\parallel\boldsymbol{y}_{L}-\mathbf{H}_{L}\mathbf{M}_{L}(\overline{\boldsymbol{x}}_{L-1|L-1}^{\mathrm{s}}+\boldsymbol{\Sigma}_{L-1|L-1}^{\mathrm{s}}\boldsymbol{w})\parallel_{\mathbf{R}_{L}}^{2}$
(3.106a)
$\displaystyle=\frac{1}{2}\parallel\boldsymbol{w}\parallel^{2}+\frac{1}{2}\parallel\boldsymbol{y}_{L}-\mathbf{H}_{L}\overline{\boldsymbol{x}}_{L}^{\mathrm{f}}-\mathbf{H}_{L}\boldsymbol{\Sigma}_{L}^{\mathrm{f}}\boldsymbol{w}\parallel_{\mathbf{R}_{L}}^{2}.$
(3.106b)

Let $\overline{\boldsymbol{w}}$ now denote the argmin of Eq. (3.106). It is
important to recognise that

$\displaystyle\boldsymbol{x}_{L}:=$
$\displaystyle\mathbf{M}_{L}\left(\overline{\boldsymbol{x}}_{L-1|L-1}^{\mathrm{s}}+\boldsymbol{\Sigma}_{L-1|L-1}^{\mathrm{s}}\boldsymbol{w}\right)$
(3.107a) $\displaystyle=$
$\displaystyle\overline{\boldsymbol{x}}_{L}^{\mathrm{f}}+\boldsymbol{\Sigma}_{L}^{\mathrm{f}}\boldsymbol{w},$
(3.107b)

such that the argmin for the smoothing problem $\overline{\boldsymbol{w}}$ is
also the argmin for the filtering MAP analysis.

The ensemble-based approximation,

$\displaystyle\boldsymbol{x}_{L-1}$
$\displaystyle=\hat{\boldsymbol{x}}_{L-1|L-1}^{\mathrm{s}}+\mathbf{X}_{L-1|L-1}^{\mathrm{s}}\boldsymbol{w},$
(3.108a) $\displaystyle\mathcal{J}_{\mathrm{EnKS}}(\boldsymbol{w})$
$\displaystyle=\frac{1}{2}(N_{e}-1)\parallel\boldsymbol{w}\parallel^{2}+\frac{1}{2}\parallel\boldsymbol{y}_{L}-\mathbf{H}_{L}\hat{\boldsymbol{x}}^{\mathrm{f}}_{L}-\mathbf{H}_{L}\mathbf{X}_{L}^{\mathrm{f}}\boldsymbol{w}\parallel^{2},$
(3.108b)

to the exact smoother cost function in Eq. (3.106) yields the retrospective
analysis of the EnKS as

$\displaystyle\hat{\boldsymbol{w}}$
$\displaystyle:=\boldsymbol{0}-\boldsymbol{\Xi}_{\mathcal{J}_{\mathrm{EnKS}}}^{-1}\nabla\mathcal{J}_{\mathrm{EnKS}}|_{\boldsymbol{w}=\boldsymbol{0}},$
(3.109a) $\displaystyle\mathbf{T}$
$\displaystyle:=\boldsymbol{\Xi}_{\mathcal{J}_{\mathrm{EnKS}}}^{-\frac{1}{2}},$
(3.109b)
$\displaystyle\begin{split}\mathbf{E}_{L-1|L}^{\mathrm{s}}&=\hat{\boldsymbol{x}}_{L-1|L-1}^{\mathrm{s}}\boldsymbol{1}^{\top}+\mathbf{X}_{L-1|L-1}^{\mathrm{s}}\left(\hat{\boldsymbol{w}}\boldsymbol{1}^{\top}+\sqrt{N_{e}-1}\mathbf{T}\mathbf{U}\right),\\\
&=\mathbf{E}_{L-1|L-1}^{\mathrm{s}}\boldsymbol{\Psi}_{L}.\end{split}$ (3.109c)

where $\boldsymbol{\Psi}_{L}$ is the ensemble transform as defined for the
filtering update as in Eq. (3.96).

The above equations generalise for arbitrary indices $k|L$ over the DAW,
providing the complete description of the inner-loop between each filter cycle
of the EnKS. After each new observation is assimilated with the ETKF analysis
step, a smoother inner-loop makes a backwards pass over the DAW applying the
transform and the weights of the ETKF filter update to each past ensemble
state stored in memory. This analysis easily generalises to the case where
there is a shift of the DAW with $S>1$, though the EnKS alternating forwards
and backwards filtering passes must be performed in sequence over the
observations, ordered-in-time, rather than making a global analysis over
$\boldsymbol{y}_{L:L-S+1}$. Finally, this easily extends to accommodate a
nonlinear observation model by using the MLEF filtering step to obtain the
optimal ensemble transform, and by applying this recursively backwards-in-time
to the lagged ensemble states.

Figure 3.3: $L=4$ (lag), $S=1$ (shift) EnKS. Observations are assimilated
sequentially via the filter cost function and a retrospective re-analysis is
applied to all ensemble states within the lag window stored in memory. Source:
[41], adapted from [4].

A schematic of the EnKS cycle for a lag of $L=4$ and a shift of $S=1$ is
pictured in Fig. 3.3. Time moves forwards from left to right in the horizontal
axis with a step size of $\Delta_{t}$. At each analysis time, the ensemble
forecast from the last filter pdf is combined with the observation to produce
the ensemble transform update. This transform is then utilised to produce the
posterior estimate for all lagged ensemble states, conditioned on the new
observation. The information in the posterior estimate thus flows in reverse
time to the lagged states stored in memory, but the information flow is
unidirectional in this scheme. This type of retrospective analysis maintains
the computational cost of the fixed-lag EnKS described above at a comparable
level to the EnKF, with the only significant additional cost being the storage
of the ensemble at lagged times to be re-analysed.

###  3.4.4 The Iterative Ensemble Kalman Filter and Smoother

While the EnKS is computationally efficient, its retrospective analysis has
some drawbacks compared to the 4D analysis in terms of forecast accuracy. The
4D fixed-lag smoothing analysis re-initialises each cycle with a re-analysed
estimate for the initial data, transmitting the observations’ information
forwards-in-time through the nonlinear dynamics. The challenge with the 4D
analysis in the absence of the tangent linear and adjoint models is in
devising how to efficiently and accurately compute the gradient of the 4D cost
function. Building on [92, 58], the analysis of the _iterative ensemble Kalman
filter_ (IEnKF) and the _iterative ensemble Kalman smoother_ (IEnKS) extends
the ensemble transform method of the ETKF to iteratively optimise the 4D cost
function.

Recall the quadratic cost function at the basis of incremental 4D-VAR, Eq.
(3.72) – making an ensemble-based approximation for the background mean and
error covariance, and writing the model state as an ensemble perturbation of
the $i$-th proposal for the ensemble mean, one has

$\displaystyle\begin{split}\mathcal{J}_{\mathrm{IEnKS}}(\boldsymbol{w}):=&\frac{1}{2}\parallel\hat{\boldsymbol{x}}^{i}_{0}-\hat{\boldsymbol{x}}^{i}_{0}-\mathbf{X}_{0}\boldsymbol{w}\parallel_{\mathbf{P}_{0}}^{2}\\\
&+\frac{1}{2}\sum_{k=1}^{L}\parallel\boldsymbol{y}_{k}-\mathcal{H}_{k}\circ\mathcal{M}_{k:1}\left(\hat{\boldsymbol{x}}^{i}_{0}\right)-\mathbf{H}_{k}\mathbf{M}_{k:1}\mathbf{X}_{0}\boldsymbol{w}\parallel_{\mathbf{R}_{k}}^{2}\end{split}$
(3.110a)
$\displaystyle\begin{split}=&\frac{1}{2}\left(N_{e}-1\right)\parallel\boldsymbol{w}\parallel^{2}\\\
&+\frac{1}{2}\sum_{k=1}^{L}\parallel\boldsymbol{y}_{k}-\mathcal{H}_{k}\circ\mathcal{M}_{k:1}\left(\hat{\boldsymbol{x}}^{i}_{0}\right)-\mathbf{H}_{k}\mathbf{M}_{k:1}\mathbf{X}_{0}\boldsymbol{w}\parallel_{\mathbf{R}_{k}}^{2}.\end{split}$
(3.110b)

Applying the finite differences approximation as given in Eq. (3.101), but
with respect to the composition of the nonlinear observation operator with the
nonlinear state model

$\displaystyle\mathcal{H}_{k}\circ\mathcal{M}_{k:1}\left(\hat{\boldsymbol{x}}_{0}^{i}\right)\approx\hat{\boldsymbol{y}}_{k}$
$\displaystyle:=\mathcal{H}_{k}\circ\mathcal{M}_{k:1}\left(\hat{\boldsymbol{x}}^{i}_{0}+\mathbf{X}_{0}\epsilon\right)\boldsymbol{1}/N_{e},$
(3.111a)
$\displaystyle\mathbf{H}_{k}\mathbf{M}_{k:1}\mathbf{X}_{0}\approx\widetilde{\mathbf{Y}}_{k}$
$\displaystyle:=\frac{1}{\epsilon}\mathcal{H}_{k}\circ\mathcal{M}_{k:1}\left(\hat{\boldsymbol{x}}^{i}_{0}+\mathbf{X}_{0}\epsilon\right)\left(\mathbf{I}_{N_{e}}-\boldsymbol{1}\boldsymbol{1}^{\top}/N_{e}\right),$
(3.111b)

this sketches the ‘bundle’ formulation of the IEnKS [13, 14]. The above
indexing refers to the case where the smoothing problem is performed offline,
with a fixed DAW, though a similar indexing to Eq. (3.65) gives the sequential
form over newly incoming observations.

The sequential form of the IEnKS can be treated as a nonlinear sequential
filter, which is the purpose that the method was originally devised for.
Indeed, setting the number of lagged states $L=1$, this provides a direct
extension of the EnKF cost function, Eq. (3.92b), but where there is an
additional dependence on the initial conditions for the ensemble forecast.
This lag-1 iterative filtering scheme is called the IEnKF [80, 12], which
formed the original basis for the IEnKS. Modern forms of the IEnKF/S analysis
furthermore include the treatment of model errors, like weak-constraint 4D-VAR
[77, 78, 36]. Alternative formulations of this analysis, based on the original
stochastic, perturbed observation EnKF [30, 17], also have a parallel
development in the _ensemble randomised maximum likelihood method_ (EnRML) of
[43, 21, 72]. A similar ensemble-variational estimator based on the EnKF
analysis is the _ensemble Kalman inversion_ (EKI) of [49, 81, 54].

###  3.4.5 The Single-Iteration Ensemble Kalman Smoother

The sequential smoothing analysis so far has presented two classical
approaches: (i) a 4D approach using the global analysis of all new
observations available within a DAW at once, optimising an initial condition
for the lagged model state; and (ii) the 3D approach based upon alternating
forwards and backwards filtering passes over the DAW. Each of these approaches
has strengths and weaknesses in a computational cost / forecast accuracy trade
off. Particularly, the 4D approach, as in the IEnKS, benefits from the
improved estimate of the initial condition when producing the subsequent
forecast statistics in a shifted DAW; however, each step of the iterative
optimisation comes at the cost of simulating the ensemble forecast over the
entire DAW in the fully nonlinear state model, which is typically the greatest
numerical expense in geophysical DA. On the other hand, the 3D approach of the
classical EnKS benefits from a low computational cost, requiring only a single
ensemble simulation of the nonlinear forecast model over the DAW, while the
retrospective analysis of lagged states is performed with the filtering
transform without requiring additional model simulations (though at a
potentially large memory storage cost). It should be noted that in the
perfect, linear-Gaussian model Bayesian analysis, both approaches produce
equivalent estimates of the joint posterior. However, when nonlinearity is
present in the DA cycle the approaches produce distinct estimates – for this
reason, the source of the nonlinearity in the DA cycle is of an important
practical concern.

Consider the situation in which the forecast horizon $\Delta_{t}$ is short, so
that the ensemble time-evolution is weakly nonlinear, and where the ensemble
mean is a good approximation for the mean of the filtered distribution. In
this case, the model forecast dynamics are well-approximated by the tangent
linear evolution of a perturbation about the ensemble mean, i.e.,

$\displaystyle\boldsymbol{x}_{k}=\mathcal{M}_{k}\left(\boldsymbol{x}_{k-1}\right)\approx\mathcal{M}_{k}\left(\hat{\boldsymbol{x}}_{k-1}\right)+\mathbf{M}_{k}\boldsymbol{\delta}_{k-1}.$
(3.112)

For such short-range forecasts, nonlinearity in the observation-analysis-
forecast cycle may instead be dominated by the nonlinearity in the observation
operator $\mathcal{H}_{k}$, or in the optimisation of hyper-parameters of the
filtering cost function, and not by the model forecast itself. In this
situation, an iterative simulation of the forecast dynamics as in the 4D
approach may not produce a cost-effective reduction in the analysis error as
compared to, e.g., the MLEF filter analysis optimising the filtering cost
function alone. However, one may still obtain the benefits of re-
initialisation of the forecast with a re-analysed prior by using a simple
hybridisation of the 3D and 4D smoothing analyses.

Specifically, in a given fixed-lag smoothing cycle, one may iteratively
optimise the sequential filtering cost functions for a given DAW corresponding
to the new observations at times $\\{t_{L-S+1},\cdots,t_{L}\\}$ as with the
MLEF. These filtering ensemble transforms not only condition the ensemble at
the corresponding observation time instance, but also produce a retrospective
analysis of a lagged state as in Eq. (3.109c) with the EnKS. However, when the
DAW itself is shifted, one does not need to produce a forecast from the latest
filtered ensemble – instead one can initialise the next ensemble forecast with
the lagged, retrospectively smoothed ensemble at the beginning of the last
DAW. This hybrid approach utilising the retrospective analysis, as in the
classical EnKS, and the ensemble simulation over the lagged states while
shifting the DAW, as in the 4D analysis of the IEnKS, was recently developed
by [41], and is called the _single-iteration ensemble Kalman smoother_
(SIEnKS). The SIEnKS is named as such because it produces its forecast, filter
and re-analysed smoother statistics with a single iteration of the ensemble
simulation over the DAW in a fully consistent Bayesian analysis. By doing so,
it seeks to minimise the leading order cost of EnVAR smoothing, i.e., the
ensemble simulation in the nonlinear forecast model. However, the estimator is
free to iteratively optimise the filter cost function for any single
observation vector without additional iterations of the ensemble simulation.
In observation-analysis-forecast cycles in which the forecast error dynamics
are weakly nonlinear, yet other aspects of the cycle are moderately to
strongly nonlinear, this scheme is shown to produce a forecast accuracy
comparable to, and at times better than, the 4D approach but with an overall
lower leading-order computational burden [41].

Figure 3.4: $L=4$ (lag), $S=2$ (shift) SIEnKS diagram. An initial condition
from the last smoothing cycle initialises a forecast simulation over the
current DAW of $L=4$ states. New observations entering the DAW are assimilated
sequentially via the filter cost function. After each filter analysis, a
retrospective re-analysis is applied to the initial ensemble, and this re-
analysed initial condition is evolved via the model $S$ analysis times
forwards to begin the next cycle. Source: [41].

A schematic of the SIEnKS cycle for a lag of $L=4$ and a shift of $S=2$ is
pictured in Fig. 3.4. This demonstrates how the sequential analysis of the
filter cost function and sequential, retrospective re-analysis for each
incoming observation differs from the global analysis of the 4D approach of
the IEnKS. Other well-known DA schemes combining a retrospective re-analysis
and re-initialisation of the ensemble forecast include the _Running In Place_
(RIP) smoother of [52, 91] and the _One Step Ahead_ (OSA) smoother of [28] and
[2]. It can be shown that, with an ETKF style filter analysis, a single
iteration of the ensemble over the DAW, a perfect model assumption and a lag
of $L=S=1$, the SIEnKS, RIP and OSA smoothers all coincide [41].

##  3.5 Machine Learning and Data Assimilation in a Bayesian Perspective

This final section demonstrates how the Bayesian DA framework can be extended
to estimate more than the state vector, including key parameters of the model
and the observation-analysis-forecast cycle. In particular, the Bayesian
framework can be used to formulate techniques for learning both the state
vector and part of, if not the full, dynamical model. If this objective was
always in the scope of classical DA, it was actually made possible,
significantly beyond linear regression, by the introduction of _Machine
Learning_ (ML) techniques to supplement traditional DA schemes.

Incorporating ML into DA algorithms was suggested quite early by [47], who
clearly advocated the use of _Neural Networks_ (NNs) in a variational DA
framework. More recently this was put forward and illustrated by [1] and [9]
with a derivation of the generalised cost function from Bayes’ law, showing
that the classical ML cost function for learning a surrogate model of the
dynamical model is a limiting case of the DA framework, as in Eq. (3.117).
This formalism was furthermore generalised by [10], using classical DA
notation, for learning the state vector, the dynamical model and the error
statistics attached to this full retrieval. This section follows this latest
paper, showing how to derive the cost function of the generalised estimation
problem. Note that going even beyond this approach, attempting to learn the
optimisation scheme of the cost function, or even the full DA procedure, an
approach called _end-to-end_ in ML, is a subject of active investigations [33,
69].

###  3.5.1 Prior Error Statistics

For the sake of simplicity, again assume Gaussian statistics for the
observation errors,
$p(\boldsymbol{y}_{k}|\boldsymbol{x}_{k})=n(\boldsymbol{y}_{k}|\boldsymbol{x}_{k},\mathbf{R}_{k})$,
where the observation error covariance matrices
$\mathbf{R}_{L:1}:=\left\\{\mathbf{R}_{L},\mathbf{R}_{L-1},\cdots,\mathbf{R}_{1}\right\\}$
are supposed to be known. The dynamical model is meant to be learned or
approximated and thus stands as a surrogate model for the unknown true
physical dynamics. Assuming that the model does not explicitly depend on time,
its resolvent is defined by

$\displaystyle\boldsymbol{x}_{k}:=\mathbf{F}^{k}_{\mathbf{A}}\left(\boldsymbol{x}_{k-1}\right)+\boldsymbol{\eta}_{k},$
(3.113)

depending on a (possibly very large) set of parameters $\mathbf{A}$.
Prototypically, $\mathbf{A}$ represents the weights and biases of a NN, which
are learned from the observations alongside the state vectors within the DAW.
The distribution for model error, i.e., $\boldsymbol{\eta}_{k}$ in Eq.
(3.113), is also assumed Gaussian such that

$p(\boldsymbol{x}_{k}|\boldsymbol{x}_{k-1},\mathbf{A},\mathbf{Q}_{k}):=n\left(\boldsymbol{x}_{k}|\mathbf{F}_{\mathbf{A}}^{k}(\boldsymbol{x}_{k-1}),\mathbf{Q}_{k}\right),$
(3.114)

where
$\mathbf{Q}_{L:1}:=\left\\{\mathbf{Q}_{L},\mathbf{Q}_{L-1},\cdots,\mathbf{Q}_{1}\right\\}$
are not necessarily known. Further assume that these Gaussian errors are
white-in-time and that the observation and model errors are mutually
independent.

Note the intriguing status of $\mathbf{Q}_{L:1}$ since it depends, a
posteriori, on how well the surrogate model is estimated. This calls for an
adaptive estimation of the model error statistics $\mathbf{Q}_{L:1}$ as the
surrogate model, parametrised by $\mathbf{A}$, is better approximated.

###  3.5.2 Joint Estimation of the Model, its Error Statistics, and the State
Trajectory

Following the Bayesian formalism, one may form a MAP estimate for the joint
pdf in $\mathbf{A}$ and $\boldsymbol{x}_{L:0}$ conditioned on the
observations, together with the model error statistics. This generalised
conditional pdf is expressed in the hierarchy

$\displaystyle p(\mathbf{A},\mathbf{Q}_{L:1},\boldsymbol{x}_{L:0}|$
$\displaystyle\boldsymbol{y}_{k:1},\mathbf{R}_{K:0})=$
$\displaystyle\frac{p(\boldsymbol{y}_{L:1}|\boldsymbol{x}_{L:0},\mathbf{R}_{L:0})p(\boldsymbol{x}_{L:0}|\mathbf{A},\mathbf{Q}_{L:1})p(\mathbf{A},\mathbf{Q}_{L:1})}{p(\boldsymbol{y}_{L:1},\mathbf{R}_{L:0})},$
(3.115)

where the mutual independence of the observation and model error is used. Once
again, this remarkably stresses how powerful and general the Bayesian
framework can be.

The first term in the numerator of the right-hand-side is the usual likelihood
of the observations. The second term in the numerator is the prior on the
trajectory, given a known model and known model error statistics. The final
term of the numerator is the joint prior of the model and the model error
statistics, described as a _hyperprior_. The associated cost function is
derived proportional to

$\displaystyle\mathcal{J}_{\mathrm{DA}\text{-}\mathrm{ML}}(\mathbf{A},\boldsymbol{x}_{L:0},\mathbf{Q}_{L:1})=$
$\displaystyle-\log\,p(\mathbf{A},\mathbf{Q}_{L:1},\boldsymbol{x}_{L:0}|\boldsymbol{y}_{L:1},\mathbf{R}_{L:1})$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{k=1}^{L}\left\\{\left\|\boldsymbol{y}_{k}-\mathcal{H}_{k}(\boldsymbol{x}_{k})\right\|^{2}_{\mathbf{R}_{k}}+\log\left(\left|\mathbf{R}_{k}\right|\right)\right\\}$
$\displaystyle+\frac{1}{2}\sum_{k=1}^{L}\left\\{\left\|\boldsymbol{x}_{k}-\mathbf{F}^{k}_{\mathbf{A}}(\boldsymbol{x}_{k-1})\right\|^{2}_{\mathbf{Q}_{k}}+\log\left(\left|\mathbf{Q}_{k}\right|\right)\right\\}$
$\displaystyle-\log\,p(\boldsymbol{x}_{0},\mathbf{A},\mathbf{Q}_{L:1}).$
(3.116)

Note the resemblance of (3.116) with the weak-constraint 4D-VAR cost function
of classical DA [88]. Very importantly, this Bayesian formulation allows for a
rigorous treatment of partial and noisy observations. The classical ML cost
function that uses noiseless, complete observations of the physical system is
derived from Eq. (3.116), assuming that $\mathbf{Q}_{k}$ is known,
$\mathbf{H}_{k}\equiv\mathbf{I}_{N_{x}}$, and setting $\mathbf{R}_{k}$ to go
to $\mathbf{0}$. Associating the initial data as
$\boldsymbol{y}_{0}\longleftarrow\boldsymbol{x}_{0}$,
$\mathcal{J}_{\mathrm{DA}\text{-}\mathrm{ML}}(\mathbf{A},\boldsymbol{x}_{L:0},\mathbf{Q}_{L:1})$
becomes in its limit

$\displaystyle\mathcal{J}_{\mathrm{ML}}(\mathbf{A})=\frac{1}{2}\sum_{k=1}^{L}\left\|\boldsymbol{y}_{k}-\mathbf{F}^{k}_{\mathbf{A}}(\boldsymbol{y}_{k-1})\right\|^{2}_{\mathbf{Q}_{k}}-\log
p(\boldsymbol{y}_{0},\mathbf{A},\mathbf{Q}_{L:1}).$ (3.117)

Such connections between ML and DA were first highlighted by [47, 1, 9].

Solving the combined DA/ML problem, i.e. minimising Eq. (3.116) leads to
several key remarks. First, as mentioned earlier, this formalism allows one to
learn a surrogate model of the true dynamics using partial, possibly sparse
and noisy observations, as opposed to off-the-shelf ML techniques. This is
obviously critical for geophysical systems. Second, in this framework, one
fundamentally looks for a stochastically additive surrogate model of the form
in Eq. (3.113) rather than a deterministic surrogate model, since
$\boldsymbol{\eta}_{k}$ is drawn from the normal distribution of covariance
matrices $\mathbf{Q}_{k}$. Third, there are many ways to carry out this
minimisation, as discussed by [9, 10]. Because the state vectors and the model
parameters $\mathbf{A}$ are of fundamentally different nature, and yet are
statistically interdependent, one idea is to minimise Eq. (3.116) through a
coordinate descent, i.e. alternating minimisations on $\mathbf{A}$ and
$\boldsymbol{x}_{L:0}$, as illustrated by Fig. 3.5. This was first suggested
and successfully implemented by [16], by using an EnKF for the assimilation
step (i.e. minimising on $\boldsymbol{x}_{L:0}$), and a _Deep Learning_ (DL)
optimiser for the ML step (i.e., minimising on $\mathbf{A}$). This work was
extended to 4D variational analysis using a 4D-VAR assimilation step with a DL
optimiser for the ML step by [35].

Figure 3.5:  Estimation of both model and state trajectory using coordinate
descent by alternately optimising on the state trajectory using DA and on the
NN model parameters using ML/DL. The iterative loop stops when an accuracy
criterion is met. Source: [35]

###  3.5.3 Joint Estimation of the Model and the Error Statistics

A slightly different objective is to obtain a MAP estimate for the surrogate
model, irrespective of any model state realisation, i.e., one is interested in
the MAP of the marginal conditional pdf

$p(\mathbf{A},\mathbf{Q}_{L:1}|\boldsymbol{y}_{L:1},\mathbf{R}_{L:1})=\int\\!p(\mathbf{A},\mathbf{Q}_{L:1},\boldsymbol{x}_{L:0}|\boldsymbol{y}_{L:1},\mathbf{R}_{L:1})\,\mathrm{d}\boldsymbol{x}_{L:0},$
(3.118)

which is theoretically obtained by minimising

$\mathcal{J}(\mathbf{A},\mathbf{Q}_{L:1})=-\log\,p(\mathbf{A},\mathbf{Q}_{L:1}|\boldsymbol{y}_{L:1},\mathbf{R}_{L:1}).$
(3.119)

As pointed out by [9], the marginal pdf (3.118) can be approximately related
to the joint pdf through a Laplace approximation of the integral. Here,
however, one is interested in the full solution to this problem. This can be
solved numerically by using the _Expectation-Maximisation_ (EM) statistical
algorithm [27] jointly with the variational numerical solution of Eq. (3.116).
This was suggested in [37, 65] and implemented and validated by [10].

##  3.6 Conclusions

Unifying techniques from statistical estimation, nonlinear optimisation and
even machine learning, the Bayesian approach to DA provides a consistent
treatment of a variety of topics in the classical DA problem, discussed
throughout this tutorial. Furthermore, the Bayesian approach can be reasonably
extended to treat additional challenges not fully considered in this tutorial,
such as estimating process model parameters, handling significant modelling
errors, and including the optimisation of various hyper-parameters of the
observation-analysis-forecast cycle. Estimation may include even learning the
dynamical process model itself, as in the emerging topic of DA/ML hybrid
algorithms. In this context especially, a Bayesian analysis provides a
coherent treatment of the problem where either DA or ML techniques themselves
may be insufficient. This is particularly relevant where surrogate ML models
are used to augment traditional physics-based dynamical models for, e.g.,
simulating unresolved dynamics at scales too fine to be dynamically
represented. The hybrid approach in the Bayesian analysis provides a means to
combine the dynamical and surrogate simulations with real-world observations,
and to produce an analysis of the state and parameters for which subsequent
simulations depend. This tutorial thus presents one framework for interpreting
the DA problem, both in its classical formulation and in the directions of the
current state-of-the-art. In surveying a variety of widely used DA schemes,
the key message of this tutorial is how the Bayesian analysis provides a
consistent framework for the estimation problem and how this allows one to
formulate its solution in a variety of ways to exploit the operational
challenges in the geosciences.

## Acknowledgements

This tutorial benefited from the opportunity to lecture on data assimilation
at: (i) the CIMPA Research School on _Mathematics of Climate Science_ in
Kigali, Rawanda (28 June - 9 July, 2021); (ii) Utrecht University’s Research
school on _Data science and beyond: data assimilation with elements of machine
learning_ in Utrecht, Netherlands (23 August - 27 August, 2021); (iii) the
Joint ICTP-IUGG Workshop on _Data Assimilation and Inverse Problems in
Geophysical Sciences_ in Trieste, Italy (18 October - 29 October 2021); and
(iv) at the University of Nevada, Reno, in Fall Semester 2021, where the
combined lecture notes formed the basis of this manuscript. The authors would
like to thank everyone, but especially the students, the named reviewer Polly
Smith, and the editor Alik Ismail-Zadeh, who supported the development of
these notes with their many and highly relevant suggestions. Finally, this
tutorial also benefited from many fruitful discussions with Alberto Carrassi,
Alban Farchi, Patrick Raanes, Julien Brajard, and Laurent Bertino. CEREA is a
member of Institut Pierre-Simon Laplace (IPSL).

## Bibliography

  * [1] H. D. I. Abarbanel, P. J. Rozdeba and S. Shirman  “Machine Learning: Deepest Learning as Statistical Data Assimilation Problems”  In _Neural Computation_ 30, 2018, pp. 2025–2055  DOI: 10.1162/neco“˙a“˙01094
  * [2] B. Ait-El-Fquih and I. Hoteit  “Filtering with One-Step-Ahead Smoothing for Efficient Data Assimilation”  In _Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications (Vol. IV)_ Springer, Cham, 2022, pp. 69–96  DOI: 10.1007/978-3-030-77722-7“˙1
  * [3] B. D. O. Anderson and J. B. Moore  “Optimal Filtering”  Courier Corporation, 2012
  * [4] M. Asch, M. Bocquet and M. Nodet  “Data Assimilation: Methods, Algorithms, and Applications”, Fundamentals of Algorithms  SIAM, Philadelphia, 2016, pp. 324  DOI: 10.1137/1.9781611974546
  * [5] R. N. Bannister  “A review of operational methods of variational and ensemble-variational data assimilation”  In _Q. J. R. Meteorol. Soc._ 143, 2017, pp. 607–633  DOI: 10.1002/qj.2982
  * [6] A. G. Baydin, B. A. Pearlmutter, A. A. Radul and J. M. Siskind  “Automatic differentiation in machine learning: a survey”  In _Journal of machine learning research_ 18, 2018, pp. 1–43  DOI: 0.5555/3122009.3242010
  * [7] C. H. Bishop, B. J. Etherton and S. J. Majumdar  “Adaptive Sampling with the Ensemble Transform Kalman Filter. Part I: Theoretical Aspects”  In _Mon. Wea. Rev._ 129, 2001, pp. 420–436  DOI: 10.1175/1520-0493(2001)129¡0420:ASWTET¿2.0.CO;2
  * [8] M. Bocquet  “Ensemble Kalman filtering without the intrinsic need for inflation”  In _Nonlin. Processes Geophys._ 18, 2011, pp. 735–750  DOI: 10.5194/npg-18-735-2011
  * [9] M. Bocquet, J. Brajard, A. Carrassi and L. Bertino  “Data assimilation as a learning tool to infer ordinary differential equation representations of dynamical models”  In _Nonlin. Processes Geophys._ 26, 2019, pp. 143–162  DOI: 10.5194/npg-26-143-2019
  * [10] M. Bocquet, J. Brajard, A. Carrassi and L. Bertino  “Bayesian inference of chaotic dynamics by merging data assimilation, machine learning and expectation-maximization”  In _Foundations of Data Science_ 2, 2020, pp. 55–80  DOI: 10.3934/fods.2020004
  * [11] M. Bocquet, P. N. Raanes and A. Hannart  “Expanding the validity of the ensemble Kalman filter without the intrinsic need for inflation”  In _Nonlin. Processes Geophys._ 22, 2015, pp. 645–662  DOI: 10.5194/npg-22-645-2015
  * [12] M. Bocquet and P. Sakov  “Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems”  In _Nonlin. Processes Geophys._ 19, 2012, pp. 383–399  DOI: 10.5194/npg-19-383-2012
  * [13] M. Bocquet and P. Sakov  “Joint state and parameter estimation with an iterative ensemble Kalman smoother”  In _Nonlin. Processes Geophys._ 20, 2013, pp. 803–818  DOI: 10.5194/npg-20-803-2013
  * [14] M. Bocquet and P. Sakov  “An iterative ensemble Kalman smoother”  In _Q. J. R. Meteorol. Soc._ 140, 2014, pp. 1521–1535  DOI: 10.1002/qj.2236
  * [15] M. Bonavita, L. Isaksen and E. Hólm  “On the use of EDA background error variances in the ECMWF 4D-Var”  In _Q. J. R. Meteorol. Soc._ 138, 2012, pp. 1540–1559  DOI: 10.1002/qj.1899
  * [16] J. Brajard, A. Carrassi, M. Bocquet and L. Bertino  “Combining data assimilation and machine learning to emulate a dynamical model from sparse and noisy observations: a case study with the Lorenz 96 model”  In _J. Comput. Sci._ 44, 2020, pp. 101171  DOI: 10.1016/j.jocs.2020.101171
  * [17] G. Burgers, P. J. Leeuwen and G. Evensen  “Analysis scheme in the ensemble Kalman filter”  In _Mon. Wea. Rev._ 126, 1998, pp. 1719–1724  DOI: 10.1175/1520-0493(1998)126¡1719:ASITEK¿2.0.CO;2
  * [18] A. Carrassi, M. Bocquet, L. Bertino and G. Evensen  “Data Assimilation in the Geosciences: An overview on methods, issues, and perspectives”  In _WIREs Climate Change_ 9, 2018, pp. e535  DOI: 10.1002/wcc.535
  * [19] A. Carrassi et al.  “Data assimilation for chaotic dynamics”  In _Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications (Vol. IV)_ Springer, Cham, 2022, pp. 1–42  DOI: 10.1007/978-3-030-77722-7“˙1
  * [20] A. Carrassi, M. Bocquet, A. Hannart and M. Ghil  “Estimating model evidence using data assimilation”  In _Q. J. R. Meteorol. Soc._ 143, 2017, pp. 866–880  DOI: 10.1002/qj.2972
  * [21] Y. Chen and D. S. Oliver  “Ensemble Randomized Maximum Likelihood Method as an Iterative Ensemble Smoother”  In _Math. Geosci._ 44, 2012, pp. 1–26  DOI: 10.1007/s11004-011-9376-z
  * [22] S. E. Cohn, N. S. Sivakumaran and R. Todling  “A Fixed-Lag Kalman Smoother for Retrospective Data Assimilation”  In _Mon. Wea. Rev._ 122, 1994, pp. 2838–2867  DOI: 10.1175/1520-0493(1994)122¡2838:AFLKSF¿2.0.CO;2
  * [23] M. Corazza et al.  “Use of the breeding technique to estimate the structure of the analysis “errors of the day””  In _Nonlin. Processes Geophys._ 10, 2003, pp. 233–243  DOI: 10.5194/npg-10-233-2003
  * [24] E. Cosme et al.  “Smoothing problems in a Bayesian framework and their linear Gaussian solutions”  In _Mon. Wea. Rev._ 140, 2012, pp. 683–695  DOI: 10.1175/MWR-D-10-05025.1
  * [25] P. Courtier, J.-N. Thépaut and A. Hollingsworth  “A strategy for operational implementation of 4D-Var, using an incremental approach”  In _Q. J. R. Meteorol. Soc._ 120, 1994, pp. 1367–1387  DOI: 10.1002/qj.49712051912
  * [26] R. Daley  “Atmospheric Data Analysis”  Cambridge University Press, New-York, 1991
  * [27] A. P. Dempster, N. M. Laird and D. B. Rubin  “Maximum likelihood from incomplete data via the EM algorithm”  In _Journal of the Royal Statistical Society. Series B_ , 1977, pp. 1–38  DOI: 10.1111/j.2517-6161.1977.tb01600.x
  * [28] F. Desbouvries, Y. Petetin and B. Ait-El-Fquih  “Direct, prediction-and smoothing-based Kalman and particle filter algorithms”  In _Signal processing_ 91, 2011, pp. 2064–2077  DOI: 10.1016/j.sigpro.2011.03.013
  * [29] G. Desroziers, J.-T. Camino and L. Berre  “4DEnVar: link with 4D state formulation of variational assimilation and different possible implementations”  In _Q. J. R. Meteorol. Soc._ 140, 2014, pp. 2097–2110  DOI: 10.1002/qj.2325
  * [30] G. Evensen  “Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics”  In _J. Geophys. Res._ 99, 1994, pp. 10143–10162  DOI: 10.1029/94JC00572
  * [31] G. Evensen  “The Ensemble Kalman Filter: Theoretical Formulation and Practical Implementation”  In _Ocean Dynamics_ 53, 2003, pp. 343–367  DOI: 10.1007/s10236-003-0036-9
  * [32] G. Evensen and P. J. Leeuwen  “An Ensemble Kalman Smoother for Nonlinear Dynamics”  In _Mon. Wea. Rev._ 128, 2000, pp. 1852–1867  DOI: 10.1175/1520-0493(2000)128¡1852:AEKSFN¿2.0.CO;2
  * [33] R. Fablet et al.  “Learning Variational Data Assimilation Models and Solvers”  In _J. Adv. Model. Earth Syst._ 13, 2021, pp. e2021MS002572  DOI: 10.1029/2021MS002572
  * [34] A. Farchi and M. Bocquet  “On the efficiency of covariance localisation of the ensemble Kalman filter using augmented ensembles”  In _Front. Appl. Math. Stat._ 5, 2019, pp. 3  DOI: 10.3389/fams.2019.00003
  * [35] A. Farchi, P. Laloyaux, M. Bonavita and M. Bocquet  “Using machine learning to correct model error in data assimilation and forecast applications”  In _Q. J. R. Meteorol. Soc._ 147.739, 2021, pp. 3067–3084  DOI: https://doi.org/10.1002/qj.4116
  * [36] A. Fillion et al.  “An iterative ensemble Kalman smoother in presence of additive model error”  In _SIAM/ASA J. Uncertainty Quantification_ 8, 2020, pp. 198–228  DOI: 10.1137/19M1244147
  * [37] Z. Ghahramani and S. T. Roweis  “Learning nonlinear dynamical systems using an EM algorithm”  In _Advances in neural information processing systems_ 11 MIT Press, 1999, pp. 431–437  DOI: 0.5555/340534.340696
  * [38] M.S. Grewal and A.P. Andrews  “Kalman filtering: Theory and Practice with MATLAB”  John Wiley & Sons, 2014  DOI: 10.1002/9781118984987
  * [39] A. Griewank  “On automatic differentiation”  In _Mathematical Programming: recent developments and applications_ 6, 1989, pp. 83–107
  * [40] A. Griewank and A. Walther  “Introduction to Automatic Differentiation”  In _PAMM_ 2.1, 2003, pp. 45–49  DOI: https://doi.org/10.1002/pamm.200310012
  * [41] C. Grudzien and M. Bocquet  “A fast, single-iteration ensemble Kalman smoother for sequential data assimilation”  In _Geosci. Model Dev._ Copernicus GmbH, 2021, pp. 1–62  DOI: https://doi.org/10.5194/gmd-2021-306
  * [42] C. Grudzien, A. Carrassi and M. Bocquet  “Chaotic dynamics and the role of covariance inflation for reduced rank Kalman filters with model error”  In _Nonlin. Processes Geophys._ 25, 2018, pp. 633–648  DOI: 10.5194/npg-25-633-2018
  * [43] Y. Gu and D. S. Oliver  “An Iterative Ensemble Kalman Filter for Multiphase Fluid Flow Data Assimilation”  In _SPE Journal_ 12, 2007, pp. 438–446  DOI: 10.2118/108438-PA
  * [44] T. M. Hamill and C. Snyder  “A Hybrid Ensemble Kalman Filter-3D Variational Analysis Scheme”  In _Mon. Wea. Rev._ 128, 2000, pp. 2905–2919  DOI: 10.1175/1520-0493(2000)128¡2905:AHEKFV¿2.0.CO;2
  * [45] W.K. Härdle, O. Okhrin and Y. Okhrin  “Basic elements of computational statistics”  Springer, 2017  DOI: https://doi.org/10.1007/978-3-319-55336-8
  * [46] L. Hascoët  “An introduction to inverse modelling and parameter estimation for atmosphere and ocean sciences”  In _Advanced data assimilation for geosciences_ Les Houches school of physics: Oxford University Press, 2014, pp. 349–369  DOI: 10.1093/acprof:oso/9780198723844.001.0001
  * [47] W. W. Hsieh and B. Tang  “Applying Neural Network Models to Prediction and Data Analysis in Meteorology and Oceanography”  In _Bull. Amer. Meteor. Soc._ 79, 1998, pp. 1855–1870  DOI: 10.1175/1520-0477(1998)079¡1855:ANNMTP¿2.0.CO;2
  * [48] B. R. Hunt, E. J. Kostelich and I. Szunyogh  “Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter”  In _Physica D_ 230, 2007, pp. 112–126  DOI: 10.1016/j.physd.2006.11.008
  * [49] M. A. Iglesias, K. J. H. Law and A. M. Stuart  “Ensemble Kalman methods for inverse problems”  In _Inverse Problems_ 29, 2013, pp. 045001  DOI: 10.1088/0266-5611/29/4/045001
  * [50] T. Janjić et al.  “On the representation error in data assimilation”  In _Q. J. R. Meteorol. Soc._ 144, 2018, pp. 1257–1278  DOI: 10.1002/qj.3130
  * [51] E. Kalnay  “Atmospheric Modeling, Data Assimilation and Predictability”  Cambridge University Press, Cambridge, 2003  DOI: doi:10.1017/CBO9780511802270
  * [52] E. Kalnay and S.-C. Yang  “Accelerating the spin-up of Ensemble Kalman Filtering”  In _Q. J. R. Meteorol. Soc._ 136, 2010, pp. 1644–1651  DOI: 10.1002/qj.652
  * [53] A.N. Kolmogorov  “Foundations of the theory of probability: Second English Edition”  Courier Dover Publications, 2018
  * [54] N. B. Kovachki and A. M Stuart  “Ensemble Kalman inversion: a derivative-free technique for machine learning tasks”  In _Inverse Problems_ 35, 2019, pp. 095005  DOI: 10.1088/1361-6420/ab1c3a
  * [55] P. Laloyaux et al.  “Towards an unbiased stratospheric analysis”  In _Q. J. R. Meteorol. Soc._ 146, 2020, pp. 2392–2409  DOI: 10.1002/qj.3798
  * [56] W. G. Lawson and J. A. Hansen  “Implications of Stochastic and Determinisitic Filters as Ensemble-Based Data Assimilation Methods in Varying Regimes of Error Growth”  In _Mon. Wea. Rev._ 132, 2004, pp. 1966–1981  DOI: 10.1175/1520-0493(2004)132¡1966:IOSADF¿2.0.CO;2
  * [57] F.-X. Le Dimet and O. Talagrand  “Variational algorithms for analysis and assimilation of meteorological observations: Theoretical aspects”  In _Tellus A_ 38, 1986, pp. 97–110  DOI: 10.3402/tellusa.v38i2.11706
  * [58] C. Liu, Q. Xiao and B. Wang  “An Ensemble-Based Four-Dimensional Variational Data Assimilation Scheme. Part I: Technical Formulation and Preliminary Test”  In _Mon. Wea. Rev._ 136, 2008, pp. 3363–3373  DOI: 10.1175/2008MWR2312.1
  * [59] A. C. Lorenc  “Analysis methods for numerical weather prediction”  In _Q. J. R. Meteorol. Soc._ 112, 1986, pp. 1177–1194  DOI: 10.1002/qj.49711247414
  * [60] A. C. Lorenc  “The potential of the ensemble Kalman filter for NWP - a comparison with 4D-Var”  In _Q. J. R. Meteorol. Soc._ 129, 2003, pp. 3183–3203  DOI: 10.1256/qj.02.132
  * [61] P.C. Mahalanobis  “On the generalized distance in statistics”  In _National Institute of Science of India_ , 1936
  * [62] B. Ménétrier and T. Auligné  “An overlooked issue of variational data assimilation”  In _Mon. Wea. Rev._ 143, 2015, pp. 3925–3930  DOI: 10.1175/MWR-D-14-00404.1
  * [63] C. D. Meyer  “Matrix Analysis and Applied Linear Algebra”  Philadelphia, PA: Society for IndustrialApplied Mathematics, 2000
  * [64] R. N. Miller, M. Ghil and F. Gauthiez  “Advanced Data Assimilation in Strongly Nonlinear Dynamical Systems”  In _J. Atmos. Sci._ 51, 1994, pp. 1037–1056  DOI: 10.1175/1520-0469(1994)051¡1037:ADAISN¿2.0.CO;2
  * [65] V. D. Nguyen, S. Ouala, L. Drumetz and R. Fablet  “EM-like Learning Chaotic Dynamics from Noisy and Partial Observations”  In _arXiv preprint arXiv:1903.10335_ , 2019
  * [66] J. Nocedal and S. Wright  “Numerical optimization”  Springer Science & Business Media, 2006  DOI: https://doi.org/10.1007/978-0-387-40065-5
  * [67] Y. Pawitan  “In all likelihood: statistical modelling and inference using likelihood”  Oxford University Press, 2001
  * [68] S. G. Penny  “Mathematical foundations of hybrid data assimilation from a synchronization perspective”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 27, 2017, pp. 126801  DOI: 10.1063/1.5001819
  * [69] M. Peyron et al.  “Latent space data assimilation by using deep learning”  In _Q. J. R. Meteorol. Soc._ 147, 2021, pp. 3759–3777  DOI: 10.1002/qj.4153
  * [70] P. N. Raanes  “On the ensemble Rauch-Tung-Striebel smoother and its equivalence to the ensemble Kalman smoother”  In _Q. J. R. Meteorol. Soc._ 142, 2016, pp. 1259–1264  DOI: 10.1002/qj.2728
  * [71] P. N. Raanes, M. Bocquet and A. Carrassi  “Adaptive covariance inflation in the ensemble Kalman filter by Gaussian scale mixtures”  In _Q. J. R. Meteorol. Soc._ 145, 2019, pp. 53–75  DOI: 10.1002/qj.3386
  * [72] P. N. Raanes, A. S. Stordal and G. Evensen  “Revising the stochastic iterative ensemble smoother”  In _Nonlin. Processes Geophys._ 26, 2019, pp. 325–338  DOI: 10.5194/npg-26-325-2019
  * [73] S. Reich and C. Cotter  “Probabilistic Forecasting and Bayesian Data Assimilation”  Cambridge University Press, 2015  DOI: 10.1017/CBO9781107706804
  * [74] R. Rojas  “Neural networks: a systematic introduction”  Springer Science & Business Media, 2013  DOI: https://doi.org/10.1007/978-3-642-61068-4
  * [75] S.M. Ross  “Introduction to probability models”  Academic press, 2014  DOI: https://doi.org/10.1016/C2012-0-03564-8
  * [76] P. Sakov and L. Bertino  “Relation between two common localisation methods for the EnKF”  In _Comput. Geosci._ 15, 2011, pp. 225–237  DOI: 10.1007/s10596-010-9202-6
  * [77] P. Sakov and M. Bocquet  “Asynchronous data assimilation with the EnKF in presence of additive model error”  In _Tellus A_ 70, 2018, pp. 1414545  DOI: 10.1080/16000870.2017.1414545
  * [78] P. Sakov, J.-M. Haussaire and M. Bocquet  “An iterative ensemble Kalman filter in presence of additive model error”  In _Q. J. R. Meteorol. Soc._ 144, 2018, pp. 1297–1309  DOI: 10.1002/qj.3213
  * [79] P. Sakov and P. R. Oke  “Implications of the Form of the Ensemble Transformation in the Ensemble Square Root Filters”  In _Mon. Wea. Rev._ 136, 2008, pp. 1042–1053  DOI: 10.1175/2007MWR2021.1
  * [80] P. Sakov, D. S. Oliver and L. Bertino  “An iterative EnKF for strongly nonlinear systems”  In _Mon. Wea. Rev._ 140, 2012, pp. 1988–2004  DOI: 10.1175/MWR-D-11-00176.1
  * [81] C. Schillings and A. M. Stuart  “Convergence analysis of ensemble Kalman inversion: the linear, noisy case”  In _Applicable Analysis_ 97, 2018, pp. 107–123  DOI: 10.1080/00036811.2017.1386784
  * [82] J.M. Tabeart et al.  “The conditioning of least-squares problems in variational data assimilation”  In _Numerical Linear Algebra with Applications_ 25, 2018, pp. e2165  DOI: 10.1002/nla.2165