# Conditional Gaussian Nonlinear System: a Fast Preconditioner and a Cheap
Surrogate Model For Complex Nonlinear Systems

Nan Chen  Department of Mathematics, University of Wisconsin-Madison, Madison,
WI 53705, USA<EMAIL_ADDRESS><EMAIL_ADDRESS> Yingda Li  Department
of Mathematics, University of Wisconsin-Madison, Madison, WI 53705, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS> Honghu Liu  Department of
Mathematics, Virginia Tech, Blacksburg, VA 24061, USA<EMAIL_ADDRESS>
###### Abstract

Developing suitable approximate models for analyzing and simulating complex
nonlinear systems is practically important. This paper aims at exploring the
skill of a rich class of nonlinear stochastic models, known as the conditional
Gaussian nonlinear system (CGNS), as both a cheap surrogate model and a fast
preconditioner for facilitating many computationally challenging tasks. The
CGNS preserves the underlying physics to a large extent and can reproduce
intermittency, extreme events and other non-Gaussian features in many complex
systems arising from practical applications. Three interrelated topics are
studied. First, the closed analytic formulae of solving the conditional
statistics provide an efficient and accurate data assimilation scheme. It is
shown that the data assimilation skill of a suitable CGNS approximate forecast
model outweighs that by applying an ensemble method even to the perfect model
with strong nonlinearity, where the latter suffers from filter divergence.
Second, the CGNS allows the development of a fast algorithm for simultaneously
estimating the parameters and the unobserved variables with uncertainty
quantification in the presence of only partial observations. Utilizing an
appropriate CGNS as a preconditioner significantly reduces the computational
cost in accurately estimating the parameters in the original complex system.
Finally, the CGNS advances rapid and statistically accurate algorithms for
computing the probability density function and sampling the trajectories of
the unobserved state variables. These fast algorithms facilitate the
development of an efficient and accurate data-driven method for predicting the
linear response of the original system with respect to parameter perturbations
based on a suitable CGNS preconditioner.

##  1 Introduction

Complex nonlinear systems are ubiquitous in many areas, including geophysics,
climate science, engineering, neuroscience and material science [118, 112,
122, 110, 49]. Mathematical modeling plays an important role in characterizing
and discovering the underlying physics of these complex systems [38, 78]. With
suitable mathematical models in hand, effective parameter inference, state
estimation and data assimilation become fundamental tasks that serve as the
prerequisites for analyzing these systems [3, 61, 86, 70, 50]. Accurate
forecast of future states and successful prediction of the system response due
to external perturbations are also central topics that have many practical
implications [75, 103, 83, 116, 71].

However, there exist quite a few mathematical and computational challenges in
analyzing and simulating complex nonlinear systems. First, the intrinsic
nonlinearity in these complex nonlinear systems often triggers strongly
chaotic or turbulent behavior [105, 36, 99]. As a consequence, intermittency,
extreme events and non-Gaussian probability density functions (PDFs) are some
of the typical features in these systems [41, 117, 92, 76, 90], which impede
the use of many traditional mathematical tools to analyze the model
properties. Second, due to the nonlinear interactions between state variables
across different scales, many of these complex nonlinear systems are high
dimensional and have multiscale spatiotemporal structures [122, 79, 115, 84].
Therefore, developing new efficient numerical algorithms to accelerate the
computational efficiency becomes essential. Particularly, enhancing the
computational efficiency by reducing the complexity of these systems via
effective stochastic parameterizations is practically important and is widely
used in for example climate sciences. Third, it is often the case in practice
that only partial observations of the state variables are available [61, 69],
which result in additional difficulties for model calibration, state
estimation and prediction where systematic uncertainty quantification needs to
be addressed [32, 33, 37, 80, 81].

Since many complex dynamical systems are too expensive to be handled directly,
it is of practical importance to develop suitable approximate models, which
capture certain features of nature and are easier to deal with. Exploiting
systematic reduced order modeling strategies and effective (stochastic)
parameterizations is often a pre-requisite for the development of approximate
models. Linear regression models are arguably the simplest class of
approximate models [125, 43], which can already provide certain skill for
short-term forecasts although they usually suffer in characterizing the
underlying nonlinear physics. Physics-constrained regression models are a set
of nonlinear approximate models [87, 55, 65], which take into account the
energy conserving nonlinear interactions in the model development that
guarantees the well-posednss of long-term behavior of the system. Another
commonly used approach to developing approximate models is to project the
starting complex nonlinear system to the leading a few energetic modes in
light of the Galerkin proper orthogonal decomposition methods [59] or other
empirical basis functions such as the principal interaction patterns [56, 68]
and the dynamic mode decomposition [104, 109]. With a careful design of the
closure terms to compensate the truncation error, these reduced order models
are skillful in resolving certain problems in fluids and turbulence [8, 97,
113, 124].

Meanwhile, many data-driven approximate modeling strategies have recently been
developed [1, 16, 14, 72, 94, 102, 57, 111]. One of them is the sparse
identification of nonlinear dynamical systems (SINDy) [7], which leads to
nonlinear regression models with parsimonious structures via sparse regression
and compressed sensing. Many other approximate modeling approaches have also
been designed for specific scientific purposes. For example, the past noise
forecasting method [17] was developed as a data-driven forecast model for
stochastic climate processes that exhibit low-frequency variability. Reduced-
space Gaussian process regression forecast [120] was designed for data-driven
probabilistic forecast of chaotic dynamical systems. Small-scale
parameterization based on a data-informed optimal homotopic deformation of
invariant manifolds was developed to design low-dimensional models for both
deterministic chaotic systems and stochastic systems [15, 13]. Physically
consistent data-driven weather forecasting techniques were proposed and
applied to operational models [11, 10]. Recently, a strong link between the
stochastic parameterization approach based on perturbation expansions of the
Koopman operator [123] and the data-driven empirical model reduction (EMR)
methodology [66] was established in [107]. In addition, machine learning
methods nowadays have been extensively incorporated into the reduced order
models to further improve the approximation and forecast skill [106, 12, 19,
101, 93].

The objective of this paper is to explore the skill of a rich class of
nonlinear stochastic models, known as the “conditional Gaussian nonlinear
system” (CGNS) [20], as approximate models for complex nonlinear systems. The
CGNS includes many physics-constrained nonlinear stochastic models (e.g., the
stochastic versions of various Lorenz models, low-order models of Charney-
DeVore flows, and a paradigm model for topographic mean flow interaction),
quite a few stochastically coupled reaction-diffusion models in neuroscience
and ecology (e.g., stochastically coupled FitzHugh-Nagumo (FHN) models and
stochastically coupled SIR epidemic models), and several large-scale dynamical
models in engineering and geophysical flows (e.g., the Boussinesq equations
with noise and stochastically forced rotating shallow water equation). See the
article [20] for a gallery of examples of the CGNS. The CGNS has also been
applied to modeling and forecasting several important climate phenomena, such
as the Madden-Julian oscillation and the monsoon [25, 26], and has been
utilized for Lagrangian data assimilation [27]. Yet, most of the previous work
focused on perfect model scenarios, where utilizing the CGNS as an approximate
model has not been systematically studied.

The CGNS has several unique features that allows it to be distinct from many
existing approximate modeling strategies. First, the CGNS aims at preserving
the underlying physical mechanism to the greatest extent. Specifically, the
nonlinearity involving the large-scale or slow variables are by design
retained, which includes not only the self-interactions among the large-scale
variables but also the cross-scale interactions between large- and small-scale
variables, while suitable approximations are imposed primarily on the
nonlinear self-interactions between small-scale, fast-varying or unresolved
state variables via effective stochastic parameterizations. This is
fundamentally different from many purely data-driven nonlinear regression
models, which may miss certain crucial underlying physics of the original
complex nonlinear systems. Second, the stochastic parameterizations of the
self-interactions between small-scale variables lead to an important feature
of the CGNS. That is, the distribution of the small-scale variables
conditioned on the large-scale ones is Gaussian. One remarkable consequence is
that the associated conditional Gaussian distribution can be calculated using
closed analytic formulae [73], which considerably facilitate the mathematical
analysis and numerical simulations of the CGNS. In fact, the closed analytic
formulae of the conditional Gaussian distributions allow the development of
efficient and statistically accurate algorithms for parameter estimation, data
assimilation and ensemble forecast in light of only partial observations. Note
that, despite the conditional Gaussianity, the joint and marginal
distributions of the CGNS remain highly non-Gaussian. Thus, the intermittency,
extreme events and turbulent features can all be preserved in a suitably
designed CGNS. Third, the CGNS is also adaptable to many data-driven
scenarios. Physics constraints, localizations and sparse identification
together with many other mathematical and computational strategies can be
possibly incorporated into the CGNS. Finally, information theory [63, 77] can
be applied to quantify the uncertainty and the statistical error of the CGNS
in approximating the original complex nonlinear systems.

The specific goal of this paper is two-fold. First, by taking advantage of its
analytically solvable properties, the CGNS can be served as a fast
preconditioner for facilitating many computationally challenging tasks
associated with the original complex nonlinear system. Important applications
include estimating the parameters of the original system in the presence of
only partial observations and recovering the non-Gaussian PDFs as a crucial
intermediate step for computing the model sensitivity and response. Second,
the CGNS is exploited as a surrogate model by exploiting systematic reduced
order modeling strategies and suitable stochastic parameterizations, aiming at
spending a much lower computational cost to create comparably accurate results
as those obtained from the original complex system. This includes for example
the state estimation of unobserved variables and the statistical forecast.
Figure 1 shows a schematic illustration of utilizing the CGNS as a fast
preconditioner and a cheap surrogate model for general nonlinear system.

Figure 1: A schematic illustration of utilizing the CGNS as a fast
preconditioner and a cheap surrogate model for general nonlinear system.

The rest of the paper is organized as follows. The general mathematical
framework of CGNS and its analytic properties are described in Section 2.
Several systematic strategies for the development of the CGNS are included in
Section 3. Sections 4 – 6 consist of three important tasks in complex
nonlinear systems, showing the roles of the CGNS both as a surrogate model and
a preconditioner for the original complex system. Specifically, Section 4
focuses on the data assimilation and ensemble forecast, Section 5 aims at
efficient parameter estimation, and Section 6 illustrates the use of CGNS in
facilitating the study of model sensitivity and response theory. The paper is
concluded in Section 7.

##  2 General Mathematical Framework of the CGNS

Let us start with the general formulation of the turbulent dynamical systems
motivated from fluid and geophysical applications [118, 105, 61, 79],

$\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}=\left(L+D\right)\mathbf{u}+B\left(\mathbf{u},\mathbf{u}\right)+\mathbf{F}\left(t\right)+\bm{\sigma}\left(\mathbf{u},t\right)\dot{\mathbf{W}}\left(t\right),$
(1)

where the state variable $\mathbf{u}\in\mathbb{C}^{N}$ is in a high
dimensional phase space. In (1), the first two components,
$\left(L+D\right)\mathbf{u}$, represent linear dispersion and dissipation
effects, where $L^{*}=-L$ is a skew-symmetric operator (with $\cdot^{*}$ being
the complex conjugate transpose); and $D$ is a negative-definite matrix. The
nonlinear effect is introduced through an energy-conserving quadratic form,
$B\left(\mathbf{u},\mathbf{u}\right)$. Besides, the system is subject to
external forcing effects that are decomposed into a deterministic component,
$\mathbf{F}\left(t\right)$, and a stochastic component represented by a
Gaussian random process,
$\bm{\sigma}\left(\mathbf{u},t\right)\dot{\mathbf{W}}\left(t\right)$, where
$\bm{\sigma}\in\mathbb{C}^{N\times K}$ is the noise matrix and
$\dot{\mathbf{W}}\in\mathbb{C}^{K}$ is the white noise. The two components
$\left(L+D\right)\mathbf{u}+B\left(\mathbf{u},\mathbf{u}\right)+\mathbf{F}\left(t\right)$
and $\bm{\sigma}\left(\mathbf{u},t\right)$ on the right hand side of (1) are
also known as the drift part and the diffusion coefficients, respectively.

###  2.1 The CGNS

Despite being highly nonlinear and possessing strongly non-Gaussian statistics
in both the marginal and joint distributions of the state $\mathbf{u}$, many
complex nonlinear dynamical systems (1) have or can be approximated by the
following nonlinear system with conditional Gaussian structures. The general
mathematical framework of the CGNS is as follows [73, 21, 20],

$\displaystyle\frac{{\,\rm d}\mathbf{X}}{{\,\rm d}t}$
$\displaystyle=\Big{[}\mathbf{A}_{\mathbf{0}}(\mathbf{X},t)+\mathbf{A}_{\mathbf{1}}(\mathbf{X},t)\mathbf{Y}(t)\Big{]}+\mathbf{B}_{\mathbf{1}}(\mathbf{X},t)\dot{\mathbf{W}}_{\mathbf{1}}(t),$
(2a) $\displaystyle\frac{{\,\rm d}\mathbf{Y}}{{\,\rm d}t}$
$\displaystyle=\Big{[}\mathbf{a}_{\mathbf{0}}(\mathbf{X},t)+\mathbf{a}_{\mathbf{1}}(\mathbf{X},t)\mathbf{Y}(t)\Big{]}+\mathbf{b}_{\mathbf{2}}(\mathbf{X},t)\dot{\mathbf{W}}_{\mathbf{2}}(t),$
(2b)

where the original model state $\mathbf{u}$ is decomposed into multi-
dimensional state variables $\mathbf{X}\in\mathbb{C}^{N_{1}}$ and
$\mathbf{Y}\in\mathbb{C}^{N_{2}}$, with $N_{1}+N_{2}=N$. In (2),
$\mathbf{A}_{0},\mathbf{a}_{0},\mathbf{A}_{1},\mathbf{a}_{1},\mathbf{B}_{1}$
and $\mathbf{b}_{2}$ are vectors or matrices that can depend nonlinearly on
the state variables $\mathbf{X}$ and time $t$ while
$\dot{\mathbf{W}}_{\mathbf{1}}$ and $\dot{\mathbf{W}}_{\mathbf{2}}$ are
independent white noise sources that can have different dimensions from
$\mathbf{X}$ and $\mathbf{Y}$.

The name “conditional Gaussian” comes from the fact that once a time series of
$\mathbf{X}(s)$ for $s\leq t$ is given, then the conditional distribution
$p(\mathbf{Y}(t)|\mathbf{X}(s\leq t))$ is Gaussian. This can be seen by
noticing that, with a given $\mathbf{X}$, the process of $\mathbf{Y}$ is
linear (with respect to the variable $\mathbf{Y}$ itself since $\mathbf{X}$
has been given) with Gaussian white noises. It is worthwhile to highlight
that, from the general form of the complex turbulent system (1) to the CGNS
(2), the nonlinear interactions between $\mathbf{X}$ itself and those between
$\mathbf{X}$ and $\mathbf{Y}$ in (1) are both completely retained. The only
simplification in the CGNS is to approximate the nonlinear self-interactions
between $\mathbf{Y}$ by a combination of nonlinear functions of $\mathbf{X}$,
conditional linear functions of $\mathbf{Y}$ and effective stochastic noises.
Nevertheless, if $\mathbf{Y}$ represents small-scale or fast variables, then
such a manipulation is expected to be an effective approximation that
preserves the underlying physics to a large extent.

It should be noted that the CGNS in (2) is still highly nonlinear due to the
nonlinearity in $\mathbf{A}_{0},\mathbf{a}_{0},\mathbf{A}_{1}$ and
$\mathbf{a}_{1}$ as well as the nonlinear coupling between the latter two with
$\mathbf{Y}$. Such nonlinearities preserve the non-Gaussian statistics in (1)
and allow to reproduce many observed features in nature such as extreme events
with the more tractable conditional Gaussian structure. A gallery of examples
of the CGNS, including many physics-constrained nonlinear stochastic models,
quite a few stochastically coupled reaction-diffusion models in neuroscience
and ecology, and some large-scale dynamical models in engineering and
geophysical flows can be found in [20].

Despite being highly nonlinear and non-Gaussian, one of the important features
of the CGNS (2) is that the conditional distribution of $\mathbf{Y}$ given one
realization of the time series $\mathbf{X}$ can be solved via closed analytic
formulae. Such a unique analytic property significantly facilitates the
analysis and calculations of state estimation, data assimilation and forecast.
This feature also makes the CGNS to be quite different from general nonlinear
or non-Gaussian systems. For the latter, particle methods have to be applied
for state estimation [9, 58, 98], in which many empirical tunings have to be
involved to mitigate the numerical sampling errors.

###  2.2 Closed analytic formulae for computing the conditional statistics and
data assimilation

####  2.2.1 Nonlinear filter

For the CGNS (2), given one realization of the time series $\mathbf{X}(s)$ for
$s\in[0,t]$, the conditional distribution

$p(\mathbf{Y}(t)|\mathbf{X}(s),s\leq
t)\sim\mathcal{N}(\bm{\mu}_{\mathbf{f}}(t),\mathbf{R}_{\mathbf{f}}(t))$ (3)

becomes Gaussian, where the conditional mean $\bm{\mu}_{\mathbf{f}}$ and the
conditional covariance $\mathbf{R}_{\mathbf{f}}$ are given by the following
explicit formulae [73]

$\displaystyle\frac{{\,\rm d}\bm{\mu}_{\mathbf{f}}}{{\,\rm d}t}$
$\displaystyle=(\mathbf{a}_{\mathbf{0}}+\mathbf{a}_{\mathbf{1}}\bm{\mu}_{\mathbf{f}})+(\mathbf{R}_{\mathbf{f}}\mathbf{A}_{\mathbf{1}}^{*})(\mathbf{B}_{1}\mathbf{B}_{1}^{*})^{-1}\left(\frac{{\,\rm
d}\mathbf{X}}{{\,\rm
d}t}-(\mathbf{A}_{\mathbf{0}}+\mathbf{A}_{\mathbf{1}}\bm{\mu}_{\mathbf{f}})\right),$
(4a) $\displaystyle\frac{{\,\rm d}\mathbf{R}_{\mathbf{f}}}{{\,\rm d}t}$
$\displaystyle=\mathbf{a}_{\mathbf{1}}\mathbf{R}_{\mathbf{f}}+\mathbf{R}_{\mathbf{f}}\mathbf{a}_{\mathbf{1}}^{*}+\mathbf{b}_{2}\mathbf{b}_{2}^{*}-(\mathbf{R}_{\mathbf{f}}\mathbf{A}_{\mathbf{1}}^{*})(\mathbf{B}_{1}\mathbf{B}_{1}^{*})^{-1}(\mathbf{A}_{\mathbf{1}}\mathbf{R}_{\mathbf{f}}),$
(4b)

with $\cdot^{*}$ being the complex conjugate transpose. The subscript
‘$\mathbf{f}$’ in the conditional mean $\bm{\mu}_{\mathbf{f}}$ and conditional
covariance $\mathbf{R}_{\mathbf{f}}$ is an abbreviation for ‘filter’. The
explicit formulae in (3)-(4) correspond to the optimal nonlinear filter
solution of the state variable $\mathbf{Y}(t)$ given a realization of the
observed time series $\mathbf{X}(s)$ for $s\in[0,t]$. Thus,
$\bm{\mu}_{\mathbf{f}}$ and $\mathbf{R}_{\mathbf{f}}$ in (4) are also known as
the filter posterior mean and the filter posterior covariance. The classical
Kalman-Bucy filter [60] is the simplest special example of (4).

The closed analytic formula (4) provides an efficient algorithm for the
nonlinear data assimilation of the CGNS, which avoids using the ensemble or
particle methods that may suffer from sampling errors. In Section 4, the
closed analytic data assimilation formula (4) will be used for an accurate
state estimation of the initial value that facilitates effective ensemble
forecast. It also allows the development of an efficient Gaussian mixture
algorithm for calculating the non-Gaussian PDFs of the CGNS (see Section 2.3),
which overcomes the curse of dimensionality. Such non-Gaussian PDFs are
crucial in analyzing the model sensitivity and predicting the system response,
the details of which will be discussed in Section 6.

####  2.2.2 Nonlinear smoother

Filtering exploits the observational information up to the current time
instant for an online state estimation. On the other hand, given the
observational time series within an entire time interval, the state estimation
can become more accurate. This is known as the smoother [108].

Given one realization of the observed variable $\mathbf{X}(t)$ for
$t\in[0,T]$, the optimal smoother estimate
$p(\mathbf{Y}(t)|\mathbf{X}(s),s\in[0,T])$ of the CGNS (2) is also Gaussian
[18],

$p(\mathbf{Y}(t)|\mathbf{X}(s),s\in[0,T])\sim\mathcal{N}(\bm{\mu}_{\mathbf{s}}(t),\mathbf{R}_{\mathbf{s}}(t)),$
(5)

where the conditional mean $\bm{\mu}_{\mathbf{s}}(t)$ and conditional
covariance $\mathbf{R}_{\mathbf{s}}(t)$ of the smoother at time $t$ satisfy
the following backward equations

$\displaystyle\frac{\overleftarrow{{\,\rm d}\bm{\mu}_{\mathbf{s}}}}{{\,\rm
d}t}$
$\displaystyle=-\mathbf{a}_{\mathbf{0}}-\mathbf{a}_{\mathbf{1}}\bm{\mu}_{\mathbf{s}}+(\mathbf{b}_{2}\mathbf{b}_{2}^{*})\mathbf{R}_{\mathbf{f}}^{-1}(\bm{\mu}_{\mathbf{f}}-\bm{\mu}_{\mathbf{s}}),$
(6a) $\displaystyle\frac{\overleftarrow{{\,\rm
d}\mathbf{R}_{\mathbf{s}}}}{{\,\rm d}t}$
$\displaystyle=-(\mathbf{a}_{\mathbf{1}}+(\mathbf{b}_{2}\mathbf{b}_{2}^{*})\mathbf{R}_{\mathbf{f}}^{-1})\mathbf{R}_{\mathbf{s}}-\mathbf{R}_{\mathbf{s}}(\mathbf{a}_{\mathbf{1}}^{*}+(\mathbf{b}_{2}\mathbf{b}_{2}^{*})\mathbf{R}_{\mathbf{f}})+\mathbf{b}_{2}\mathbf{b}_{2}^{*},$
(6b)

with $\bm{\mu}_{\mathbf{f}}$ and $\mathbf{R}_{\mathbf{f}}$ given by (4). Here,
the subscript ‘$\mathbf{s}$’ in the conditional mean $\bm{\mu}_{\mathbf{s}}$
and conditional covariance $\mathbf{R}_{\mathbf{s}}$ is an abbreviation for
‘smoother’, which should not be confused by the time variable $s$ in
$\mathbf{X}(s)$. The notation $\overleftarrow{{\,\rm d}\cdot}/{\,\rm d}t$
corresponds to the negative of the usual derivative, which means that the
system (6) is solved backward over $[0,T]$ with
$(\bm{\mu}_{\mathbf{s}}(T),\mathbf{R}_{\mathbf{s}}(T))=(\bm{\mu}_{\mathbf{f}}(T),\mathbf{R}_{\mathbf{f}}(T))$
with the starting value of the nonlinear smoother
$(\bm{\mu}_{\mathbf{s}}(T),\mathbf{R}_{\mathbf{s}}(T))$ being the same as the
filter estimate $(\bm{\mu}_{\mathbf{f}}(T),\mathbf{R}_{\mathbf{f}}(T))$.

The nonlinear smoother plays an important role for an unbiased state
estimation and postprocessing of the data. It is also able to quantify the
uncertainty in the unobserved variables in the parameter estimation given only
partial observations, which will be a topic to be studied in Section 5. In
addition, the nonlinear smoother is the basis to the development of a
nonlinear sampling formula, which will be shown below and is a necessary step
in analyzing the model sensitivity and predicting the system response in
Section 6.

####  2.2.3 Nonlinear sampling formula

Associated with the nonlinear smoother, a nonlinear conditional sampling
formula can be derived. In addition to satisfying the point-wise optimal
estimate (6), i.e., a distribution formed by conditional mean and conditional
distribution at each time stamp, the conditional sampled trajectories further
take into account the path-wise temporal correlation. These sampled
trajectories in the CGNS framework can be regarded as the analogs of the
ensemble members in the ensemble Kalman smoother [40], but the former can be
obtained via a closed analytic formula.

Conditioned on one realization of the observed variable $\mathbf{X}(s)$ for
$s\in[0,T]$, the optimal strategy of sampling the trajectories associated with
the unobserved variable $\mathbf{Y}$ satisfies the following explicit formula
[24],

$\frac{\overleftarrow{{\,\rm d}\mathbf{Y}}}{{\,\rm
d}t}=\frac{\overleftarrow{{\,\rm d}\bm{\mu}_{\mathbf{s}}}}{{\,\rm
d}t}-\big{(}\mathbf{a}_{\mathbf{1}}+(\mathbf{b}_{2}\mathbf{b}_{2}^{*})\mathbf{R}_{\mathbf{f}}^{-1}\big{)}(\mathbf{Y}-\bm{\mu}_{\mathbf{s}})+\mathbf{b}_{2}\dot{\mathbf{W}}_{\mathbf{Y}}(t),$
(7)

where $\dot{\mathbf{W}}_{\mathbf{Y}}(t)$ is a random noise that is independent
from $\dot{\mathbf{W}}_{2}(t)$ in (2). The conditional sampling formula is
another necessary component in analyzing the model sensitivity and predicting
the system response in Section 6.

###  2.3 Semi-analytic and statistically accurate formulae for solving the
non-Gaussian PDFs via mixtures

The closed analytic formula (4) in calculating the conditional distribution
$p(\mathbf{Y}(t)|\mathbf{X}(s),s\leq t)$ in (3) also provides an extremely
useful way to compute the marginal distribution $p({\mathbf{Y}(t)})$. In fact,
assuming there are $L$ trajectories of $\mathbf{X}(s\leq t)$, denoted by
$\mathbf{X}^{\text{obs}}_{l}(s\leq t)$ for $l=1,\ldots,L$, then in the limit
with $L\to\infty$, the marginal distribution $p({\mathbf{Y}(t)})$ is given by

$p({\mathbf{Y}(t)})=\lim_{L\to\infty}\frac{1}{L}\sum_{l=1}^{L}p\Big{(}{\mathbf{Y}(t)}|{\mathbf{X}^{\text{obs}}_{l}(s\leq
t)}\Big{)}.$ (8)

While the above identity is in the asymptotic form with $L\rightarrow\infty$,
it has been shown that the error bound in approximating $p({\mathbf{Y}(t)})$
with a finite $L$ does not depend on the dimension of $\mathbf{Y}$ [28]. In
other words, fundamentally different from the traditional Monte Carlo
simulations, the method in (8) avoids the curse of dimensionality. If it is
further assumed that the dimension of $\mathbf{X}$ is low, then the following
efficient and statistically accurate approach can be utilized to compute the
joint PDF at any transient phase $p(\mathbf{X}(t),\mathbf{Y}(t))$ [23, 22],

$p(\mathbf{X}(t),\mathbf{Y}(t))=\lim_{L\to\infty}\frac{1}{L}\sum_{l=1}^{L}\Big{(}K_{\mathbf{H}}(\mathbf{X}(t)-\mathbf{X}^{\text{obs}}_{l}(t))\,p(\mathbf{Y}(t)|\mathbf{X}^{\text{obs}}_{l}(s\leq
t))\Big{)}.$ (9)

Here the distribution of $\mathbf{X}$ is approximated by a mixture
distribution that is solved via a kernel density estimation. One of the
simplest choices, which is also the one used in this paper, is that each
mixture component is a Gaussian function centered at
$\mathbf{X}^{\text{obs}}_{l}(t)$. The same bandwidth $\mathbf{H}$ is utilized
for different Gaussian mixture components and it is determined via the “solve-
the-equation” method [5] that is an appropriate method for approximating non-
Gaussian distributions. In addition to the advantage of applying (9) in
solving high-dimensional PDF (especially when the dimension of $\mathbf{Y}$ is
large), the semi-analytic formula in (9) also allows a smoothed PDF, which
reduces the sampling error compared with other approaches even for systems
with moderate or low dimensions.

It is important to note that if the underlying system contains model error,
then the PDF provided by (9) is very different from the one created by simply
running the imperfect model. This is because the model error in the marginal
distribution $p({\mathbf{Y}(t)})$ is mitigated with the help from the
observations $\mathbf{X}^{\text{obs}}_{l}(s\leq t)$. Such a unique feature
plays a crucial role in improving the results for computing the model response
and sensitivity analysis, where the perfect model is seldom known in practice.
The details will be illustrated in Section 6. As a final remark, only the
equilibrium PDF is required in many applications, including the study of the
model response. Therefore, if the system is ergodic, then only a single
(sufficiently long) trajectory of $\mathbf{X}(0\leq t\leq T)$, denoted by
$\mathbf{X}^{\text{obs}}(0\leq t\leq T)$, is needed in computing the
equilibrium distribution $p_{\textrm{eq}}(\mathbf{X},\mathbf{Y})$

$p_{\textrm{eq}}(\mathbf{X},\mathbf{Y})=\lim_{J\to\infty}\frac{1}{J}\sum_{j=1}^{J}\Big{(}K_{\mathbf{H}}(\mathbf{X}-\mathbf{X}^{\text{obs}}(t_{j}))\,p(\mathbf{Y}|\mathbf{X}^{\text{obs}}(s\leq
t_{j}))\Big{)},$ (10)

where $[t_{1},\ldots,t_{J}]$ is a partition of the time interval $[T_{0},T]$
with some burn-in time $T_{0}$.

##  3 Strategies of Developing CGNS

The goals of developing approximate models for solving different problems in
practice are often distinct to each other. For example, some applications
require a skillful forecast model while others seek for a suitable model to
reproduce certain statistics. While there is not a universal criterion to
build the “optimal” CGNS as an approximate model for all applications, a few
potentially useful strategies are provided below for constructing CGNS, which
will be applied in the following sections.

###  3.1 Fast wave averaging

Recall the general form of the complex systems with quadratic nonlinearity
(1). Writing it into the form of state variables $(\mathbf{X},\mathbf{Y})$
yields

$\begin{split}\\!\\!\\!\frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}&\\!=\\!L_{11}\mathbf{X}+L_{12}\mathbf{Y}+B^{1}_{11}(\mathbf{X},\\!\mathbf{X})+B^{1}_{12}(\mathbf{X},\\!\mathbf{Y})+B^{1}_{22}(\mathbf{Y},\\!\mathbf{Y})+\mathbf{F}_{1}(t)+\bm{\sigma}_{1}(\mathbf{X},\\!\mathbf{Y},t)\dot{\mathbf{W}}_{1}(t),\\\
\\!\\!\\!\frac{\mathrm{d}\mathbf{Y}}{\mathrm{d}t}&\\!=\\!L_{21}\mathbf{X}+L_{22}\mathbf{Y}+B^{2}_{11}(\mathbf{X},\\!\mathbf{X})+B^{2}_{12}(\mathbf{X},\\!\mathbf{Y})+B^{2}_{22}(\mathbf{Y},\\!\mathbf{Y})+\mathbf{F}_{2}(t)+\bm{\sigma}_{2}(\mathbf{X},\\!\mathbf{Y},t)\dot{\mathbf{W}}_{2}(t).\end{split}$
(11)

In (11), $L_{ij}$ are constant matrices while $B^{k}_{ij}$ are vector
functions with the entries being quadratic functions of the state variables.
In some applications, there exists a scale separation of the state variables,
where $\mathbf{X}$ and $\mathbf{Y}$ represent slow and fast variables,
respectively. In such a case, it is natural to apply a fast wave average such
that the terms representing the self-interaction between the fast variables,
i.e., $B^{1}_{22}\left(\mathbf{Y},\mathbf{Y}\right)$ and
$B^{2}_{22}\left(\mathbf{Y},\mathbf{Y}\right)$, are approximated by stochastic
damping and noise [89, 88]. By further approximating the diffusion
coefficients $\bm{\sigma}_{1}\left(\mathbf{X},\mathbf{Y}\right)$ and
$\bm{\sigma}_{2}\left(\mathbf{X},\mathbf{Y}\right)$ by functions of only
$\mathbf{X}$, the resulting system becomes

$\begin{split}\frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}&=\widetilde{L}_{11}\mathbf{X}+\widetilde{L}_{12}\mathbf{Y}+B^{1}_{11}\left(\mathbf{X},\mathbf{X}\right)+B^{1}_{12}\left(\mathbf{X},\mathbf{Y}\right)+\mathbf{F}_{1}(t)+\widetilde{\bm{\sigma}}_{1}\left(\mathbf{X},t\right)\dot{\mathbf{W}}_{1}(t),\\\
\frac{\mathrm{d}\mathbf{Y}}{\mathrm{d}t}&=\widetilde{L}_{21}\mathbf{X}+\widetilde{L}_{22}\mathbf{Y}+B^{2}_{11}\left(\mathbf{X},\mathbf{X}\right)+B^{2}_{12}\left(\mathbf{X},\mathbf{Y}\right)+\mathbf{F}_{2}(t)+\widetilde{\bm{\sigma}}_{2}\left(\mathbf{X},t\right)\dot{\mathbf{W}}_{2}(t),\end{split}$
(12)

which belongs to the CGNS (2). Note that, if the scale separation is not
strong enough to apply the fast wave averaging, then
$B^{1}_{22}\left(\mathbf{Y},\mathbf{Y}\right)$ and
$B^{2}_{22}\left(\mathbf{Y},\mathbf{Y}\right)$ can be approximated by
additional closure terms [106, 95] that include nonlinear functions of
$\mathbf{X}$ and bilinear functions of $\mathbf{X}$ and $\mathbf{Y}$ to fit
the CGNS framework.

###  3.2 Stochastic parameterizations

The fast wave averaging or closure approximations are suitable approaches to
build CGNS if the starting complex nonlinear system is completely known.
However, in many practical applications, the information of the perfect model
is not entirely available. Specifically, while the large-scale dynamics of
nature is often accessible, the details of the small- or unresolved-scale
features are not fully understood in many applications. In such a situation,
suitable stochastic parameterizations can be adopted to approximate the
processes of the unobserved variables $\mathbf{Y}$ such that the feedback from
small/unresolved to large/resolved scales are well characterized and the
parameterized system follows the CGNS structure.

One of the simplest strategies is to apply a linear stochastic model with
Gaussian white noise to describe each component of the hidden processes of
$\mathbf{Y}$ while the processes of $\mathbf{X}$ remain highly nonlinear. This
parameterization strategy has been utilized in data assimilation and short-
term statistical prediction [46, 45, 6]. The model with such a simple
stochastic parameterization automatically fits the CGNS as there is no
quadratic function of $\mathbf{Y}$ involved. A more sophisticated stochastic
parameterization is to incorporate the physics-constrained nonlinear
regression model [87, 55] into the CGNS. This allows a more accurate way in
characterizing the nonlinear dynamics of the small-scale features in
$\mathbf{Y}$, influenced by the large-scale variables $\mathbf{X}$. In
addition, the coupled system with physics constraints also prevents finite-
time blow up of the solutions and facilitates a skillful medium- to long-range
forecast.

###  3.3 System augmentation

Another strategy to derive approximate CGNS is via a simple system
augmentation technique detailed below. As has been discussed in Section 3.1
that the most significant difference between the general nonlinear system (1)
and the CGNS (2) is the quadratic nonlinear self-interactions of $\mathbf{Y}$,
namely the terms $B^{1}_{22}\left(\mathbf{Y},\mathbf{Y}\right)$ and
$B^{2}_{22}\left(\mathbf{Y},\mathbf{Y}\right)$ appearing in (11). After
suitable regrouping of the terms in (11) and assuming for simplicity that the
diffusion coefficients depend only on $X$, we can rewrite this system (11)
into the following form,

$\displaystyle\frac{{\,\rm d}\mathbf{X}}{{\,\rm d}t}$
$\displaystyle=\Big{[}\mathbf{A}_{\mathbf{0}}(\mathbf{X},t)+\mathbf{A}_{\mathbf{1}}(\mathbf{X},t)\mathbf{Y}+B^{1}_{22}(\mathbf{Y},\mathbf{Y})\Big{]}+\bm{\sigma}_{1}(\mathbf{X},t)\dot{\mathbf{W}}_{\mathbf{1}}(t),$
(13a) $\displaystyle\frac{{\,\rm d}\mathbf{Y}}{{\,\rm d}t}$
$\displaystyle=\Big{[}\mathbf{a}_{\mathbf{0}}(\mathbf{X},t)+\mathbf{a}_{\mathbf{1}}(\mathbf{X},t)\mathbf{Y}+B^{2}_{22}(\mathbf{Y},\mathbf{Y})\Big{]}+\bm{\sigma}_{2}(\mathbf{X},t)\dot{\mathbf{W}}_{\mathbf{2}}(t),$
(13b)

which differs from the CGNS (2) by the two quadratic terms
$B^{1}_{22}(\mathbf{Y},\mathbf{Y})$ and $B^{2}_{22}(\mathbf{Y},\mathbf{Y})$.
Instead of approximating these two terms directly via fast wave averaging as
proposed in Section 3.1, we consider here the situation that the scale
separation between $\mathbf{X}$ and $\mathbf{Y}$ is not pronounced.

We will still handle $B^{2}_{22}(\mathbf{Y},\mathbf{Y})$ in (13b) via suitable
stochastic parameterization in terms of $\mathbf{X}$, leading to an
approximate equation for $\mathbf{Y}$ of the form

$\frac{{\,\rm d}\mathbf{Y}}{{\,\rm
d}t}=\Big{[}\widetilde{\mathbf{a}}_{\mathbf{0}}(\mathbf{X},t)+\widetilde{\mathbf{a}}_{\mathbf{1}}(\mathbf{X},t)\mathbf{Y}\Big{]}+\widetilde{\bm{\sigma}}_{2}(\mathbf{X},t)\dot{\mathbf{W}}_{\mathbf{2}}(t).$
(14)

The term $B^{1}_{22}(\mathbf{Y},\mathbf{Y})$ is then dealt with through a
system augmentation strategy as explained below. Denote by
$\mathbf{Z}=((Y_{1})^{2},Y_{1}Y_{2},Y_{1}Y_{3},\ldots)^{\top}$ the
$N_{2}(N_{2}+1)/2$-dimensional vector whose components consist of all possible
quadratic monomials involving the components of the $N_{2}$-dimensional small-
scale variable $\mathbf{Y}$. Using Itô’s formula [44] and (14), we can obtain
the corresponding equation for $\mathbf{Z}$, which takes the following form

$\begin{gathered}\frac{{\,\rm d}\mathbf{Z}}{{\,\rm
d}t}=\Big{[}\mathbf{c}_{\mathbf{0}}(\mathbf{X},t)+\mathbf{c}_{\mathbf{1}}(\mathbf{X},t)\mathbf{Y}+\mathbf{c}_{\mathbf{2}}(\mathbf{X},t)\mathbf{Z}\Big{]}+\bm{\sigma}_{3}(\mathbf{X},\mathbf{Y},t)\dot{\mathbf{W}}_{\mathbf{2}}(t).\end{gathered}$
(15)

Note that the quadratic term $B^{1}_{22}(\mathbf{Y},\mathbf{Y})$ in (13a) can
be rewritten as a linear function of $\mathbf{Z}$, denoted by $L\mathbf{Z}$,
thanks to the very choice of $\mathbf{Z}$. Thus, if we further approximate
$\mathbf{Y}$ in the diffusion coefficient
$\bm{\sigma}_{3}(\mathbf{X},\mathbf{Y},t)$ of (15) by, e.g., its global mean
$\overline{\mathbf{Y}}$, the augmented system for
$(\mathbf{X},\mathbf{Y},\mathbf{Z})$ consisting of (13a), (14), and (15) fits
into the form of CGNS given by (2) with $(\mathbf{Y},\mathbf{Z})$ here playing
the role of $\mathbf{Y}$ in (2), after replacing
$B^{1}_{22}(\mathbf{Y},\mathbf{Y})$ in (13a) by $L\mathbf{Z}$ and
approximating $\bm{\sigma}_{3}(\mathbf{X},\mathbf{Y},t)$ in (15) by
$\bm{\sigma}_{3}(\mathbf{X},\overline{\mathbf{Y}},t)$.

With the CGNS structure available through this augmented system, the
conditional distribution such as $p(\mathbf{Y}(t)|\mathbf{X}(s\leq t))$ can
now be approximated by first computing the conditional distribution
$p(\mathbf{Y}(t),\mathbf{Z}(t)|\mathbf{X}(s\leq t))$ for the augmented system
using formulae provided in Section 2.2 and then marginalizing it over
$\mathbf{Z}$. The essence of this simple approach is thus to increase the
dimension of the system in exchange of analytic formulae of the conditional
distributions, which would otherwise be computationally expensive to compute
even for models with moderate dimensions.

When the dimension, $N_{2}$, of $\mathbf{Y}$ is very high, one may prefer to
identify only a subset of $\mathbf{Y}$ to apply the system augmentation
technique in order not to inflate too much the number of variables, since the
auxiliary variable $\mathbf{Z}$ has dimension $N_{2}(N_{2}+1)/2$. Such an
extension goes beyond the scope of the current article, and will be addressed
in a separate communication. Note also that if
$B^{2}_{22}(\mathbf{Y},\mathbf{Y})\equiv 0$, then Eq. (14) is reduced to
(13b). In this case, there is no approximation involved in the drift part of
(15).

We demonstrate now the efficiency of this approach in the context of data
assimilation and ensemble forecast using a low-dimensional truncation of a
stochastic Burgers-type equation.

##  4 Data Assimilation and Ensemble Forecast

Data assimilation concerns the problem of estimating the state variables of a
given, usually nonlinear and possibly stochastic, dynamical system when
observations of certain related output variables are available [40, 61, 70,
86]. One major challenge in data assimilation is the strong nonlinearity and
the associated non-Gaussian statistics in the underlying dynamics, in which a
direct application of the particle methods may be inaccurate especially in the
high dimensional situations. The development of cheap and effective
approximate models that capture the main characteristics of the underlying
dynamics is thus an important topic in state estimation and data assimilation.
Since the data assimilation solution corresponds to the initialization of the
subsequent forecast, an efficient and accurate data assimilation scheme is
also essential to advancing the forecast skill. Note that there is usually a
stronger demand in developing suitable approximate models for data
assimilation than the subsequent short- or medium-range forecast since the
former often involves many numerical or sampling issues in the presence of
strong nonlinearity and non-Gaussianity.

In this section, a particular type of approximate models for this purpose
obtained by the method of system augmentations presented in Section 3.3 is
studied. The resulting approximate model has the form of a CGNS. Thus, the
associated data assimilation solutions can be calculated using the closed
analytic formula (4) as was discussed in Section 2.2. To simplify the
presentation, the idea is illustrated using a low-dimensional stochastic
differential equation (SDE) with energy-conserving quadratic terms. The data
assimilation results from the CGNS, which is an approximate model, will be
compared with that by applying a classical ensemble data assimilation method
directly to the perfect model. The goal is to illustrate the efficiency and
accuracy of the data assimilation scheme using the CGNS, especially in
avoiding the potential sampling and other numerical issues that appear in the
ensemble-based approaches.

###  4.1 A truncated stochastic quadratic system and its CGNS approximation
through system augmentation

The model considered here is the following three-dimensional SDEs with energy-
conserving quadratic nonlinear terms and subject to additive white noise
forcing

$\displaystyle\frac{{\,\rm d}x}{{\,\rm d}t}=\beta_{x}x+\alpha xy+\alpha
yz+\sigma_{x}\dot{W}_{x},$ (16a) $\displaystyle\frac{{\,\rm d}y}{{\,\rm
d}t}=\beta_{y}y-\alpha x^{2}+2\alpha xz+\sigma_{y}\dot{W}_{y},$ (16b)
$\displaystyle\frac{{\,\rm d}z}{{\,\rm d}t}=\beta_{z}z-3\alpha
xy+\sigma_{z}\dot{W}_{z}.$ (16c)

Here, the coefficients for the linear terms are chosen such that $\beta_{x}$
is positive to introduce linear instability into the system, while $\beta_{y}$
and $\beta_{z}$ are negative, representing linear damping effects. The
coefficient $\alpha>0$ controls the strength of the nonlinearity; and the
noise strength coefficients $\sigma_{x}$, $\sigma_{y}$, and $\sigma_{z}$ are
positive constants. This system can for instance be obtained as a Fourier-
Galerkin projection of the stochastic Burgers-Sivashinsky equation

$\frac{\partial u}{\partial t}=\big{(}\nu\partial_{xx}u+\lambda
u-u\partial_{x}u\big{)}+\dot{W}(t,x)$

posed on a bounded interval $x\in(0,L)$ subject to homogeneous Dirichlet
boundary conditions. In this context, $\beta_{x}$, $\beta_{y}$, and
$\beta_{z}$ are simply the three largest eigenvalues of the linear operator
and $\alpha$ is linked to the domain size $L$ via
$\alpha=\pi/(\sqrt{2}L^{3/2})$; see e.g. [15, Chapter 6].

In the following, the largest scale variable $x$ is treated as the observed
variable while there is no direct observations for the state variables
$(y,z)$. Under this splitting of the state variables, system (16) does not
have the conditional Gaussian structure due to the quadratic nonlinear term
$\alpha yz$ between the unobserved variables that appears in (16a). Following
the idea presented in Section 3.3, in order to obtain a CGNS to approximate
the system (16), three auxiliary variables are introduced for the possible
quadratic interactions between the two unobserved variables:

$p=y^{2},\quad q=yz,\quad r=z^{2}.$ (17)

Using (16) and apply Itô’s formula yields

$\begin{split}&\frac{{\,\rm d}p}{{\,\rm
d}t}=(\sigma_{y})^{2}+2\big{(}\beta_{y}p-\alpha x^{2}y+2\alpha
xq\big{)}+2\sigma_{y}y\dot{W}_{y},\\\ &\frac{{\,\rm d}q}{{\,\rm
d}t}=(\beta_{y}+\beta_{z})q-\alpha x^{2}z-3\alpha xp+2\alpha
xr+\sigma_{y}z\dot{W}_{y}+\sigma_{y}y\dot{W}_{z},\\\ &\frac{{\,\rm d}r}{{\,\rm
d}t}=(\sigma_{z})^{2}+2\big{(}\beta_{z}r-3\alpha
xq\big{)}+2\sigma_{z}z\dot{W}_{z}.\end{split}$ (18)

Assume that the global mean values of the unobserved variables $y$ and $z$ are
accessible (from a period of training data). Then, in combination with (16),
and after replacing $y$ and $z$ in the state-dependent noise terms of (18) by
their respective global mean, the following augmented system is arrived at:

$\begin{split}&\frac{{\,\rm d}x}{{\,\rm d}t}=\beta_{x}x+\alpha xy+\alpha
q+\sigma_{x}\dot{W}_{x},\\\ &\frac{{\,\rm d}y}{{\,\rm d}t}=\beta_{y}y-\alpha
x^{2}+2\alpha xz+\sigma_{y}\dot{W}_{y},\\\ &\frac{{\,\rm d}z}{{\,\rm
d}t}=\beta_{z}z-3\alpha xy+\sigma_{z}\dot{W}_{z},\\\ &\frac{{\,\rm d}p}{{\,\rm
d}t}=(\sigma_{y})^{2}+2\big{(}\beta_{y}p-\alpha x^{2}y+2\alpha
xq\big{)}+2\sigma_{y}\overline{y}\dot{W}_{y},\\\ &\frac{{\,\rm d}q}{{\,\rm
d}t}=(\beta_{y}+\beta_{z})q-\alpha x^{2}z-3\alpha xp+2\alpha
xr+\sigma_{y}\overline{z}\dot{W}_{y}+\sigma_{z}\overline{y}\dot{W}_{z},\\\
&\frac{{\,\rm d}r}{{\,\rm d}t}=(\sigma_{z})^{2}+2\big{(}\beta_{z}r-3\alpha
xq\big{)}+2\sigma_{z}\overline{z}\dot{W}_{z},\end{split}$ (19)

where $yz$ in (16a) becomes the state variable $q$. This augmented system (19)
fits into the CGNS form of (2) with now the unobserved variables taken to be
$\mathbf{Y}=(y,z,p,q,r)^{\top}$.

Although the dimension of the approximate system is increased compared with
the original system, closed analytic equations are now accessible for the
evolution of the corresponding conditional statistics for the data
assimilation solutions (see equation (4) in Section 2.2). As will be shown
below, the approximate system (19) can provide a significantly more accurate
estimation of $(y,z)$ compared with another conditional Gaussian approximation
obtained by simply removing the term $\alpha yz$ in (16a), called the bare
truncation (BT) system below. The skill of the proposed method is comparable
and sometimes even more accurate than the ensemble Kalman-Bucy filter (EnKBF)
[4], while being more efficient thanks to the availability of analytic
formulae.

###  4.2 Dynamical regimes and numerical setup

In the following, we consider two dynamical regimes:

$\begin{split}\text{Regime
I:}&\quad\sigma_{x}=1,\sigma_{y}=1,\sigma_{z}=2,\beta_{x}=0.1,\beta_{y}=-0.5,\beta_{z}=-1,\alpha=\pi/\sqrt{2},\\\
\text{Regime
II:}&\quad\sigma_{x}=0.1,\sigma_{y}=1,\sigma_{z}=2,\beta_{x}=0.1,\beta_{y}=-0.5,\beta_{z}=-1,\alpha=\pi/\sqrt{2}.\end{split}$
(20)

In particular, we have a relatively strong nonlinear effects with
$\alpha\approx 2.2$, and a relatively small spectral gap between the observed
and the hidden variables with $\beta_{x}-\beta_{y}=0.6$. Both of the two
hidden variables are subject to strong noise perturbations. The two regimes
differ only in the value of the noise strength $\sigma_{x}$ in the
$x$-equation.

The same numerical setup is adopted for both parameter regimes. The true
signal is obtained by integrating the original SDE system (16) for an
arbitrarily fixed noise path using the Euler-Maruyama scheme with a uniform
time step size $\delta t=5\times 10^{-4}$ and initialized at
$(x,y,z)=(0,0,0)$.

For the state estimation of the unobserved variables $(y,z)$, we compare three
methods:

  * Method 1: Apply the nonlinear filtering formulae (4) for the general CGNS (2) to the augmented system (19), with $\mathbf{X}=x$ and $\mathbf{Y}=(y,z,p,q,r)^{\top}$. This method will be referred as the CG method below.

  * Method 2: Apply the nonlinear filtering formulae (4) to a bare truncation of (16) in which we simply remove the term $\alpha yz$ in (16a) to obtain a CGNS, with $\mathbf{X}=x$ and $\mathbf{Y}=(y,z)^{\top}$. This method will be referred as the BT method below.

  * Method 3: Apply the ensemble Kalman-Bucy filtering (EnKBF) method to (16). See (24) below for its formulation.

The data assimilation for each of the above methods is performed over the time
window $[0,400]$ with the same time step size $\delta t$ as the true signal.
For the CG method, the global mean values $\overline{y}$ and $\overline{z}$ in
(19) are taken to be the mean values of the corresponding true signal over the
interval $[0,200]$. For both CG and BT, the initial values of the conditional
mean and conditional covariance are taken to be zero. For EnKBF, the size of
ensemble is taken to be $N=100$ and the unobserved variables are initialized
at $(y,z)=(0,0)$. For the sake of clarity, we provide below some details about
the EnKBF applied to (16). We introduce the following notations for the drift
part of the system (16):

$\begin{split}&g(x,y,z)=\beta_{x}x+\alpha xy+\alpha yz,\\\
&f_{1}(x,y,z)=\beta_{y}y-\alpha x^{2}+2\alpha xz,\\\
&f_{2}(x,y,z)=\beta_{z}z-3\alpha xy.\end{split}$ (21)

Denote by $\mathbf{y}=(y_{1},y_{2},\ldots,y_{N})^{\top}$ and
$\mathbf{z}=(z_{1},z_{2},\ldots,z_{N})^{\top}$ the collection of all the $N$
ensemble members. We define also

$\begin{split}&\mathcal{N}_{1}(x_{\text{obs}}(t),\mathbf{y},\mathbf{z})=\frac{1}{\sigma_{x}^{2}(N-1)}\sum_{j=1}^{N}(y_{j}-\overline{\mathbf{y}}(t))(g(x_{\text{obs}}(t),y_{j},z_{j})-\overline{g}(x_{\text{obs}}(t),\mathbf{y},\mathbf{z})),\\\
&\mathcal{N}_{2}(x_{\text{obs}}(t),\mathbf{y},\mathbf{z})=\frac{1}{\sigma_{x}^{2}(N-1)}\sum_{j=1}^{N}(z_{j}-\overline{\mathbf{z}}(t))(g(x_{\text{obs}}(t),y_{j},z_{j})-\overline{g}(x_{\text{obs}}(t),\mathbf{y},\mathbf{z})),\\\
\end{split}$ (22)

where

$\overline{\mathbf{y}}(t)=\frac{1}{N}\sum_{\ell=1}^{N}y_{\ell}(t),\quad\overline{\mathbf{z}}(t)=\frac{1}{N}\sum_{\ell=1}^{N}z_{\ell}(t),\quad\overline{g}(x_{\text{obs}}(t),\mathbf{y},\mathbf{z})=\frac{1}{N}\sum_{\ell=1}^{N}g(x_{\text{obs}}(t),y_{\ell},z_{\ell}).$
(23)

Then, each ensemble member $(y_{i},z_{i})$, $i=1,2,\ldots,N$, of the EnKBF is
computed using

$\begin{split}\frac{{\,\rm d}y_{i}}{{\,\rm
d}t}&=f_{1}(x_{\text{obs}}(t),y_{i},z_{i})+\sigma_{y}\dot{W}_{y,i}\\\
&\qquad\qquad\qquad-\mathcal{N}_{1}(x_{\text{obs}}(t),\mathbf{y},\mathbf{z})\big{[}g(x_{\text{obs}}(t),y_{i},z_{i})-\dot{x}_{\text{obs}}(t)+\sigma_{x}\dot{W}_{x,i}\big{]},\\\
\frac{{\,\rm d}z_{i}}{{\,\rm
d}t}&=f_{2}(x_{\text{obs}}(t),y_{i},z_{i})+\sigma_{z}\dot{W}_{z,i}\\\
&\qquad\qquad\qquad-\mathcal{N}_{2}(x_{\text{obs}}(t),\mathbf{y},\mathbf{z})\big{[}g(x_{\text{obs}}(t),y_{i},z_{i})-\dot{x}_{\text{obs}}(t)+\sigma_{x}\dot{W}_{x,i}\big{]},\end{split}$
(24)

where $W_{x,i}$, $W_{y,i}$, and $W_{z,i}$, $i=1,2,\ldots,N$, are all mutually
independent one-dimensional Brownian motions, and $x_{\text{obs}}$ is the
observed signal of $x$.

For Regime II, we will also compare the ensemble forecast skills. The forecast
is performed over the time window $[200,400]$, which is chosen to avoid
overlap with the training window $[0,200]$ from which the global mean values
of $y$ and $z$ appearing in (19) are computed. The forecast model is taken to
be the true SDE system (16), and the initial conditions (IC) of $(y,z)$ are
drawn from multivariate Gaussian distributions with mean and covariance
estimated respectively from BT, CG and EnKBF described above. For $x$, its
initial value is taken to be that of the true signal at the corresponding time
instant. We will also compute the results when the forecast is initialized
with the true signal for all the three variables, which serves as the
reference of the theoretic forecast/predictability limit and will be referred
as the case with the perfect IC. The time locations at which to issue the
forecasts are equally spaced over the chosen time interval, with a gap of 0.01
between two adjacent forecasts, leading thus to a total of $2\times 10^{4}$
forecast locations. Each forecast is computed up to a lead time of 1 time
unit, and a total of 40 ensemble members are generated at each forecast
location. This procedure is repeated for each of the methods used to construct
the IC.

###  4.3 Numerical results

We present now the results obtained based on the numerical procedure described
above. For the two regimes given by (20), due to the larger noise strength
parameter $\sigma_{x}$ used in Regime I for the observed variable, the
corresponding DA exercise is less challenging and will be presented first.

Results for Regime I. As is shown in Figure 2 for Regime I, the dynamics of
$x$ exhibits intermittent behavior with relatively quiescent episodes
punctuated by large excursion events. Due to the relatively small spectral
gap, the dynamics of $y$ also exhibits highly nonlinear oscillations sustained
by noise. In contrast, the dynamics of $z$ is mainly a damped oscillation
sustained by noise due to the relatively strong linear stabilizing effects,
and it is the variable that decays the fastest.

Figure 2: Solution of (16) in Regime I given by (20) for one arbitrarily fixed
realization of the noise, and the corresponding probability density functions
(PDFs) and autocorrelation functions (ACFs). This solution trajectory is taken
to be the true signal. See Section 4.2 for details about the numerical setup.
The ACFs and PDFs are estimated based on the solution trajectory over the time
window $[0,10^{4}]$, corresponding to $2\times 10^{7}$ data points for the
time step used. Figure 3: Panel (b): the filtered posterior mean of $y$ for
Regime I obtained from CG (green) and EnKBF (orange); the corresponding true
signal previously shown in Figure 2 is plotted in blue. Panel (c): analogue of
Panel (b) for $z$. Figure 4: Panel (a): The filtered posterior mean of $y$
obtained from BT (red) and the true signal of $y$ (blue). Panel (b): The true
signal of $yz$. The deterioration of the estimated $y$ using BT occurs over
time windows when the magnitude of $yz$ gets large.

The posterior mean states of $(y,z)$ for Regime I obtained by CG and EnKBF are
fairly close to each other as is shown in Figure 3. The conditional covariance
matrices of $(y,z)$ estimated by these two methods are also close to each
other for this regime (not shown). In contrast, for BT, there are prolonged
time windows over which the posterior mean of $y$ deviates significantly from
the true signal as is shown in Panel (a) of Figure 4. An inspection of the
time series of $yz$ shown in Panel (b) of Figure 4 reveals that such deviation
typically occurs when the value of $yz$ is large, which is expected, since the
omission of the term $\alpha yz$ in (16a) is the only difference between the
BT system and the full system. The posterior mean of $z$ obtained by BT is
similar to those obtained by CG and EnKBF shown in Panel (c) of Figure 3,
which is thus not presented.

Results for Regime II. The dynamics of the true system (16) in Regime II
exhibits similar features as in Regime I shown in Figures 2, although the
amplitude of each variable is slightly reduced due to the smaller noise
intensity $\sigma_{x}$ used for this regime. The shape of both the PDFs and
ACFs of all the three variables are similar to those shown in Figures 2,
except that the decorrelation time of $x$ becoming comparable with that of $y$
in this regime. The analogue of Figures 2 for Regime II is thus omitted.

Figure 5: Panel (b): the filtered posterior mean of $y$ for Regime II obtained
from CG (green) and EnKBF (orange) with the corresponding true signal shown in
blue. Panel (c): the analogue of Panel (b) for $z$. Panel (d): PDFs of the
filtered posterior mean of $y$ from CG (green) and EnKBF (orange) compared
with that of the true signal. Panel (e): the analogue of Panel (d) for $z$.

For this regime, CG and EnKBF still provide comparable posterior mean state of
$z$ as shown in Panel (c) of Figure 5, although the PDF of the posterior mean
obtained by EnKBF approximates slightly better the PDF of the true signal of
$z$. However, CG is significantly more skillful in estimating the conditional
mean of $y$ (Figure 5, Panels (b) and (d)).

The deterioration of the skill from EnKBF in this regime is associated with a
false bimodal behavior appearing in the posterior mean of the $y$ variable
(see Panel (d) of Figure 5), whereas the true signal is unimodal, skewed
towards negative values. It has also been checked that the bimodality is
always there for EnKBF by further increasing the total number of ensemble
members $N$ or decreasing the numerical integration time step $\delta t$. Such
a pathological behavior is associated with the filter divergence [52, 62],
which often occurs for the ensemble-based filters when the noise in the
observational process is small and the observational process is highly
nonlinear. In fact, the small noise in the observational process $x$ makes the
filter trusts more towards the information provided by the observations.
However, the strong nonlinear and non-Gaussian features of $x$ make it very
difficult to accurately recover the states of $y$ and $z$ by inferring mainly
from the $x$ process. In contrast, CG tracks well the modulations of the true
signal, leading to a much better reproduction of the PDF of the true signal of
$y$; see again Panel (d) of Figure 5. It is worth pointing out that the
original system (16) can also exhibit bimodal dynamics in a broad range of
dynamical regimes, even though bimodality is not observed in the true signals
of $(x,y,z)$ for neither of the two regimes considered here. This bimodality
that can occur in the dynamics of (16) is induced by the additive noise that
drives the system to switch from the two locally stable steady states of the
corresponding deterministic system produced from a supercritical pitchfork
bifurcation, although when occurs, the bimodality is mainly visible in the $x$
variable.

Regarding BT, compared with the corresponding result shown in Figure 4 for
Regime I, its performance here is even worse and is thus not shown. In
particular, the posterior mean state of $y$ not only deviates significantly
from the true signal but also has spurious fast oscillations presenting
throughout the whole time window. Such fast oscillations also appear in the
filtered posterior mean of $z$, although to a lesser extent.

For this regime, we also compared the skills of ensemble forecast with the
forecast model simply taken to be the true SDE system (16), and the initial
conditions (IC) of $(y,z)$ drawn from multivariate Gaussian distributions in
which the mean and covariance are estimated respectively from BT, CG and EnKBF
described above; see Section 4.2 for further details. In addition, we compute
the results when the forecast is initialized with the true signal for all the
three variables, which serves as the reference of the theoretic forecast
limit.

Figure 6: Panels (a)-(c): RMSE of the forecast skills for Regime II when the
ICs are chosen either to be the perfect values from the true signals, or are
drawn from multivariate Gaussian distributions in which the mean and
covariance are estimated respectively from BT, CG and EnKBF; the RMSE are
normalized by dividing by the standard deviation of the respective true
signal. Panels (d)-(f): the correlation coefficients of the forecast skills.
The vertical dashed line corresponds to lead time $\tau=0.4$ for which the
corresponding forecasted ensemble mean time series are shown in Figure 7.

In Figure 6, we presented the normalized root-mean-square error (RMSE;
normalized by the standard deviation of the true signal) and correlation
coefficients of the forecasts for all the methods used. As is expected, the
better skills of CG at the DA stage carries over to the ensemble forecast as
well. BT performs the worst due to large spurious oscillations appearing in
its DA stage. While the RMSE is a convenient way of ranking the performance,
to provide a better visualization of the skills, we also show in Figure 7 the
forecasted ensemble mean trajectories at a given lead time, chosen here to be
$\tau=0.4$ time unit, which corresponds roughly to one half of the
decorrelation time of the $y$-variable. The results in Figure 7 show that the
forecast based on initial conditions provided by CG (middle row) actually
performs fairly well for all the three variables compared with those when
perfect initial conditions are used (top row), although the uncertainty in the
$y$ variable is slightly higher for CG. The results for EnKBF are correlated
with its performance at the DA stage, with significant error in the $x$ and
$y$ variables over the time windows when the filtered posterior mean of $y$
deviates from the true signal as was previously shown in the middle panel of
Figure 5. For BT, large spurious oscillations appear in the forecasted
ensemble mean time series for all the variables, especially for $x$ and $y$.
Such oscillations are inherited from those appearing in the assimilated $y$
variables, which propagate to the other two variables due to nonlinear
interactions.

Figure 7: The forecasted ensemble mean time series for Regime II at lead time
$\tau=0.4$ time unit, when the ICs are from either the true signal (top row),
or assimilated from the CG method (middle row), or the EnKBF (bottom row). The
gray area on each panel marks the spread between 5 percentile and 95
percentile of the corresponding ensemble forecast. The true signals are
plotted in black. For BT, large spurious oscillations appear in the forecasted
ensemble mean time series for all the variables; and the corresponding results
are thus not shown here.

##  5 Parameter Estimation

Parameter estimation is an important topic and a necessary precursor for
effective state estimation, data assimilation, and prediction. Maximum
likelihood estimation (MLE) and maximum a posteriori (MAP) estimation are
often adopted to infer the parameter values if the observed time series for
all the state variables are accessible [96, 67, 30]. However, only partial
observations are available in many complex nonlinear systems. In such a
situation, data augmentation is widely used to simultaneously estimate the
model parameter and recover the unobserved state variables [114]. In
particular, data augmentation has been extensively incorporated into the
Markov Chain Monte Carlo (MCMC) algorithms for improving the Bayesian
inference [121, 51, 39, 100]. In contrast to targeting the global optimal
solution based on the MCMC, many other parameter estimation algorithms seek
locally optimal solutions. One of the widely used local optimal parameter
estimation approaches is the expectation-maximization (EM) algorithm [48, 47,
35, 91]. The EM algorithm is an iteration method that aims to find the
parameter values that maximize the likelihood function and compute certain
statistical expectations of the unobserved state variables in an alternating
fashion. Note that, due to the local optimality property, the EM algorithm
often requires fewer iterations than the MCMC algorithm. Unfortunately, both
methods are computationally expensive for general complex nonlinear systems
since neither the data augmentation in MCMC nor the computation of the
statistical expectation in the EM can be easily obtained.

Unlike the general complex nonlinear systems, the closed analytic formulae
offered by the CGNS facilitate the acceleration of the computational
efficiency in parameter estimation. In particular, the conditional sampling
formula (7) has the potential to allow a rapid data augmentation in the MCMC
method while the analytic state estimation formula (6) offers an exact and
accurate way to compute the statistical expectation, which is an essential
component in the EM algorithm. Note that appropriately incorporating the CGNS
into the MCMC may require several additional manipulations, which deserves a
separate topic to study. Therefore, the focus below is on applying the CGNS to
accelerate the parameter estimation utilizing the EM algorithm, where detailed
mathematical justifications are more accessible. The CGNS in this entire
procedure acts as a preconditioner to seek a suitable approximation of the
statistical expectation of the unobserved state variables given the
observational time series via the EM algorithm. With such a statistical
expectation available from the CGNS, the original complex nonlinear system is
only essential in computing the maximum likelihood solution at the last
iteration step, leading thus to a significant speedup of the computational
time required.

###  5.1 Accelerating the EM algorithm with a CGNS preconditioner

As before, denote by $(\mathbf{X},\mathbf{Y})$ the state variables of a
complex nonlinear system of the form (11), where only a time series of
$\mathbf{X}$ is observed while there is no direct observations for the state
variable $\mathbf{Y}$. Let $\bm{\theta}$ be a vector consisting of all the
model parameters. Denote by
$\widehat{\mathbf{X}}=\\{\mathbf{X}^{0},\ldots,\mathbf{X}^{j},\ldots,\mathbf{X}^{J}\\}$
and
$\widehat{\mathbf{Y}}=\\{\mathbf{Y}^{0},\ldots,\mathbf{Y}^{j},\ldots,\mathbf{Y}^{J}\\}$
a discrete approximation of the continuous time series of $\mathbf{X}$ and
$\mathbf{Y}$, respectively, within the time interval $t\in[0,T]$, where
$T=J\Delta{t}$, $\mathbf{X}^{j}=\mathbf{X}(t_{j})$ and
$\mathbf{Y}^{j}=\mathbf{Y}(t_{j})$, with $t_{j}=j\Delta t$.

The goal is to seek an optimal estimation of the unknown parameters
$\bm{\theta}$ by maximizing the log-likelihood function. Since only the time
series of $\mathbf{X}$ is observed (denoted by the discrete sequence
$\widehat{\mathbf{X}}$ above), the log-likelihood estimate is obtained by
averaging over the state variable $\mathbf{Y}$ at the corresponding instants,

$\mathcal{L}(\bm{\theta})=\log
q(\widehat{\mathbf{X}}|\bm{\theta})=\log\int_{\widehat{\mathbf{Y}}}p(\widehat{\mathbf{X}},\widehat{\mathbf{Y}}|\bm{\theta}){\,\rm
d}\widehat{\mathbf{Y}}.$ (25)

Due to the unknown state variable $\mathbf{Y}$, there is in general no simple
formula for a direct calculation of the log-likelihood. To find the optimal
parameter $\bm{\theta}$ that maximizes $\mathcal{L}$, the EM iteration
algorithm [18, 35, 64] operates on the following function instead which is
obtained as a lower bound of $\mathcal{L}(\bm{\theta})$ by applying Jensen’s
inequality (cf. [18, Section 3.1])

$\int_{\widehat{\mathbf{Y}}}Q(\widehat{\mathbf{Y}})\log
p(\widehat{\mathbf{Y}},\widehat{\mathbf{X}}|\bm{\theta}){\,\rm
d}\widehat{\mathbf{Y}}-\int_{\widehat{\mathbf{Y}}}Q(\widehat{\mathbf{Y}})\log
Q(\widehat{\mathbf{Y}}){\,\rm d}\widehat{\mathbf{Y}},$ (26)

where $Q(\widehat{\mathbf{Y}})$ is a distribution over the unobserved variable
$\widehat{\mathbf{Y}}$. It alternates between performing an expectation (E)
step to update $Q(\widehat{\mathbf{Y}})$ and a maximization (M) step to update
$\bm{\theta}$ until the estimated parameter $\bm{\theta}$ converges. Denote by
$\bm{\theta}_{k}$ the updated parameters after the $k$-th iteration. The EM
algorithm at step $k+1$ is the following:

  * E-Step.

Computing the conditional distribution
$p(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}_{k})$ using the
previously estimated parameters $\bm{\theta}_{k}$. In fact, the maximization
in the E-Step is reached when $Q(\widehat{\mathbf{Y}})$ is exactly the
conditional distribution of $\widehat{\mathbf{Y}}$ corresponding to the
smoother estimates.

  * M-Step.

Updating the parameters $\bm{\theta}_{k+1}$ by maximizing the cost function
$\mathcal{Q}$ defined by

$\mathcal{Q}(\bm{\theta};\bm{\theta}_{k})=\int_{\widehat{\mathbf{Y}}}p(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}_{k})\log
p(\widehat{\mathbf{Y}},\widehat{\mathbf{X}}|\bm{\theta}){\,\rm
d}\widehat{\mathbf{Y}}.$ (27)

That is,

$\bm{\theta}_{k+1}=\arg\max_{\bm{\theta}}\mathcal{Q}(\bm{\theta};\bm{\theta}_{k}).$
(28)

Note that in (27),
$p(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}_{k})$ is treated as a
known distribution that is solved in the E-step. Therefore,
$\mathcal{Q}(\bm{\theta};\bm{\theta}_{k})$ is a function of $\bm{\theta}$
only. The distribution
$p(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}_{k})$ can be regarded
as the weight function for computing the average (i.e., the integration) of
$\log p(\widehat{\mathbf{Y}},\widehat{\mathbf{X}}|\bm{\theta})$.

In many situations, the M-Step usually involves solving a quadratic
optimization problem, the analytic formula of which is available. However, for
general nonlinear systems, the conditional distribution
$p(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}_{k})$ in the E-Step
is extremely difficult to solve. Note that such a conditional distribution is
precisely the smoother estimate of the complex nonlinear system. Particle
methods can be applied. Yet, repeatedly using these particle methods through
the iteration procedure can be computationally expensive, and careful tuning
is required, especially for systems with large dimensions.

To overcome this most significant barrier of computing the conditional
distribution $p(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}_{k})$ in
the above EM algorithm, we exploit a suitable CGNS model as a preconditioner
to accelerate this calculation in the E-Step. To that end, denote by
$\bm{\theta}^{M}$ the collection of the parameters in the chosen CGNS model.
In stark contrast to
$p(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}_{k})$, we can now use
the closed analytic formulae in (6) to efficiently compute the conditional
distribution
$p^{M}(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}^{M}_{k})$ for the
CGNS at each EM iteration.111We should have used
$p^{M}(\widehat{\mathbf{Y}}^{M}|\widehat{\mathbf{X}},\bm{\theta}^{M}_{k})$ to
denote the conditional distribution for the CGNS given the observation
$\widehat{\mathbf{X}}$, where $\widehat{\mathbf{Y}}^{M}$ is the analogue of
$\widehat{\mathbf{Y}}$ for the CGNS. We used
$p^{M}(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}^{M}_{k})$ instead
to avoid excessive notations.  When $\bm{\theta}^{M}_{k}$ is converged after a
sufficient number of EM iterations based on the CGNS, to obtain an
approximation of the optimal $\bm{\theta}$ for the original system, we simply
perform another EM cycle, but this time using $\log
p(\widehat{\mathbf{Y}},\widehat{\mathbf{X}}|\bm{\theta})$ instead of $\log
p^{M}(\widehat{\mathbf{Y}},\widehat{\mathbf{X}}|\bm{\theta}^{M})$ in the
M-Step. The above procedure is summarized in Algorithm 5.1; and we refer to
Appendix A for further technical details. Throughout this section, a quantity
with the superscript $M$ indicates that it is a quantity related to the
approximate CGNS model instead of the original model (see also footnote 1 with
the exception of the notation $\widehat{\mathbf{Y}}$). Note also that
$\bm{\theta}^{M}$ can be different from $\bm{\theta}$, since the CGNS can
involve new parameters not appearing in the original system and vice versa.

[h]  EM algorithm for nonlinear systems with CGNS as a preconditioner Start
with a given realization of the (partial) observations $\widehat{\mathbf{X}}$
for (11);

2 Propose an approximate CGNS model of the form (2);

3 Assign an initial guess for the CGNS model’s parameters
$\bm{\theta}^{M}_{0}$;

4 for  _$k=1:K$_ do

5 E-step: compute the conditional distribution
$p^{M}(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}^{M}_{k-1})$ using
the previously estimated parameters $\bm{\theta}^{M}_{k-1}$;

6 M-step: update the parameters $\bm{\theta}^{M}_{k}$ with
$\bm{\theta}^{M}_{k}=\arg\max_{\bm{\theta}^{M}}\widetilde{\mathcal{Q}}(\bm{\theta}^{M};\bm{\theta}^{M}_{k-1})$,
where
$\widetilde{\mathcal{Q}}(\bm{\theta}^{M};\bm{\theta}^{M}_{k-1})=\int_{\widehat{\mathbf{Y}}}p^{M}(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}^{M}_{k-1})\log
p^{M}(\widehat{\mathbf{Y}},\widehat{\mathbf{X}}|\bm{\theta}^{M}){\,\rm
d}\widehat{\mathbf{Y}}$.

7E-step: compute the conditional distribution
$p^{M}(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}^{M}_{K})$;

M-step: compute the optimal parameters $\bm{\theta}_{\text{opt}}$ for (11)
with
$\bm{\theta}_{\text{opt}}=\arg\max_{\bm{\theta}}\mathcal{Q}(\bm{\theta};\bm{\theta}^{M}_{K})$,
where
$\mathcal{Q}(\bm{\theta};\bm{\theta}^{M}_{K})=\int_{\widehat{\mathbf{Y}}}p^{M}(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}^{M}_{K})\log
p(\widehat{\mathbf{Y}},\widehat{\mathbf{X}}|\bm{\theta}){\,\rm
d}\widehat{\mathbf{Y}}$.

Note that $\bm{\theta}_{\text{opt}}$ obtained from Algorithm 5.1 should be
viewed as an approximation of the true optimal parameters for the original
nonlinear system. Its quality relies obviously on the quality of the CGNS in
approximating the true dynamics. It is also worth mentioning that additional
(physics) constraints can be naturally incorporated into the CGNS for
parameter estimation, while still preserving the closed analytic formulae in
the corresponding EM algorithm; and under certain conditions, a block
decomposition of the conditional covariance can be devised to further reduce
the computational efforts when high dimensional systems are considered [18].
We illustrate now the approach on the classical two-layer Lorenz 1996 model.

###  5.2 A multiscale turbulent test model

In this subsection, a two-layer inhomogeneous Lorenz model is utilized to
demonstrate that a suitable CGNS can serve as both a preconditioner and a
surrogate model in parameter estimation to speed up the computation. First, we
show that a simple approximate model that belongs to CGNS can accelerate the
EM algorithm as a preconditioner by following Algorithm 5.1. Then, we show
that this approximate model itself can be used as a surrogate model for
prediction once the involved parameters, $\bm{\theta}^{M}$, are optimized
through the EM cycles in lines 4-6 of Algorithm 5.1.

####  5.2.1 The perfect model

The two-layer Lorenz 96 (L96) model [74] is a conceptual representation of
geophysical turbulence that is commonly used as a testbed for various
stochastic parameterization and dimension reduction techniques [29, 42, 2, 31,
119, 85, 54, 53, 84]. The model mimics a coarse discretization of atmospheric
flow on a latitude circle. It supports complex wave-like and chaotic behavior,
and the two-layer structure schematically depicts the interactions between
small-scale fluctuations and large-scale motions. The stochastic version of
the model subject to additive noise forcing reads

$\displaystyle\frac{du_{i}}{dt}$
$\displaystyle=\left(-u_{i-1}\left(u_{i-2}-u_{i+1}\right)-u_{i}+f-\frac{hc_{i}}{J}\sum^{J}_{j=1}v_{i,j}\right)+\sigma_{u_{i}}\dot{W}_{u_{i}},\quad
i=1,\dots,I,$ (29a) $\displaystyle\frac{dv_{i,j}}{dt}$
$\displaystyle=\left(-bc_{i}v_{i,j+1}\left(v_{i,j+2}-v_{i,j-1}\right)-c_{i}v_{i,j}+\frac{hc_{i}}{J}u_{i}\right)+\sigma_{v_{i,j}}\dot{W}_{v_{i,j}},\quad
j=1,\dots,J,$ (29b)

where $I$ denotes the total number of large-scale variables, $J$ the number of
small-scale variables corresponding to each large-scale variable, $f$, $h$,
$c_{i}$, $b$, $\sigma_{u_{i}}$ and $\sigma_{v_{i,j}}$ are given scalar
parameters while $\dot{W}_{u_{i}}$ and $\dot{W}_{v_{i,j}}$ are white noise.
The large-scale variables $u_{i}$ are periodic in $i$ with
$u_{i+I}=u_{i-I}=u_{i}$. The corresponding small-scale variables $v_{i,j}$ are
periodic in $i$ with $v_{i+I,j}=v_{i-I,j}=v_{i,j}$ and satisfy the following
cyclic conditions in $j$: $v_{i,j+J}=v_{i+1,j}$, and $v_{i,j-J}=v_{i-1,j}$.

The model discussed here uses variables $u_{i}$ to describe large-scale or
slow movements which are resolved; small scales or rapid fluctuations
represented by $v_{i,j}$ are often unresolved ones. The coupling of fast and
slow variables is regulated by the parameter $h$. The parameter $f$ controls
the magnitude of external large-scale forcing, while $b$ determines the
amplitude of nonlinear interactions between the fast variables. The parameter
$c_{i}$ specifies how quickly the fast variables are damped in comparison to
the slow variables. We take $I=40$, corresponding to a discretization of the
latitude circle into a total of 40 sites/sectors, and choose $J=4$ small-scale
variables associated with each $u_{i}$. The constant forcing is set to be
$f=4$, which makes the system chaotic for the parameter regime chosen here.
The parameters $h$, $c_{i}$, and $b$ are chosen in such a way that the small-
scale variables have a comparatively significant impact on the large-scale
ones. In other words, the perfect model only has a weak scale separation. The
reason that we consider such a weak scale separation is that it better mimics
the real atmosphere with chaotic/turbulent behavior. The parameter $c_{i}$
varies across the spatial sites, which aims to mimic the fact that the
coupling across the variables above the ocean is weaker than that above the
land since the latter usually have stronger friction or dissipation. In this
sense, the model is inhomogeneous. Finally, additional stochastic noise is
added to the system, representing the contribution of the variables that are
not explicitly modeled. The noise also interacts with the deterministic part
via nonlinear terms, introducing additional complexity that mimics nature. To
summarize, the parameters used in the perfect model (29) are as follows,

$\displaystyle I=40,\quad J=4,\quad h=2,\quad c_{i}=2+0.7\cos(2\pi i/I),\quad
b=2,\quad f=4,$ (30)
$\displaystyle\sigma_{u_{i}}=\sigma_{u}=0.2,\quad\sigma_{v_{i,j}}=\sigma_{v}=1.$

####  5.2.2 The approximate model

Since in general the perfect model is not always fully known, or it is too
complicated to be used in practice, it is essential to develop a simple and
computationally tractable approximate model, which is nevertheless able to
capture the key nonlinear feedback from the unobserved variables ($v_{i,j}$
here) to the observed variables ($u_{i}$ here). As was discussed in Section 3,
stochastic parameterization is widely used in describing chaotic signals [21],
which replaces the nonlinear eddy terms by quasilinear stochastic processes on
formally infinite embedded domains where the stochastic processes are Gaussian
conditional to the large scale mean flow. In addition, physics-constraint is
adopted in designing approximate models, which also includes the effects from
the large-scale $u_{i}$ to the small-scale variables $v_{i,j}$. Therefore,
such approximate models can potentially be used as surrogate models of the
perfect model.

The approximate model that we introduce for (29) is as follows

$\displaystyle\frac{du_{i}}{dt}$
$\displaystyle=\left(-u_{i-1}\left(u_{i-2}-u_{i+1}\right)-u_{i}+\hat{f}_{i}-\hat{a}_{i}\sum^{J}_{j=1}v_{i,j}\right)+\hat{\sigma}_{u_{i}}\dot{W}_{u_{i}},\quad
i=1,\dots,I,$ (31a) $\displaystyle\frac{dv_{i,j}}{dt}$
$\displaystyle=-\hat{d}_{i,j}v_{i,j}+\hat{v}_{i,j}+\hat{c}_{i}u_{i}+\hat{\sigma}_{v_{i,j}}\dot{W}_{v_{i,j}},\quad
j=1,\dots,J,~{}$ (31b)

where $\hat{f}_{i}$, $\hat{a}_{i}$, $\hat{\sigma}_{u_{i}}$, $\hat{d}_{i,j}$,
$\hat{v}_{i,j}$, $\hat{c}_{i}$, and $\hat{\sigma}_{v_{i,j}}$ are unknown
constants. Compared with the original system (29), the main simplification
here is in the small-scale equations, where we have replaced the small-scale
nonlinear interaction and linear damping terms
$-bc_{i}v_{i,j+1}(v_{i,j+2}-v_{i,j-1})-c_{i}v_{i,j}$ in (29b) by a simple
linear term $-\hat{d}_{i,j}v_{i,j}+\hat{v}_{i,j}$. The equations for $u_{i}$
are essentially the same as before, although $\hat{f}_{i}$ is allowed now to
vary from site to site. In the numerical experiments below, we will enforce
$\hat{a}_{i}=\hat{c}_{i}$ for the consistency of the dynamical property with
the original system (29).

One desirable feature of this approximate model is that the direct coupling of
the state variables only involves $u_{i}$ and the corresponding $v_{i,j}$ for
each fixed $i$. This is different from the original system (29) where $u_{i}$
can have direct interactions with $v_{i+1,1}$ and $v_{i-1,J}$ due to the
cyclic boundary conditions of the small scale $v_{i,j}$ over $j$. Such a
property allows to use a block decomposition of the covariance matrix of the
smoother estimate during both E-step and M-step [18]. The entire state space
for all the variables $\\{u_{i},v_{i,j}\;|\;i=1,\dots,I,j=1,\dots J\\}$ can be
decomposed into $I$ subspaces, where each subspace can be dealt with in
parallel.

####  5.2.3 Setup of the numerical simulations

The true signal is obtained by integrating the inhomogeneous L96 model (29)
using the Euler-Maruyama scheme with the parameters given by (30), a uniform
time step size $\delta t=2\times 10^{-3}$, and zero initial condition. The
same time step size and initial condition are adopted for the approximate
model (31) as well as the identified model, where the latter takes the same
form as (29) but with estimated parameters. The true signals of $u_{i}$ with
100 time units are used as the observations for parameter estimation while
longer data with 1000 time units are used to compute their statistics. This
latter length is also adopted when computing the statistics of the variables
in the approximate and the identified models. The total number of the EM loops
with CGNS is fixed to be 200.

####  5.2.4 CGNS as a preconditioner for identifying parameters in the perfect
model

We first discuss the results of applying Algorithm 5.1 to estimate the
parameters, $\bm{\theta}^{M}$, in the approximate model (31). Figure 8 shows
that the EM algorithm with the approximate model (31) provides a very accurate
approximation of the parameters corresponding to the Gaussian fit of $v_{i,j}$
in the true signal. The quantities shown in Figure 8 are the trace plots at
the site $i=2$ corresponding to the variables $u_{2}$ and $v_{2,j}$ (in the
first and the third rows) and the final estimation of the parameters in (31)
(in the second and the fourth rows). The trace plots at other spatial
locations have similar behavior; the parameters involved in the
$u_{i}$-equations all converge quickly (within 10 iteration steps), and those
involved in the $v_{i,j}$-equations converge at a relatively slower speed, but
all stabilized after about 100 iterations. The black curves are shown as a
reference, where $\hat{f}_{i}$, $\hat{a}_{i}$ ($=hc_{i}/J$),
$\hat{\sigma}_{u_{i}}$, and $\hat{c}_{i}$ are taken to be those in the perfect
model and $\hat{d}_{i,j}=\hat{d}_{i}$, $\hat{v}_{i,j}=\hat{v}_{i}$, and
$\hat{\sigma}_{v_{i,j}}=\hat{\sigma}_{v_{i}}$ are calibrated by the true
statistics, i.e., the mean, the variance, and the decorrelation time, in the
perfect model of $v_{i,j}$ averaged over $j$. More precisely, denoting by
$\mu_{i}$, $r_{i}$, and $\tau_{i}$ these averaged mean, variance, and
decorrelation time in the perfect model of $v_{i,j}$. the parameters
$\hat{v}_{i}$, $\hat{d}_{i}$, and $\hat{\sigma}_{v_{i}}$ are then obtained via

$\tau_{i}=\frac{1}{\hat{d}_{i}},\qquad
r_{i}=\frac{\hat{\sigma}_{v_{i}}^{2}}{2\hat{d}_{i}},\qquad\mu_{i}=\frac{\hat{v}_{i}+hc_{i}/J\bar{u}_{i}}{\hat{d}_{i}},$
(32)

where $\bar{u}_{i}$ is the mean value of $u_{i}$.

Figure 8: Learning the parameters of the approximate L96 model (31). The
panels (a)–(d) and (i)–(k) in the 1st and 3rd rows show the trace plots of the
estimated parameters of the EM algorithm 5.1 at site $i=2$. The panels (e)–(h)
and (l)–(n) in the 2nd and 4th rows show the final estimated parameters
$\bm{\theta}^{M}_{K}$ (blue) after the $K$-th step of the EM algorithm with
$K$ taken here to be 200; also shown are the initial guesses of the parameters
(red) as well as the reference parameter values (back) chosen as follows:
$\hat{f}_{i}$, $\hat{a}_{i}$ ($=hc_{i}/J$), $\hat{\sigma}_{u_{i}}$, and
$\hat{c}_{i}$ are set to be the corresponding true values used in the perfect
model, while $\hat{d}_{i,j}=\hat{d}_{i}$, $\hat{v}_{i,j}=\hat{v}_{i}$, and
$\hat{\sigma}_{v_{i,j}}=\hat{\sigma}_{v_{i}}$ are calibrated by the true
statistics according to (32).

The conditional distribution of the approximate model (31) with the estimated
parameters is shown in Figure 9. Panel (a) and (c) show the true signal of two
large-scale modes $u_{10}$ and $u_{20}$, and Panel (c) and (d) show the true
signal, and smoother estimate of two hidden modes. The smoother mean is given
by the dashed-black curves, and one, two, and three standard derivations of
the uncertainty are shown by the light, moderate, and dark shading areas. The
shading areas cover most of the true signal, which indicates appropriate
amount of uncertainties are obtained by combining the true observations of the
large scale variables $u_{i}$ and the optimized approximate model. In fact,
characterizing appropriate amount of uncertainty plays an important role in
the EM algorithm when the hidden process contains large uncertainty. If the
uncertainty is totally ignored in computing the optimization in the M-Step,
for example replacing
$p(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}^{M}_{k})$ by its mean
state, then the estimated parameters can be very biased. In fact, the solution
of the estimation even blows up in the test model used here.

Figure 10 shows the estimated parameters $\bm{\theta}_{\text{opt}}$ for the
perfect model, where the CGNS (31) is utilized as a preconditioner following
Algorithm 5.1. Due to the intrinsic model error of the approximate CGNS model
where the hidden variables are fully decomposed, in the sense that the
correlation between the small scales corresponding to different large scales
are omitted, the estimated value for the $bc_{i}$ term is zero. However, other
parameters, for example $f$ and $c_{i}$, are adjusted accordingly. It is shown
in Figure 11 and Figure 12 that the identified full model with the estimated
parameters $\bm{\theta}_{\text{opt}}$ produces dynamics that resembles the
truth to a remarkable extent. Indeed, Figure 11 shows that the hovmoller plot
of the large scale variables from the identified model (Panel (b)) is almost
the same as that from the true signal (Panel (a)) over the given time window.
Time series comparison as well as PDF and ACF comparisons are also shown in
the first four row of Figure 12. Both the dynamical properties and the
statistics are recovered with high accuracy. Of course, the time series from
the identified model should not be expected to follow the true time series at
all time instants since the original model is placed in a regime with chaotic
dynamics.

Figure 9: Smoother estimate of the approximate model in the $K$-th iteration
in Algorithm 5.1, with $K=200$. The blue curves: true signals; the black
dashed curves: the smoother mean time series of the hidden variable; the
light, moderate, and dark shading areas show the one, two, and three standard
derivations (STDs) of the uncertainty in the smoother estimate. Panel (a):
true signal of $u_{10}$; Panel (b): true signal of $v_{10,1}$ with one, two,
and three, standard derivations of the uncertainty. Panel (c)–(d): Same as
Panel (a)–(b) but for $u_{20}$ and $v_{20,1}$. Figure 10: The estimated
parameter $\bm{\theta}_{\text{opt}}$ using Algorithm 5.1 for the perfect L96
model (29). Here, $\bm{\theta}_{\text{opt}}$ is a vector consisting of the
parameters $f$, $hc_{i}/J$, $\sigma_{u}$, $bc_{i}$, $c_{i}$, and $\sigma_{v}$
in (29). The true values given by (30) are marked by the dashed black curves,
and the estimated parameters are shown by the blue curves.

####  5.2.5 CGNS as a surrogate model

Finally, we mention that the approximate model (31) with the estimated
parameters itself can be exploited as a surrogate model of the perfect system,
which can be applied for ensemble forecast and other tasks. Indeed, (31) with
the optimally estimated parameters $\bm{\theta}^{M}_{K}$ from Algorithm 5.1
recovers the dynamical properties of the true model to an extent that is
almost the same as the full model with the identified parameters
$\bm{\theta}_{\text{opt}}$ as shown in Figure 11 and the last four rows in
Figure 12.

Figure 11: Hovmoller diagram of the large scale variables $u_{i}$ from
different models. Panel (a): perfect model (29) with parameters (30); Panel
(b): perfect model (29) with estimated parameters $\bm{\theta}_{\text{opt}}$
after line 8 using Algorithm 5.1; Panel (c): approximate model (31) with
parameters $\bm{\theta}^{M}_{K}$ with $K=200$. Figure 12: Comparison of time
series and statistics obtained from the full two-layer L96 model (29) with the
true parameters (30) (blue), the identified model (29) with the estimated
parameters $\bm{\theta}_{\text{opt}}$ from Algorithm 5.1 (green), and the
approximate CGNS model (31) with the estimated parameters
$\bm{\theta}^{M}_{K}$ (orange) after the $K$-th iteration in Algorithm 5.1
where $K=200$. Panel (a)–(c): trajectory, PDF, and ACF, respectively.

##  6 Predicting the Statistical Response

Yet another important topic in studying complex nonlinear systems is to
predict the model response to the perturbation of the external forcing.
Developing efficient and accurate approaches to study such a key issue
facilitates understanding model sensitivity, regime-switching behavior, and
nonlinear interactions across different scales. Resolving this problem using
advanced mathematical tools also has significant practical implications, such
as coping with the climate change scenario. Due to various uncertainties from
the internal instability and external forcing, a probabilistic description is
more suitable for characterizing the complex turbulent systems (1).

However, there exist several challenges in predicting the statistical response
of complex nonlinear systems. First, solving the high-dimensional Fokker-
Planck equation is the prerequisite of obtaining the model statistics, which,
however, often suffers from the curse of dimensionality. Second, since
simulating the perfect system is not always computationally feasible in
practice, the predicted model response may become inaccurate when an
approximate model is utilized. It is then important to take advantage of a
suitable combination of partial observations with the approximate model to
mitigate error in calculating the model statistics. Third, for general
nonlinear systems, computing the statistical response in terms of the
perturbations with different strengths and categories requires repeatedly
solving the Fokker-Planck equation. As a consequence, even with a relatively
fast solver of the Fokker-Planck equation, the total computational cost can
still remain significant.

The advantage of the exact and statistically accurate solver of the
equilibrium PDF (9) makes the CGNS a natural framework for the development of
approximate models for efficiently computing the model statistics and the
associated response. One important feature in finding the PDF of the CGNS
based on the formula (9) is that the available partially observed time series
$\mathbf{X}$ is used to compute the conditional distribution of $\mathbf{Y}$
given $\mathbf{X}$. As a consequence, the model error in the PDF associated
with the direct simulation of the approximate model is mitigated with the help
of observations in computing the PDF based on (9). In other words, the
resulting PDF from (9) is in general closer to the truth than that computed
purely based on the approximate model without taking into account any input
information from observations.

What remains is the third challenge mentioned above. To overcome such a
difficulty, the linear statistical response is utilized as an approximate
method for computing the exact statistical response. The linear response only
requires a linearization of the statistical equation while there is no
linearization involved in the original nonlinear dynamics. Therefore, the
nonlinear features of the underlying dynamics is preserved. In addition, the
linear response can be computed utilizing the fluctuation-dissipation theorem
(FDT) [77], which involves only a single PDF for computing different linear
responses. Note that such a PDF is the equilibrium distribution of the
unperturbed state, which can be efficiently solved using the formula (10).

###  6.1 Computing the linear statistical response via the fluctuation-
dissipation theorem (FDT)

Consider the general complex nonlinear systems in (1). Defining
$\mathbf{G}(\mathbf{u},t)=\left(L+D\right)\mathbf{u}+B\left(\mathbf{u},\mathbf{u}\right)+\mathbf{F}(t)$,
the model can be written in a concise form

$\frac{{\,\rm d}\mathbf{u}}{{\,\rm
d}t}=\mathbf{G}(\mathbf{u},t)+\bm{\sigma}(\mathbf{u},t)\dot{\mathbf{W}}.$ (33)

The equilibrium statistics of some functional $A(\mathbf{u})$ associated with
(33) is formulated as

$\langle A(\mathbf{u})\rangle=\int
A(\mathbf{u})p_{\textrm{eq}}(\mathbf{u}){\,\rm d}\mathbf{u},$ (34)

where $p_{\textrm{eq}}(\mathbf{u})$ is the equilibrium PDF of $\mathbf{u}$ in
(33).

Now consider the dynamics in (33) by a small time separable external forcing
perturbation $\delta\mathbf{w}(\mathbf{u})f(t)$, where $\delta$ is a small
scalar $\mathbf{w}$ is a general nonlinear function of $\mathbf{u}$. The
perturbed system reads

$\frac{{\,\rm d}\mathbf{u}}{{\,\rm
d}t}=\mathbf{G}(\mathbf{u},t)+\delta\mathbf{w}(\mathbf{u})f(t)+\bm{\sigma}(\mathbf{u},t)\dot{\mathbf{W}}.$
(35)

The FDT states that if $\delta$ is small enough, then the leading-order
correction to the statistics in (34) is given by

$\delta\langle
A(\mathbf{u})\rangle(t)=\delta\int_{0}^{t}\mathbf{R}(t-s)f(s){\,\rm d}s,$ (36)

where $\mathbf{R}(t)$ is the linear response operator, which is calculated
through correlation functions in the unperturbed dynamics:

$\mathbf{R}(t)=\langle
A(\mathbf{u}(t))\mathcal{B}(\mathbf{u}(0))\rangle,\qquad\mathcal{B}(\mathbf{u})=-\frac{\mbox{div}_{\mathbf{u}}(\mathbf{w(\mathbf{u})}p_{\textrm{eq}}(\mathbf{u}))}{p_{\textrm{eq}}(\mathbf{u})}.$
(37)

See [77, Chapter 2] for a rigorous derivation of (36)–(37). In particular, the
above procedure of computing the linear response via the FDT does not require
the linearization of the underlying complex nonlinear systems. Therefore, the
features of the nonlinear dynamics are preserved. Note also that if the
functional $A(\mathbf{u})$ in (36) is given by $A(\mathbf{u})=\mathbf{u}$,
then the response computed is for the statistical mean. Likewise,
$A(\mathbf{u})=(\mathbf{u}-\bar{\mathbf{u}})^{2}$ is used for computing the
response in the variance.

###  6.2 Calculating the linear statistical response via the CGNS
preconditioner

According to (37), the calculation of the linear statistical response via the
FDT requires the information of

  1. 1.

the equilibrium PDF $p_{\textrm{eq}}(\mathbf{u})$,

  2. 2.

the time series $A(\mathbf{u})$, and

  3. 3.

the correct formulation of $\mathcal{B}(\mathbf{u})$.

Even if the perfect model is known, directly solving the high-dimensional
Fokker-Planck equation is often not computationally affordable. Therefore, a
suitable CGNS, serving as a preconditioner, is utilized to find a suitable
approximation of the non-Gaussian equilibrium PDF $p_{\textrm{eq}}$ in an
efficient way. Specifically, the explicit formula in (10) is utilized to
achieve this goal. Besides, the observations also need to be incorporated in
computing the equilibrium $p_{\textrm{eq}}$ (and later in recovering the
hidden components in $A(\mathbf{u})$) to reduce model error from the
approximate model free run. We denote one realization of the observational
variable by $\mathbf{X}^{\textrm{obs}}$, the posterior distribution (filter
and smoother) given $\mathbf{X}^{\textrm{obs}}$ by $p^{M|\textrm{obs}}$ where
the superscript $M$ indicates that an approximate model is used in computing
the filter distribution (4), the explicit formula (10) becomes

$p_{\textrm{eq}}^{M|\textrm{obs}}(\mathbf{X},\mathbf{Y})=\lim_{J\to\infty}\frac{1}{J}\sum_{j=1}^{J}\Big{(}K_{\mathbf{H}}(\mathbf{X}-\mathbf{X}^{\textrm{obs}}(t_{j}))p^{M|\textrm{obs}}(\mathbf{Y}|\mathbf{X}^{\textrm{obs}}(s\leq
t_{j}))\Big{)},$ (38)

where $p_{\textrm{eq}}^{M|\textrm{obs}}$ is an efficient and effective
approximation of the true equilibrium PDF $p_{\textrm{eq}}$.

Next, in the presence of partial observations, the conditional sampling
formula (7) is exploited to calculate the unobserved component of the time
series in $A(\mathbf{u})$. Note that the approximate model is used to compute
the filter distribution (4), the smoother distribution formulae (6) and the
conditional sampling formula (7). Note that the partial observations are
involved in computing both $p_{\textrm{eq}}(\mathbf{u})$ and $A(\mathbf{u})$,
which aims to mitigate the model error in the approximate model in both the
equilibrium PDF and time series of the unobserved components. In addition,
with $p_{\textrm{eq}}(\mathbf{u})$ and $A(\mathbf{u})$ obtained from the CGNS
preconditioner, the original nonlinear system structure is utilized to form
$\mathcal{B}(\mathbf{u})$ to compute the linear response $\mathbf{R}(t)$. The
entire procedure of the FDT via CGNS is given in Algorithm 1.

1 Start with a given realization of the observations
$\mathbf{X}^{\textrm{obs}}$;

2 Propose an approximate model that belongs to CGNS (2);

3 Compute the filter posterior distribution
$p^{M|\textrm{obs}}(\mathbf{Y}|\mathbf{X}^{\textrm{obs}}(s\leq t_{i}))$ via
(4) ;

4 Form the equilibrium $p_{\textrm{eq}}^{M|\textrm{obs}}$ via equation (38);

5 Compute the smoother posterior distribution
$p^{M|\textrm{obs}}(\mathbf{Y}(t)|\mathbf{X}^{\textrm{obs}}(s),s\in[0,T])$
from (6);

6 Sample one realization of the hidden time series via (7) that is used to
approximate the unobserved component of $A(\mathbf{u})$;

Compute the response operator $\mathbf{R}(t)$ via (37) and compute the linear
response via (36).

Algorithm 1 FDT with the CGNS preconditioner

###  6.3 A 4D stochastic climate model

This section utilizes a four-mode stochastic model with key features of
atmospheric low-frequency variability to show how to use CGNS as a
preconditioner incorporated with partial observations to calculate the linear
statistical response.

####  6.3.1 The perfect model

The stochastic climate model is designed in such a way that it involves many
of the major dynamical properties of comprehensive global circulation models
(GCMs) but with only four degree of freedom [82, 77, 88, 89]. The model reads
as follows

$\displaystyle\frac{{\,\rm d}x_{1}}{{\,\rm d}t}$
$\displaystyle=\left(-x_{2}(L_{12}+a_{1}x_{1}+a_{2}x_{2})-d_{1}x_{1}+F_{1}+L_{13}y_{1}+b_{123}x_{2}y_{1}\right)+\sigma_{1}\dot{W}_{x_{1}},$
(39a) $\displaystyle\frac{{\,\rm d}x_{2}}{{\,\rm d}t}$
$\displaystyle=\left(+x_{1}(L_{12}+a_{1}x_{1}+a_{2}x_{2})-d_{2}x_{2}+F_{2}+L_{24}y_{2}+b_{213}x_{1}y_{1}\right)+\sigma_{2}\dot{W}_{x_{2}},$
(39b) $\displaystyle\frac{{\,\rm d}y_{1}}{{\,\rm d}t}$
$\displaystyle=\left(-L_{13}x_{1}+b_{312}x_{1}x_{2}+F_{3}-\gamma_{1}y_{1}\right)+\sigma_{3}\dot{W}_{y_{1}},$
(39c) $\displaystyle\frac{{\,\rm d}y_{2}}{{\,\rm d}t}$
$\displaystyle=\left(-L_{24}x_{2}+F_{4}-\gamma_{2}y_{2}\right)+\sigma_{4}\dot{W}_{y_{2}},$
(39d)

where $b_{123}+b_{213}+b_{312}=0$. Consistent with many geophysical flow
models, the model has energy-conserving quadratic nonlinear terms, a linear
operator, and external forcing terms. The linear operator contains two parts:
one is a skew-symmetric component formally related to the Coriolis effect and
topographic Rossby wave propagation; the other is a negative definite
symmetric portion conceptually analogous to dissipative processes such as
surface drag and radiative damping. The coupling in different variables is
through both linear and nonlinear terms, where the nonlinear coupling through
$b_{ijk}$ produces multiplicative noise effects. In fact, the strategies
described in Section 3.1 are applied to $y_{1}$ and $y_{2}$ that introduce the
stochastic noise and damping terms. The variables $x_{1}$ and $x_{2}$ can be
regarded as the climate variables and $y_{1}$ and $y_{2}$ represents the
weather variables. The parameters used to generate the true dynamics are as
follows

$\displaystyle d_{1}=1,\quad
d_{2}=0.4,\quad\gamma_{1}=0.5,\quad\gamma_{2}=0.5,\quad L_{12}=1,\quad
L_{13}=0.5,\quad L_{24}=0.5,$ (40) $\displaystyle a_{1}=2,\quad a_{2}=1,\quad
b_{123}=1.5,\quad b_{213}=1.5,$
$\displaystyle\sigma_{1}=0.5,\quad\sigma_{2}=2,\quad\sigma_{3}=0.5,\quad\sigma_{4}=1,\quad
F_{1}=F_{2}=F_{3}=F_{4}=0.$

One realization of the true signal is shown in black color in Figure 13. Both
climate variable $x_{1}$ and weather variable $y_{1}$ have intermittent
behavior with non-Gaussian PDFs. Note that this stochastic model is CGNS with
$\mathbf{X}=(x_{1},x_{2})^{\top}$ and $\mathbf{Y}=(y_{1},y_{2})^{\top}$. We
use a CGNS as the perfect model such that the FDT based on the perfect model
can be computed in an accurate fashion, which can be served as a reference
solution.

####  6.3.2 The approximate model

In practice, running the entire perfect model is prohibitively costly. As a
result, simpler or reduced models are commonly utilized in computing the
responses. Linear stochastic models are widely used as approximate models for
the unresolved variables [34]. Therefore, the hidden processes are replaced by
two linear Gaussian equations, the parameters of which are calibrated by the
true equilibrium statistics, i.e., the mean, the variance, and the
decorrelation time. Besides, the parameters in the observed processes are
assumed to be the same as those in the perfect model (39). The approximate
model reads,

$\displaystyle\frac{{\,\rm d}x_{1}}{{\,\rm d}t}$
$\displaystyle=\left(-x_{2}(L_{12}+a_{1}x_{1}+a_{2}x_{2})-d_{1}x_{1}+F_{1}+L_{13}y_{1}+b_{123}x_{2}y_{1}\right)+\sigma_{1}\dot{W}_{x_{1}},$
(41a) $\displaystyle\frac{{\,\rm d}x_{2}}{{\,\rm d}t}$
$\displaystyle=\left(+x_{1}(L_{12}+a_{1}x_{1}+a_{2}x_{2})-d_{2}x_{2}+F_{2}+L_{24}y_{2}+b_{213}x_{1}y_{1}\right)+\sigma_{2}\dot{W}_{x_{2}},$
(41b) $\displaystyle\frac{{\,\rm d}y_{1}}{{\,\rm d}t}$
$\displaystyle=-\hat{\gamma}_{3}\left(y_{1}-\hat{y}_{1}\right)+\hat{\sigma}_{3}\dot{W}_{y_{1}},$
(41c) $\displaystyle\frac{{\,\rm d}y_{2}}{{\,\rm d}t}$
$\displaystyle=-\hat{\gamma}_{4}\left(y_{2}-\hat{y}_{2}\right)+\hat{\sigma}_{4}\dot{W}_{y_{1}},$
(41d)

which belongs to the CGNS. Note that despite the simplicity of utilizing
linear Gaussian models to approximate the hidden processes, one major issue in
(41) is that the physics constraint is no longer satisfied in (41). Therefore,
the model in (41) can contain large errors for a long-term simulation.

####  6.3.3 Calculating linear response via FDT

In this example, we aim at calculating the linear response to the perturbation
of the external forcing and linear interaction parameters. Specifically, the
following two perturbation cases are considered:

  * Case 1: Perturbing parameters of forcing in the observed processes, i.e., $F_{1}^{\delta}=F_{2}^{\delta}=0.3$, $\delta\mathbf{w}(\mathbf{u})f(t)=(0.3,0.3,0,0)^{\top}$;

  * Case 2: Perturbing parameters in linear interaction terms, i.e., $L_{13}^{\delta}=L_{24}^{\delta}=0.1$, $\delta\mathbf{w}(\mathbf{u})f(t)=(0.1y_{1},0.1y_{2},-0.1x_{1},-0.1x_{2})^{\top}$.

Note that in the second case, the parameters appear in both the observed and
hidden processes. We compare the linear response in the following models:

  * Perfect FDT (or perfect model): The equilibrium PDF $p_{\textrm{eq}}(\mathbf{u})$, the time series $A(\mathbf{u})$ and the formulation of $\mathcal{B}(\mathbf{u})$ are all from the perfect model (39).

  * Imperfect FDT (or imperfect model): The equilibrium PDF $p_{\textrm{eq}}(\mathbf{u})$, the time series $A(\mathbf{u})$ and the formulation of $\mathcal{B}(\mathbf{u})$ are all from the approximate model (41).

  * Concatenate FDT (or concatenate model): The equilibrium PDF $p_{\textrm{eq}}(\mathbf{u})$, the time series $A(\mathbf{u})$ are from the simple concatenation of the observations, i.e., the trajectories of observed variables $x_{1}$ and $x_{2}$ from the perfect model (39), and the trajectories of the hidden variables $y_{1}$ and $y_{2}$ from the approximate model (41) free-run. The formulation of $\mathcal{B}(\mathbf{u})$ is from the approximate model (41).

  * FDT with preconditioner (or preconditioner model): The equilibrium PDF $p_{\textrm{eq}}(\mathbf{u})$ is calculated via (10) where true observations and the CGNS approximate model (41) are utilized in computing $K_{\mathbf{H}}(\mathbf{X}-\mathbf{X}^{\textrm{obs}}(t_{i}))$ and $p^{M|\textrm{obs}}(\mathbf{Y}|\mathbf{X}(s\leq t_{i}))$. The time series $A(\mathbf{u})$ are generated by equation (7) with the CGNS approximate model and the true observations. The formulation of $\mathcal{B}(\mathbf{u})$ is from the perfect model.

The details of the $\mathcal{B}(\mathbf{u})$’s forms from the perfect model
and approximate model can be found in Appendix B. In this experiment, the true
signal is obtained using Euler-Maruyama scheme with a uniform time step
$\delta t=5\times 10^{-3}$ and the total length is 1000 time units.

Before discussing the linear response using CGNS as a preconditioner, we start
by showing some pre-requisites, the equilibrium covariance, and time series
$A(\mathbf{u})$ from different models. The trajectories from the approximate
model (41) free-run (red color) are shown in Figure 13. Note that due to the
violation of the energy-conserving constraint of the nonlinear terms in (41),
the amplitude of $x_{1}$ and $x_{2}$ from the approximate model is much larger
than the one from the perfect model (39). In addition, the PDF of $y_{1}$ is
Gaussian by design, which is also different from the skewed PDF as in the
perfect model. The model error can also be found in the comparison of the
equilibrium covariance matrices of the perfect model (Panel (a)) and the
approximate model (Panel (d)) in Figure 14. Due to large error caused by the
simple parameterization of the hidden processes, one may consider
concatenating the observations (from the perfect model) and the trajectories
of the hidden variables $y_{1}$ and $y_{2}$ from the approximate model (41)
free-run to approximate the required equilibrium PDF $p_{\textrm{eq}}$ and the
unobserved time series. However, it is expected that the correlation between
the observations from the perfect model and trajectories from the approximate
model free-run is neglected. See Panel (c) of Figure 14.

In contrast, in light of the desired structures of CGNS, the partial
observations can be incorporated with both approximating the equilibrium PDF
and recovering the unobserved times series. Therefore, the model error is
significantly mitigated, and the correlation between the observed and
unobserved variables is reserved. The green curves in Figure 13 show one
sampled trajectories of $y_{1}$ and $y_{2}$ using conditional sampling formula
(7) from the approximate model (41) and observations. The overall dynamics of
the recovered sampled trajectories are very similar to those in the perfect
model. In addition, the skewed PDF of $y_{1}$ can be found in the green curve
but the appropriate model free run only brings Gaussian statistics by design.
More importantly, the correlation between the sampled trajectories and the
observations is consistent with that in the perfect model as shown in Panel
(b) of Figure 14. A final remark is that the smoother (or filter) mean time
series is widely used as a surrogate of the true signals, which, however,
underestimate the uncertainty as shown in the orange color in Figure 13.
Therefore, conditional sampling is essential in approximating the time series
of $A(\mathbf{u})$.

Figure 15 shows the linear response operator $R(t)$ in (37) for the response
of the first four moments when perturbing the external forcing parameters in
the observed processes. Here, the blue and red curves show the linear response
from the perfect model (39) and the free-run of the approximate model (41),
respectively. The magenta curves show the linear response from simply
concatenating the trajectories of observed variables $x_{1}$ and $x_{2}$ from
the perfect model (39) and the trajectories of the hidden variables $y_{1}$
and $y_{2}$ from the approximate model (41). The green curves show the linear
response from the procedure discussed in Section 6.2. The imperfect FDT
contains huge errors in capturing higher moments for the observed variables
$x_{1}$ and $x_{2}$. Concatenate FDT works well for computing the four moments
of $x_{1}$; however, it is gradually away from the perfect FDT from lower
moments to higher moments. For example, the gap between the concatenate FDT
and the perfect FDT is obvious in the last row of $x_{2}$. This is because the
strong correlation between the $x_{2}$ and $y_{1}$ (shown in Panel (a) in
Figure 14) is omitted in this simple concatenation (shown in Panel (c) in
Figure 14). In contrast, the response operator $R(t)$ of the two observed
variables from preconditioner FDT is very close to the perfect FDT. Due to the
indirect perturbation of the unresolved variables, the response is expected to
be small. The preconditioner FDT can still capture the trend of the linear
response operator as that using perfect FDT. Figure 16 shows the second
perturbation case, i.e., perturbing the linear interaction parameters
appearing in both the observed and the hidden processes. In addition to the
model error of the approximate model being mitigated and correlation between
the observations and hidden dimensions of the approximate model being
reserved, the perfect model structure is utilized to calculate the
$\mathcal{B}(\mathbf{u})$ in FDT with preconditioner. Therefore, the
performance of the FDT with preconditioner outperforms concatenate FDT and
imperfect FDT.

Figure 13: Comparison of the trajectories from the perfect 4D stochastic
climate model (39), approximate model (41), the smoother mean time series, and
the sampled trajectories based on the approximate model. Blue curves: perfect
model trajectories; red curves: approximate model trajectories; orange curves:
the smoother mean time series; green curves: sampled trajectories. Panel (a):
trajectories; Panel (b): PDFs. Figure 14: Comparison of covariance matrices in
different scenarios. Panel (a): the perfect model where the equilibrium PDF is
from the perfect model (39); Panel (b): preconditioner model where equilibrium
PDF is from (4) where true observations and approximate model are utilized in
computing $K_{\mathbf{H}}(\mathbf{X}-\mathbf{X}^{\textrm{obs}}(t_{i}))$ and
$p^{M|\textrm{obs}}(\mathbf{Y}|\mathbf{X}(s\leq t_{i}))$; Panel (c):
concatenate model where the equilibrium PDF $p_{\textrm{eq}}$ is from the
simple concatenation of the observations, i.e., the trajectories of observed
variables $x_{1}$ and $x_{2}$ from the perfect model (39) and the trajectories
of the hidden variables $y_{1}$ and $y_{2}$ from the approximate model (41)
free-run; Panel (d): the imperfect model where the equilibrium PDF is from the
approximate model (41). Figure 15: Response operator $\mathbf{R}(t)$ in (37)
for the response of the first four moments of the 4D climate model (39) when
perturbing parameters in the observed processes with
$F_{1}^{\delta}=F_{2}^{\delta}=0.3$. In each panel, the blue and red curves
show the linear response from the perfect model (39) and the free-run of the
approximate model (41), respectively. The magenta curves show linear response
from simply concatenating the trajectories of observed variables $x_{1}$ and
$x_{2}$ from the perfect model (39) and the trajectories of the hidden
variables $y_{1}$ and $y_{2}$ from the approximate model (41). The green
curves show the linear response from the procedure discussed in Section 6.2.
Panel (a)–(d): $x_{1}$, $x_{2}$, $y_{1}$, and $y_{2}$, respectively.  Figure
16: Similar to Figure 15 but perturbing parameters in the linear interaction
terms $L_{13}^{\delta}=L_{24}^{\delta}=0.1$. Panel (a)–(d): $x_{1}$, $x_{2}$,
$y_{1}$, and $y_{2}$, respectively.

##  7 Discussion and Conclusions

In this paper, the skill of a rich class of nonlinear stochastic models, known
as the “conditional Gaussian nonlinear system” (CGNS) (2), as both a cheap
surrogate model and a fast preconditioner is explored to advance many
computationally challenging tasks in complex nonlinear systems. The CGNS not
only preserves the main underlying physics of nature but reproduces the
observed intermittency, extreme events and other non-Gaussian features as
well. The closed analytic formula of solving the conditional statistics
facilitates the development of many mathematical theories and fast numerical
algorithms. Three topics are covered in this paper. First, the closed analytic
formulae of the conditional statistics of the CGNS allow an efficient and
accurate data assimilation scheme. It is shown in Section 4 that the data
assimilation skill of a suitable CGNS approximate forecast model outweighs the
EnKBF applying directly even to the perfect model in the presence of strong
nonlinear and turbulent features. The latter may suffer from filter divergence
when the observational process is highly nonlinear with small observational
uncertainties. Second, as is shown in Section 5, the CGNS allows the
development of a fast algorithm for simultaneously estimating the parameters
and the unobserved state variables with uncertainty quantification in the
presence of only partial observations. Utilizing an appropriate CGNS as a
preconditioner significantly reduces the computational cost in accurately
estimating the parameters in the original complex system. The same CGNS can
also serve as a surrogate model for reproducing the large-scale dynamics and
statistics of nature and can be applicable to ensemble forecast. Finally, the
CGNS advances rapid and statistically accurate algorithms for both computing
the probability density function and sampling the trajectories of the
unobserved state variables. Showing in Section 6, these fast algorithms
facilitate the development of an efficient and accurate data-driven method for
predicting the linear response of the original system with respect to
parameter perturbations based on a suitable CGNS preconditioner.

Several important topics are remained as future work. First, it is crucial to
systematically determine the CGNS approximate model. A promising approach is
to write down the general structure of (12) and then apply a parameter
estimation algorithm together with certain sparse identification method to
prevent the overfitting issue. Second, it is of practical significance to
develop further pathways that allow to apply the CGNS to more complicated
systems and explore the approximate errors. Incorporating the CGNS into
intermediate complicated models or even certain versions of the general
circulation models (GCMs), say for climate science or geophysics, can be
interesting and practically useful tasks. Finally, the CGNS may have the
potential to combine with machine learning algorithms to advance the ensemble
forecast. One possible direction is to exploit the analytic formula in (9) for
the ensemble forecast, where the complicated nonlinear interactions in solving
the conditional distributions can be replaced by a cheaper machine learning
architecture.

## Acknowledgments

The research of N.C. is partially funded by the Office of VCRGE at UW-Madison
and ONR N00014-21-1-2904. Y. L. is supported as a graduate student under the
these grants. The research of Y.L. was also supported in part by NSF Award
DMS-2023239 through the IFDS at UW-Madison. The work of H.L. is partially
funded by NSF Award DMS-2108856.

##  Appendix A Details of the EM algorithm for parameters estimation

###  A.1 The EM algorithm for the CGNS

In this appendix, we provide some technical details about Algorithm 5.1, which
concerns the use of CGNS as a preconditioner in the EM approach for parameters
estimation. To fix ideas, we will focus on a special case in which the
original nonlinear system (11) is subject to diagonal and additive noise,
which is sufficient for the applications considered in Section 5. For the
simplicity of discussions, we will also consider only the real-valued
variables and systems. We refer to [18] for more general settings; see also
[48].

We first generate the discrete partial observations $\widehat{\mathbf{X}}$
using for instance the Euler-Maruyama scheme [44] applied to the original
nonlinear system (11), with a sufficiently small time step size $\Delta{t}$.
Given a CGNS approximation of the form (2), we explain now how the conditional
distribution
$p^{M}(\widehat{\mathbf{Y}}|\widehat{\mathbf{X}},\bm{\theta}^{M}_{k-1})$ and
the cost function
$\widetilde{\mathcal{Q}}(\bm{\theta}^{M};\bm{\theta}^{M}_{k-1})$ in lines 5-6
of Algorithm 5.1 are formed at each EM iteration for the CGNS.

Due to the above choice of the noise terms in (11), $\mathbf{B}_{1}$ and
$\mathbf{b}_{2}$ in its corresponding CGNS approximation (2) are simply
diagonal matrices with unknown diffusion coefficient parameters appearing on
the corresponding main diagonal. In the following, we also denote by
$\bm{\xi}$ the parameters in the drift part of the CGNS. Then, the collection
$\bm{\theta}^{M}$ of the parameters in the CGNS is
$\bm{\theta}^{M}=(\bm{\xi},\mathrm{diag}(\mathbf{B}_{1}),\mathrm{diag}(\mathbf{b}_{2}))^{\top}$.
The discretization of (2) using the Euler-Maruyama scheme reads,

$\displaystyle\mathbf{X}^{j+1}=$
$\displaystyle\mathbf{X}^{j}+(\mathbf{A}_{0}(\mathbf{X}^{j},t;\bm{\xi})+\mathbf{A}_{1}(\mathbf{X}^{j},t;\bm{\xi})\mathbf{Y}^{j})\Delta
t+\mathbf{B}_{1}\sqrt{\Delta t}\bm{\varepsilon}^{j}_{1},$ (42a)
$\displaystyle\mathbf{Y}^{j+1}=$
$\displaystyle\mathbf{Y}^{j}+(\mathbf{a}_{0}(\mathbf{X}^{j},t;\bm{\xi})+\mathbf{a}_{1}(\mathbf{X}^{j},t;\bm{\xi})\mathbf{Y}^{j})\Delta
t+\mathbf{b}_{2}\sqrt{\Delta t}\bm{\varepsilon}^{j}_{2},$ (42b)

where $j=0,\ldots,J$ for some fixed positive integer $J$. Here,
$\bm{\varepsilon}^{j}_{1}$ and $\bm{\varepsilon}^{j}_{2}$ are independent and
identically distributed Gaussian white noises. They have the same dimensions
as $\mathbf{X}$ and $\mathbf{Y}$, respectively, due again to the way the noise
forcing is chosen. Assume also all the parameters in the drift part appear as
multiplicative prefactors of some functions of $\mathbf{X}^{j}$ and
$\mathbf{Y}^{j}$.

We introduce now some additional notations in order to put (42) into a more
compact form to be used below. Denote by $\mathbf{M}^{j}$ the matrix that
includes those linear/nonlinear functions in the drift part which are
multiplied by the parameters $\bm{\xi}$. Denote also by $\mathbf{S}^{j}$ those
terms that do not involve parameters such as the first terms $\mathbf{X}^{j}$
or $\mathbf{Y}^{j}$ in (42). Finally, let $\mathbf{R}$ be the covariance
matrix associated with the noise terms in (42); namely, the diagonal matrix
whose diagonal consisting of those from the diagonal matrices
$\mathbf{B}_{1}\mathbf{B}^{\top}_{1}\Delta t$ and
$\mathbf{b}_{2}\mathbf{b}^{\top}_{2}\Delta t$. Apparently, there is a one-to-
one correspondence between the diagonal of $\mathbf{R}$ and the parameters in
the diffusion terms. With these notations we can rewrite (42) into

$\mathbf{u}^{j+1}=\mathbf{M}^{j}\bm{\xi}+\mathbf{S}^{j}+\mathbf{R}^{1/2}\bm{\varepsilon}^{j},\quad
j=0,\ldots,J,$ (43)

where $\mathbf{u}^{j+1}=(\mathbf{X}^{j+1},\mathbf{Y}^{j+1})^{\top}$ and
$\bm{\varepsilon}^{j}=(\bm{\varepsilon}_{1}^{j},\bm{\varepsilon}_{2}^{j})^{\top}$.
Thus, at each time step, given $\mathbf{M}^{j}$ and $\mathbf{S}^{j}$,
$\mathbf{u}^{j+1}$ follows a Gaussian distribution with mean
$\bm{\mu}^{j}=\mathbf{M}^{j}\bm{\xi}+\mathbf{S}^{j}$ and variance given by
$\mathbf{R}$:

$p^{M}(\mathbf{u}^{j+1}|\mathbf{M}^{j},\mathbf{S}^{j},\bm{\theta}^{M})=\frac{1}{\sqrt{(2\pi)^{N}}}|\mathbf{R}|^{-\frac{1}{2}}\exp\Big{(}-\frac{1}{2}(\mathbf{u}^{j+1}-\bm{\mu}^{j})^{\top}(\mathbf{R})^{-1}(\mathbf{u}^{j+1}-\bm{\mu}^{j})\Big{)},$
(44)

where $N$ is the dimension of the phase space.

At the $k$-th EM iteration for each $k=1,\ldots,K$, the parameters
$\bm{\theta}^{M}_{k-1}$ is already computed. We can thus use (6) to compute
the optimal smoother estimate
$p^{M}(\mathbf{Y}(t)|\mathbf{X}(s),s\in[0,T],\bm{\theta}^{M}_{k-1})$ or
equivalently its “discretized” form
$p^{M}(\mathbf{Y}^{j}|\widehat{\mathbf{X}},\bm{\theta}^{M}_{k-1})$ for
$j=0,\ldots,J$ In the E-Step. On the other hand, by exploiting the
relationship in (43) for all $j$, the M-Step (line 6 in Algorithm 5.1) is
solved via minimizing the following cost function:

$\widetilde{\mathcal{L}}=\frac{1}{2}\sum_{j=J_{1}}^{J-1}\left\langle(\mathbf{u}^{j+1}-\mathbf{M}^{j}\bm{\xi}-\mathbf{S}^{j})^{\top}(\mathbf{R})^{-1}(\mathbf{u}^{j+1}-\mathbf{M}^{j}\bm{\xi}-\mathbf{S}^{j})\right\rangle+\frac{J^{\prime}}{2}\log|\mathbf{R}|,$
(45)

where the summation of $j$ starts from a certain non-zero integer $J_{1}$ to
eliminate the inaccuracy from the burn-in period and $J^{\prime}=J-J_{1}$.
Note that (45) corresponds to the negative of
$\mathcal{Q}(\bm{\theta};\bm{\theta}_{k})$ defined by (27) if $J_{1}=1$, after
dropping some constant terms independent of $\bm{\xi}$ and $\mathbf{R}$. In
(45), $\left\langle\cdot\right\rangle$ denotes the expectation over the
uncertain component of $\mathbf{u}^{j}$ and $\mathbf{u}^{j+1}$, namely
$\mathbf{Y}^{j}$ and $\mathbf{Y}^{j+1}$ while the expectations of the observed
component $\mathbf{X}^{j}$ and $\mathbf{X}^{j+1}$ are simply themselves since
the time series of them are given. Since the hidden variables $\mathbf{Y}^{j}$
appears in a linear way in the matrix $\mathbf{M}^{j}$ due to the structure of
the CGNS, only the quadratic terms of $\mathbf{Y}^{j}$, namely
$\langle\mathbf{Y}^{j+1},(\mathbf{Y}^{j+1})^{\top}\rangle$,
$\langle\mathbf{Y}^{j+1},(\mathbf{Y}^{j})^{\top}\rangle$,
$\langle\mathbf{Y}^{j},(\mathbf{Y}^{j})^{\top}\rangle$, need to be solved in
the expectation in (45). These terms can be solved via some manipulations of
the results from the closed formulae of the smoother estimates (6). Details
can be found in [18, Appendix A]. To find the minimum of
$\widetilde{\mathcal{L}}$, we set
$\frac{\partial\widetilde{\mathcal{L}}}{\partial\xi_{i}}=0$ and
$\frac{\partial\widetilde{\mathcal{L}}}{\partial R_{\ell\ell}}=0$ for each
component $\xi_{i}$ of $\bm{\xi}$ and each diagonal element $R_{\ell\ell}$ of
$\mathbf{R}$, which leads to

$\displaystyle\mathbf{R}$
$\displaystyle=\frac{1}{J^{\prime}}\sum_{j=J_{1}}^{J-1}\left\langle(\mathbf{u}^{j+1}-\mathbf{M}^{j}\bm{\xi}-\mathbf{S}^{j})(\mathbf{u}^{j+1}-\mathbf{M}^{j}\bm{\xi}-\mathbf{S}^{j})^{\top}\right\rangle,$
(46a) $\displaystyle\bm{\xi}$ $\displaystyle=\mathbf{D}^{-1}\mathbf{c},$ (46b)

where

$\mathbf{D}=\sum_{j=J_{1}}^{J-1}\left\langle(\mathbf{M}^{j})^{\top}\mathbf{R}^{-1}\mathbf{M}^{j}\right\rangle\quad\text{and}\quad\mathbf{c}=\sum_{j=J_{1}}^{J-1}\left\langle(\mathbf{M}^{j})^{\top}\mathbf{R}^{-1}(\mathbf{u}^{j+1}-\mathbf{S}^{j})\right\rangle.$
(47)

Note that we solve (46) based on an iteration method where $\bm{\xi}$ is
obtained given $\mathbf{R}$ from the previous step, and then $\mathbf{R}$ is
calculated from the updated $\bm{\xi}$.

###  A.2 The last M-Step in Algorithm 5.1

Recall Algorithm 5.1. The last M-Step (line 8 in the algorithm) requires to
solve the minimization of the cost function, which is based on the original
nonlinear system. The main difference here compared with minimizing the cost
function associated with the CGNS (line 6 in the algorithm) is that higher
order moments of $\mathbf{Y}$, resulting from the general nonlinear structure
of the original complex system, may be involved. If the nonlinearity of the
original complex system is up to quadratic, then the expectation in the
analogue of (45) for the original system may involve up to the fourth moments
of $\mathbf{Y}$.

Note that these moments are calculated based on the smoother estimate from the
E-Step (line 7), which still utilizes the CGNS. In other words, the smoother
estimate only provides a conditional Gaussian distribution. The higher order
moments are thus computed based on the quasi-Gaussian closure approximation,
which are represented by the known information from the mean and the variance
of the conditional Gaussian smoother estimate. For example, denote by $Y_{i}$
a scalar component of $\mathbf{Y}$. So does $Y_{j},Y_{k}$ and $Y_{m}$. Then
the third order moment $\langle Y_{i}Y_{j}Y_{k}\rangle$ and the fourth order
moment $\langle Y_{i}Y_{j}Y_{k}Y_{m}\rangle$ can be obtained as follows

$\displaystyle\langle Y_{i}Y_{j}Y_{k}\rangle$
$\displaystyle=\mu_{i}\mu_{j}\mu_{k}+\mu_{k}\sigma_{ij}+\mu_{i}\sigma_{kj}+\mu_{j}\sigma_{ik},$
(48) $\displaystyle\langle Y_{i}Y_{j}Y_{k}Y_{m}\rangle$ $\displaystyle=\langle
Y_{i}Y_{j}Y_{k}\rangle\mu_{m}+\mu_{i}\mu_{j}\sigma_{km}+\mu_{i}\mu_{k}\sigma_{jm}+\mu_{k}\mu_{j}\sigma_{im}+\sigma_{ij}\sigma_{km}+\sigma_{ik}\sigma_{jm}+\sigma_{jk}\sigma_{im},$

where $\mu_{i}$ and $\sigma_{ij}$ are mean and covariance of corresponding
components.

##  Appendix B Calculating $\mathcal{B}(\mathbf{u})$ in the linear response
operator

In light of (37) and (38), one has the following explicit expression of
$\mathcal{B}(u)$

$\displaystyle\mathcal{B}(\mathbf{u})$
$\displaystyle=-\frac{\textrm{div}_{\mathbf{u}}(\mathbf{w}(\mathbf{u})p_{\textrm{eq}}^{M|\textrm{obs}}(\mathbf{u}))}{p_{\textrm{eq}}^{M|\textrm{obs}}(\mathbf{u})}$
(49)
$\displaystyle=-\sum_{i=1}^{N}\frac{\partial}{\partial{\mathbf{u}_{i}}}\mathbf{w}_{i}(\mathbf{u})-\sum_{i=1}^{N}\mathbf{w}_{i}\frac{\partial}{\partial\mathbf{u}_{i}}p_{\textrm{eq}}^{M|\textrm{obs}}(\mathbf{u})\,.$

When perturbing the parameters of forcing in the observed processes, the forms
of $\mathcal{B}(\mathbf{u})$ from the perfect model and the approximate model
are the same, since the parameters $F_{1}$ and $F_{2}$ appear exactly the same
way as in both the perfect and the approximate model. The
$\mathcal{B}(\mathbf{u})$ term reads as follows

$\mathcal{B}(\mathbf{u})=-\frac{\partial}{\partial
x_{1}}p_{\textrm{eq}}^{M|\textrm{obs}}(\mathbf{u})-\frac{\partial}{\partial
x_{2}}p_{\textrm{eq}}^{M|\textrm{obs}}(\mathbf{u})\,.$ (50)

When perturbing the parameters in linear interactions terms that appear in
both the observed and hidden processes, the formulation of
$\mathcal{B}(\mathbf{u})$ from the perfect and approximate models are
different. Given the perturbation vector
$\mathbf{w}(\mathbf{u})=(y_{1},y_{2},-x_{1},-x_{2})^{\top}$, the
$\mathcal{B}(\mathbf{u})$ from the perfect model is as follows

$\mathcal{B}(\mathbf{u})=-y_{1}\frac{\partial}{\partial
x_{1}}p_{\textrm{eq}}^{M|\textrm{obs}}(\mathbf{u})-y_{2}\frac{\partial}{\partial
x_{2}}p_{\textrm{eq}}^{M|\textrm{obs}}(\mathbf{u})+x_{1}\frac{\partial}{\partial
y_{1}}p_{\textrm{eq}}^{M|\textrm{obs}}(\mathbf{u})+x_{2}\frac{\partial}{\partial
y_{2}}p_{\textrm{eq}}^{M|\textrm{obs}}(\mathbf{u})\,.$ (51)

However, since there is no $L_{13}$ and $L_{24}$ parameters in the hidden
processes of the approximate model, the formulation of
$\mathcal{B}(\mathbf{u})$ from the approximate model remains

$\mathcal{B}(\mathbf{u})=-y_{1}\frac{\partial}{\partial
x_{1}}p_{\textrm{eq}}^{M|\textrm{obs}}(\mathbf{u})-y_{2}\frac{\partial}{\partial
x_{2}}p_{\textrm{eq}}^{M|\textrm{obs}}(\mathbf{u}).$ (52)

## References

  * [1] Shady E Ahmed, Suraj Pawar, Omer San, Adil Rasheed, Traian Iliescu, and Bernd R Noack.  On closures for reduced order models–A spectrum of first-principle to machine-learned avenues.  Physics of Fluids, 33(9):091301, 2021.
  * [2] HM Arnold, IM Moroz, and TN Palmer.  Stochastic parametrizations and model uncertainty in the Lorenz’96 system.  Phil. Trans. R. Soc. A, 371(1991):20110479, 2013.
  * [3] Mark Asch, Marc Bocquet, and Maëlle Nodet.  Data assimilation: methods, algorithms, and applications.  SIAM, 2016.
  * [4] K. Bergemann and S. Reich.  An ensemble Kalman-Bucy filter for continuous data assimilation.  Meteorologische Zeitschrift, 21:213–219, 2012.
  * [5] Zdravko I Botev, Joseph F Grotowski, and Dirk P Kroese.  Kernel density estimation via diffusion.  The annals of Statistics, 38(5):2916–2957, 2010.
  * [6] Michal Branicki, Nan Chen, and Andrew J Majda.  Non-Gaussian test models for prediction and state estimation with model errors.  Chinese Annals of Mathematics, Series B, 34(1):29–64, 2013.
  * [7] Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 113(15):3932–3937, 2016.
  * [8] Kevin Carlberg, Charbel Farhat, Julien Cortial, and David Amsallem.  The GNAT method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows.  Journal of Computational Physics, 242:623–647, 2013.
  * [9] James Carpenter, Peter Clifford, and Paul Fearnhead.  Improved particle filter for nonlinear problems.  IEE Proceedings-Radar, Sonar and Navigation, 146(1):2–7, 1999.
  * [10] Ashesh Chattopadhyay, Mustafa Mustafa, Pedram Hassanzadeh, Eviatar Bach, and Karthik Kashinath.  Towards physically consistent data-driven weather forecasting: Integrating data assimilation with equivariance-preserving spatial transformers in a case study with ERA5.  Geoscientific Model Development Discussions, pages 1–23, 2021.