#  Nambu-Jona-Lasinio model correlation functions from QCD

Marco Frasca111Speaker, Corresponding author<EMAIL_ADDRESS>Via Erasmo
Gattamelata, 3, 00176 Rome (Italy)  Anish Ghoshal  INFN, Rome, Italy and
Warsaw University, Poland  Stefan Groote  University of Tartu, Estonia

###### Abstract

We treat quantum chromodynamics (QCD) using a set of Dyson-Schwinger equations
derived, in differential form, with the Bender-Milton-Savage technique. In
this way, we are able to derive the low energy limit that assumes the form of
a non-local Nambu-Jona-Lasinio model with all the parameters properly fixed by
the QCD Lagrangian and the determination of the mass gap of the gluon sector.

###### keywords:

##  1 Introduction

Adding quarks to the Yang-Mills Lagrangian makes the theory not exactly
treatable. Notwithstanding such a difficulty, full QCD can be handled with
Dyson-Schwinger equations even if some approximations are needed to get the
low-energy limit. Such approximations entail both the strong coupling limit
and ’t Hooft limit $N\rightarrow\infty,\ \ Ng^{2}=constant\gg 1$. In this way,
we will recover, in the low-energy limit, the equations for the correlation
functions of a non-local NJL model.

Yang-Mills set of Dyson-Schwinger equations can be solved through a class of
exact solutions of the 1-point function, similarly to the $\phi^{4}$ theory
[1]. The spectrum can be obtained as well and appears in agreement with
lattice data [2]. Adding quarks to the theory makes it unsolvable but, also in
this case, our approach could be proven meaningful to derive the proper low
energy limit of QCD.

Our idea is to provide a method to derive the full hierarchy of Dyson-
Schwinger equations also for QCD, retaining their full differential form. This
is possible provided we use a technique devised by Bender, Milton and Savage
[3].

Some approximations are needed to get the low-energy limit. Such
approximations entail both the strong coupling limit and ’t Hooft limit
$N\rightarrow\infty$ and $Ng^{2}=constant\gg 1$ [4, 5]. In this way, we will
recover, in the low-energy limit, the equations for the correlation functions
of a non-local Nambu-Jona-Lasinio model.

This is an interesting result in view of the fact that it permits to recover
from an error present in Ref. [6] granting anyway the conclusions. This has
been recently applied to the $g-2$ problem [7] to evaluate the hadronic vacuum
polarization for the $\pi\pi$ contribution.

##  2 Bender-Milton-Savage technique

This technique can be better explained referring to a scalar field. We
consider the following partition function

$Z[j]=\int[D\phi]e^{iS(\phi)+i\int d^{4}xj(x)\phi(x)}.$ (1)

To derive 1P-function, one has

$\left\langle\frac{\delta S}{\delta\phi(x)}\right\rangle=j(x),$ (2)

assuming

$\left\langle\ldots\right\rangle=\frac{\int[D\phi]\ldots e^{iS(\phi)+i\int
d^{4}xj(x)\phi(x)}}{\int[D\phi]e^{iS(\phi)+i\int d^{4}xj(x)\phi(x)}}$ (3)

By setting $j=0$ one obtains the equation for th 1P-function. Next, we derive
this equation again with respect to $j$ to get the equation for the
2P-function. We are taking for the nP-functions the following definition

$\langle\phi(x_{1})\phi(x_{2})\ldots\phi(x_{n})\rangle=\frac{\delta^{n}\ln(Z[j])}{\delta
j(x_{1})\delta j(x_{2})\ldots\delta j(x_{n})}.$ (4)

This implies

$\frac{\delta G_{k}(\ldots)}{\delta j(x)}=G_{k+1}(\ldots,x).$ (5)

Such a procedure can be iterated to whatever order giving, in principle, all
the hierarchy of the Dyson-Schwinger equations in PDE form. Going to higher
orders could imply complicated computations but this approach shows itself to
be very useful when some known solutions are given for 1P- and 2P-functions as
in our case.

##  3 1P and 2P functions for QCD

We work in the Landau gauge as can make some computations simpler as seen in
[1].

The Bender-Milton-Savage method yields for the 1P-functions

$\displaystyle\partial^{2}G_{1\nu}^{a}(x)+gf^{abc}(\partial^{\mu}G_{2\mu\nu}^{bc}(0)+$
$\displaystyle\partial^{\mu}G_{1\mu}^{b}(x)G_{1\nu}^{c}(x)-\partial_{\nu}G_{2\mu}^{\nu
bc}(0)$ $\displaystyle-\partial_{\nu}G_{1\mu}^{b}(x)G_{1}^{\mu c}(x))$
$\displaystyle+gf^{abc}\partial^{\mu}G_{2\mu\nu}^{bc}(0)+gf^{abc}\partial^{\mu}(G_{1\mu}^{b}(x)G_{1\nu}^{c}(x))$
$\displaystyle+g^{2}f^{abc}f^{cde}(G_{3\mu\nu}^{\mu
bde}(0,0)+G_{2\mu\nu}^{bd}(0)G_{1}^{\mu e}(x)$
$\displaystyle+G_{2\nu\rho}^{eb}(0)G_{1}^{\rho
d}(x)+G_{2\mu\nu}^{de}(0)G_{1}^{\mu b}(x)+$ $\displaystyle G_{1}^{\mu
b}(x)G_{1\mu}^{d}(x)G_{1\nu}^{e}(x))$
$\displaystyle=g\sum_{q,i}\gamma_{\nu}T^{a}S_{q}^{ii}(0)+g\sum_{q,i}{\bar{q}}_{1}^{i}(x)\gamma_{\nu}T^{a}q_{1}^{i}(x),$
(6)

and for the quarks

$(i\not{\partial}-m_{q})q_{1}^{i}(x)+g{\bm{T}}\cdot\not{\bm{G}}_{1}(x)q_{1}^{i}(x)+g{\bm{T}}\cdot\not{\bm{W}}^{i}_{q}(x,x)=0.$
(7)

Here and in the following Greek indexes ($\mu,\nu,\ldots$) are for the space-
time and Latin index ($a,b,\ldots$) for the gauge group. A usual for the
Dyson-Schwinger set, the equations for the lower order nP-functions depend on
the higher order correlation functions. We will see how to treat this aspect
in the following.

At this stage, we apply the re-mapping idea to such equations as done in [1].
So, we assume

$G_{1\nu}^{a}(x)\rightarrow\eta_{\nu}^{a}\phi(x)$ (8)

being $\phi(x)$ a scalar field. Let us introduce the $\eta$-symbols as follows

$\displaystyle\eta_{\mu}^{a}\eta^{a\mu}$ $\displaystyle=$ $\displaystyle
N^{2}-1.$ $\displaystyle\eta_{\mu}^{a}\eta^{b\mu}$ $\displaystyle=$
$\displaystyle\delta_{ab},$ $\displaystyle\eta_{\mu}^{a}\eta_{\nu}^{a}$
$\displaystyle=$ $\displaystyle\left(g_{\mu\nu}-\delta_{\mu\nu}\right)/2.$ (9)

All this permits to get the reduced equations

$\displaystyle\partial^{2}\phi(x)+2Ng^{2}\Delta(0)\phi(x)+Ng^{2}\phi^{3}(x)$
$\displaystyle=\frac{1}{N^{2}-1}\left[g\sum_{q,i}\eta^{a\nu}\gamma_{\nu}T^{a}S_{q}^{ii}(0)\right.$
$\displaystyle\left.+g\sum_{q,i}{\bar{q}}_{1}^{i}(x)\eta^{a\nu}\gamma_{\nu}T^{a}q_{1}^{i}(x)\right]$
$\displaystyle(i\not{\partial}-m_{q}^{i})q_{1}^{i}(x)+g{\bm{T}}\cdot\not{\eta}\phi(x)q_{1}^{i}(x)=0.$
(10)

We do the same for the 2P-functions. In the Landau gauge, the gluon
2P-function takes the form

$G_{2\mu\nu}^{ab}(x-y)=\left(\eta_{\mu\nu}-\frac{\partial_{\mu}\partial_{\nu}}{\partial^{2}}\right)\Delta_{\phi}(x-y)$
(11)

being $\eta_{\mu\nu}$ the Minkowski metric, and $\Delta_{\phi}(x-y)$ is the
propagator of the $\phi$ given the map between the scalar and the Yang-Mills
fields. Finally, we can write

$\displaystyle\partial^{2}\Delta_{\phi}(x-y)+2Ng^{2}\Delta_{\phi}(0)\Delta_{\phi}(x-y)+3Ng^{2}\phi^{2}(x)\Delta_{\phi}(x-y)$
$\displaystyle=g\sum_{q,i}{\bar{Q}}^{ia}_{\nu}(x-y)\gamma^{\nu}T^{a}q_{1}^{i}(x)$
$\displaystyle+g\sum_{q,i}{\bar{q}}_{1}^{i}(x)\gamma^{\nu}T^{a}Q^{ia}_{\nu}(x-y)+\delta^{4}(x-y)$
$\displaystyle\partial^{2}P^{ad}_{2}(x-y)=\delta_{ad}\delta^{4}(x-y)$
$\displaystyle(i\not{\partial}-m_{q}^{i})S^{ij}_{q}(x-y)$
$\displaystyle+g{\bm{T}}\cdot\not{\eta}\phi(x)S^{ij}_{q}(x-y)=\delta_{ij}\delta^{4}(x-y)$
$\displaystyle\partial^{2}W_{q\nu}^{ai}(x-y)+2Ng^{2}\Delta_{\phi}(0)W_{q\nu}^{ai}(x-y)+3Ng^{2}\phi^{2}(x)W_{q\nu}^{ai}$
$\displaystyle=g\sum_{j}{\bar{q}}_{1}^{j}(x)\gamma_{\nu}T^{a}S^{ji}_{q}(x-y)$
$\displaystyle(i\not{\partial}-m_{q}^{i})Q^{ia}_{\mu}(x-y)+g{\bm{T}}\cdot\not{\eta}\phi(x)Q^{ia}_{\mu}(x-y)$
$\displaystyle+gT^{a}\gamma_{\mu}\Delta_{\phi}(x-y)q_{1}^{i}(x)=0.$ (12)

##  4 ’t Hooft limit

’t Hooft limit means to solve the theory assuming [4, 5]

$N\rightarrow\infty,\qquad Ng^{2}=constant,\qquad Ng^{2}\gg 1.$ (13)

In our case, the gauge group is SU(N) being $N$ is the number of colors. To
evaluate our equations in such a limit, we need a perturbation series for a
very large coupling. We proposed such a technique in Ref.[8]. We do a
rescaling, $x\rightarrow\sqrt{Ng^{2}}x$, and write the equation for the gluon
field as follows

$\displaystyle\partial^{2}\phi(x^{\prime})+2\Delta_{\phi}(0)\phi(x^{\prime})+3\phi^{3}(x^{\prime})$
$\displaystyle=$ (14)
$\displaystyle\frac{1}{\sqrt{Ng^{2}}\sqrt{N}(N^{2}-1)}\left[\sum_{q,i}\eta\cdot\gamma\cdot
TS_{q}^{ii}(0)+\right.$
$\displaystyle\left.\sum_{q,i}{\bar{q}}_{1}^{i}(x^{\prime})\eta\cdot\gamma\cdot
Tq_{1}^{i}(x^{\prime})\right].$

In the ’t Hooft limit, we get at the leading order for the 1P-functions

$\displaystyle\partial^{2}\phi_{0}(x)+2Ng^{2}\Delta_{\phi}(0)\phi_{0}(x)+3Ng^{2}\phi_{0}^{3}(x)=0,$
$\displaystyle(i\not{\partial}-m_{q}^{i}){\hat{q}}_{1}^{i}(x)+g{\bm{T}}\cdot\not{\eta}\phi(x)q_{1}^{i}(x)=0.$
(15)

At the leading order the only effect is seen on masses. For the quark field,
this will be clearer below. We can solve the equation for the gluon field
taking

$\displaystyle\phi_{0}(x)=\sqrt{\frac{2\mu^{4}}{m^{2}+\sqrt{m^{4}+2Ng^{2}\mu^{4}}}}\times$
$\displaystyle{\rm sn}\left(p\cdot x+\chi,\kappa\right),$ (16)

being sn a Jacobi elliptical function, $\mu$ and $\chi$ arbitrary integration
constants and $m^{2}=2Ng^{2}\Delta_{\phi}(0)$. We have

$\kappa=\frac{-m^{2}+\sqrt{m^{4}+2Ng^{2}\mu^{4}}}{-m^{2}-\sqrt{m^{4}+2Ng^{2}\mu^{4}}}.$
(17)

This is true provided that the following dispersion relation holds

$p^{2}=m^{2}+\frac{Ng^{2}\mu^{4}}{m^{2}+\sqrt{m^{4}+2Ng^{2}\mu^{4}}}.$ (18)

For the equations of the 2P-functions one has

$\displaystyle\partial^{2}\Delta_{\phi}(x,y)+2Ng^{2}\Delta_{\phi}(0)\Delta(x-y)+3Ng^{2}\phi_{0}^{2}(x)\Delta_{\phi}(x-y)$
$\displaystyle=g\sum_{q,i}{\bar{Q}}^{ia}_{\nu}(x,y)\gamma^{\nu}T^{a}{\hat{q}}_{1}^{i}(x)$
$\displaystyle+g\sum_{q,i}{\bar{\hat{q}}}_{1}^{i}(x)\gamma^{\nu}T^{a}Q^{ia}_{\nu}(x,y)+\delta^{4}(x-y)$
$\displaystyle\partial^{2}P^{ad}_{2}(x-y)=\delta_{ad}\delta^{4}(x-y)$
$\displaystyle(i\not{\partial}-m_{q}^{i}){\hat{S}}^{ij}_{q}(x-y)+g{\bm{T}}\cdot\not{\eta}\phi(x)S^{ij}_{q}(x-y)=\delta_{ij}\delta^{4}(x-y)$
$\displaystyle\partial^{2}W_{q\nu}^{ai}(x,y)+2Ng^{2}\Delta_{\phi}(0)W_{q\nu}^{ai}(x,y)+3Ng^{2}\phi_{0}^{2}(x)W_{q\nu}^{ai}(x,y)$
$\displaystyle=g\sum_{j}{\bar{\hat{q}}}_{1}^{j}(x)\gamma_{\nu}T^{a}{\hat{S}}^{ji}(x-y)$
$\displaystyle(i\not{\partial}-{\hat{M}}_{q}^{i}){\hat{Q}}^{ia}_{\mu}(x,y)+gT^{a}\gamma_{\mu}\Delta_{\phi}(x-y){\hat{q}}_{1}^{i}(x)=0.$
(19)

To solve these equations, let us consider

$\displaystyle\partial^{2}\Delta_{0}(x-y)+[m^{2}+3Ng^{2}\phi_{0}^{2}(x)]\Delta_{0}(x-y)$
$\displaystyle=$ $\displaystyle\delta^{4}(x-y).$ (20)

In momenta space, the solution of this equation is given by [1, 8]

$\displaystyle\Delta_{0}(p)=M{\hat{Z}}(\mu,m,Ng^{2})\frac{2\pi^{3}}{K^{3}(\kappa)}\times$
$\displaystyle\sum_{n=0}^{\infty}(-1)^{n}\frac{e^{-(n+\frac{1}{2})\pi\frac{K^{\prime}(\kappa)}{K(\kappa)}}}{1-e^{-(2n+1)\frac{K^{\prime}(\kappa)}{K(\kappa)}\pi}}\times$
$\displaystyle(2n+1)^{2}\frac{1}{p^{2}-m_{n}^{2}+i\epsilon}$ (21)

being

$M=\sqrt{m^{2}+\frac{Ng^{2}\mu^{4}}{m^{2}+\sqrt{m^{4}+2Ng^{2}\mu^{4}}}},$ (22)

and ${\hat{Z}}(\mu,m,Ng^{2})$ a given constant. The spectrum is given by
$m_{n}$ and a proper gap equation [2].

##  5 QCD in the low energy limit

Using the technique devised in [8], the next-to-leading order term is given by

$\displaystyle\phi_{1}(x)=g\frac{1}{N^{2}-1}\int
d^{4}x^{\prime}\Delta_{0}(x-x^{\prime})\left[\sum_{q,i}\eta\cdot\gamma\cdot
TS_{q}^{ii}(0)\right.$
$\displaystyle\left.+\sum_{q,i}{\bar{q}}_{1}^{i}(x^{\prime})\eta\cdot\gamma\cdot
Tq_{1}^{i}(x^{\prime})\right].$ (23)

Given equation for the quark 1P-function

$(i\not{\partial}-m_{q})q_{1}^{i}(x)+g{\bm{T}}\cdot\not{\eta}\phi(x)q_{1}^{i}(x)=0,$
(24)

one has

$\displaystyle(i\not{\partial}-m_{q})q_{1}^{i}(x)+g{\bm{T}}\cdot\not{\eta}\phi_{0}(x)q_{1}^{i}(x)$
$\displaystyle+g^{2}\frac{1}{N^{2}-1}\int
d^{4}x^{\prime}\Delta_{0}(x-x^{\prime})\sum_{q,k}{\bar{q}}_{1}^{k}(x^{\prime}){\bm{T}}\cdot\not{\eta}q_{1}^{k}(x^{\prime})\times$
$\displaystyle{\bm{T}}\cdot\not{\eta}q_{1}^{i}(x)=0.$ (25)

’t Hooft limit implies that the $\phi_{0}$ term can be neglected with respect
to the second one and we get

$\displaystyle(i\not{\partial}-m_{q})q_{1}^{i}(x)$
$\displaystyle+g^{2}\frac{1}{N^{2}-1}\int
d^{4}x^{\prime}\Delta_{0}(x-x^{\prime})\sum_{q,k}{\bar{q}}_{1}^{k}(x^{\prime}){\bm{T}}\cdot\not{\eta}q_{1}^{k}(x^{\prime})\times$
$\displaystyle{\bm{T}}\cdot\not{\eta}q_{1}^{i}(x)=0.$ (26)

For the quark propagator one has instead

$(i\not{\partial}-m_{q})S^{ij}_{q}(x-y)+g{\bm{T}}\cdot\not{\eta}\phi(x)S^{ij}_{q}(x-y)=\delta_{ij}\delta^{4}(x-y)$
(27)

therefore

$\displaystyle(i\not{\partial}-m_{q})S^{ij}_{q}(x-y)$
$\displaystyle+g{\bm{T}}\cdot\not{\eta}\phi_{0}(x)S^{ij}_{q}(x-y)$
$\displaystyle+g^{2}\frac{1}{N^{2}-1}\int
d^{4}x^{\prime}\Delta_{0}(x-x^{\prime})\sum_{q,k}{\bar{q}}_{1}^{k}(x^{\prime}){\bm{T}}\cdot\not{\eta}q_{1}^{k}(x^{\prime})\times$
$\displaystyle{\bm{T}}\cdot\not{\eta}S^{ij}_{q}(x-y)=\delta_{ij}\delta^{4}(x-y).$
(28)

Again, by the ’t Hooft limit we can neglect the $\phi_{0}$ term with respect
to the second one and we get

$\displaystyle(i\not{\partial}-m_{q})S^{ij}_{q}(x-y)$
$\displaystyle+g^{2}\frac{1}{N^{2}-1}\int
d^{4}x^{\prime}\Delta_{0}(x-x^{\prime})\sum_{q,k}{\bar{q}}_{1}^{k}(x^{\prime}){\bm{T}}\cdot\not{\eta}q_{1}^{k}(x^{\prime})\times$
$\displaystyle{\bm{T}}\cdot\not{\eta}S^{ij}_{q}(x-y)$
$\displaystyle=\delta_{ij}\delta^{4}(x-y).$ (29)

We can recognize here the equations for the 1P- and 2P-functions of a non-
local Nambu-Jona-Lasinio model. These are not generally treatable. They should
be solved straightforwardly being already quantum averaged. In order to obtain
a gap equation, we need to recover the Nambu-Jona-Lasinio-model from which
they can be obtained doing some kind of backtracking. Only in this way a gap
equation is derived. Indeed, such a model has the Lagrangian

$\displaystyle
L_{NJL}=\sum_{i,q}\left[{\bar{q}}_{i}(x)(i\not{\partial}-m_{q})q_{i}(x)+\right.$
$\displaystyle\left.\frac{g^{2}}{N^{2}-1}{\bar{q}}_{i}(x)\int
d^{4}x^{\prime}\Delta_{0}(x-x^{\prime})\sum_{k,q^{\prime}}{\bar{q}^{\prime}}_{k}(x^{\prime}){\bm{T}}\cdot\not{\eta}q^{\prime}_{k}(x^{\prime}){\bm{T}}\cdot\not{\eta}q_{i}(x)\right].$
(30)

From this, one can get a quark gap equation that is identical to the one given
in [6] as proven in [7].

##  6 Conclusions

We derived the set of Dyson-Schwinger equations, till to 2P-functions, for QCD
with the Bender-Milton-Savage technique. We treated them in the ’t Hooft
limit. The low-energy limit is a nonlocal-Nambu-Jona-Lasinio model. This was
shown obtaining the corresponding 1P- and 2P-equations for its correlation
functions. The Lagrangian of the model is also given.

Work is ongoing to get the gap equation and to analyze its properties in view
of recent g-2 Fermilab measurement [7].

## Acknowledgements

The research was supported in part by the European Regional Development Fund
under Grant No. TK133.

## References

  * [1] M. Frasca, Eur. Phys. J. Plus 132, no.1, 38 (2017) [erratum: Eur. Phys. J. Plus 132, no.5, 242 (2017)] doi:10.1140/epjp/i2017-11321-4 [arXiv:1509.05292 [math-ph]].
  * [2] M. Frasca, Nucl. Part. Phys. Proc. 294-296, 124-128 (2018) doi:10.1016/j.nuclphysbps.2018.02.005 [arXiv:1708.06184 [hep-ph]].
  * [3] C. M. Bender, K. A. Milton and V. Savage, Phys. Rev. D 62, 085001 (2000) doi:10.1103/PhysRevD.62.085001 [arXiv:hep-th/9907045 [hep-th]].
  * [4] G. ’t Hooft, Nucl. Phys. B 72, 461 (1974) doi:10.1016/0550-3213(74)90154-0
  * [5] G. ’t Hooft, Nucl. Phys. B 75, 461-470 (1974) doi:10.1016/0550-3213(74)90088-1
  * [6] M. Frasca, Eur. Phys. J. C 80, no.8, 707 (2020) doi:10.1140/epjc/s10052-020-8261-7 [arXiv:1901.08124 [hep-ph]].
  * [7] M. Frasca, A. Ghoshal and S. Groote, [arXiv:2109.05041 [hep-ph]].
  * [8] M. Frasca, Eur. Phys. J. C 74, 2929 (2014) doi:10.1140/epjc/s10052-014-2929-9 [arXiv:1306.6530 [hep-ph]].


# Stochastic Parametrization using Compressed Sensing: Application to the
Lorenz-96 Atmospheric Model

A. MUKHERJEE1, Y. AYDOGDU1, T. RAVICHANDRAN1 and N. SRI NAMACHCHIVAYA1
1Department of Applied Mathematics, University of Waterloo, Waterloo, Ontario
N2L 3G1, Canada Corresponding author. e-mail<EMAIL_ADDRESS>
###### Abstract

Growing set of optimization and regression techniques, based upon sparse
representations of signals, to build models from data sets has received
widespread attention recently with the advent of compressed sensing. In this
paper, sparse approximations in high-dimensional spaces are used to build
models (vector fields) to emulate the behavior of the fine-scale process, so
that explicit simulations become an online benchmark for parameterization.
Observations are assimilated during the integration of low-dimensional built
model to provide predictions. We outline how the parameterization schemes
developed here and the low-dimensional filtering algorithm can be applied to
the Lorenz-96 Lorenz1996 (Lorenz, 1995) atmospheric model that mimics mid-
latitude atmospheric dynamics with microscopic convective processes.

keywords: compressed sensing, sparse regression, ensemble Kalman filter, auto-
regression

##  1 Introduction

Global weather and climate models are the best available tools for projecting
likely climate changes. However, current climate model projections are still
considerably uncertain, due to uncertainties and errors in mathematical
models, “curse of dimensionality,” model initialization, as well as inadequate
observations and computational errors. These global models represent the
atmosphere on computational grids with horizontal spacing of perhaps 100 km,
stacked in layers, which can be over a kilometer thick. These multi-layer
models are adequate to capture the dynamics of large-scale atmospheric
phenomena. In practice, low-resolution models resolve the large-scale
processes while the effects of the small-scale processes are often
parameterized in terms of the large-scale variables. In addition to higher
degrees of freedom, such large-scale models have the capacity to include the
sub-grid scale dynamics and additive noise. In this case, learning methods are
used to develop parameterization schemes that utilize the explicitly simulated
small-scale processes.

In 1995, Ed Lorenz came up with rather “crude, but a simple model” that
captured the essence of this problem with two stacked layers, which is
referred to as the Lorenz-96 Lorenz1996 (Lorenz, 1995) model. Even though this
simple model is still a long way from the Primitive Equations that we need to
eventually study, it is hoped that the parameterization schemes and filtering
methods presented here can be ultimately adapted for the realistic models that
need to be computationally solved for weather and climate predictions. We use
the Lorenz-96 model to demonstrate the sparse optimization methods and
regression techniques to build models from data. This model, as discussed in
Section 2, has nonlinearities resembling the advective nonlinearities of fluid
dynamics and a multiscale coupling of slow and fast variables similar to what
is seen in the two-layer atmospheric models. Nonlinearities of the governing
physical processes in atmospheric circulations allow energy transfer between
different scales.

The objective is to develop a data-driven method to approximate the effects of
small-scale processes with simplified functions called parameterizations.
Parameterization schemes are inherently approximate, and the development of
schemes which produce realistic simulations is a central challenge of model
development. Most parameterization schemes depend critically on various
parameters whose values cannot be determined a priori but must instead be
found through trial and error. Machine learning is a growing set of
optimization and regression techniques to build models from data sets. Recent
developments in optimization and sparse representations in high-dimensional
spaces have provided powerful new techniques. The goal of this paper is to
replace the current ad hoc parametrization with a combination of machine
learning and data assimilation, and fine-scale process modelling to produce
realistic stochastic parameterization in global weather and climate models.

This paper deals with the parameterization of the Lorenz-96 model with two
time-scales. We assume that the model that expresses the resolved atmospheric
variables of the known physics (represented by $X$ in Figure 1) or the
resolved modes are well understood. Furthermore, Lorenz in introducing this
multiscale model showed that the small-scale features of the atmosphere act
like small random forcing. The theoretical approaches to resolve the sub-scale
variables are two fold. We demonstrate how sparsity-promoting regression (SPR)
techniques based on compressed sensing enables parameterization to generalize
to chaotic systems due to the ergodic nature of the simulated time series. We
present a compressed sensing and sparse recovery approach for constructing
approximations to the Lorenz-96 model that blends physics-informed candidate
functions with functions that enforce sparse structure. In this case, learning
methods are used to develop a reduced model. The ultimate aim is to adapt the
mathematical methods developed here for the realistic models that can be
solved computationally to examine the turbulent atmospheric circulations. A
complete study of these complex systems with practical impact, involves models
of extremely unstable, chaotic dynamical systems with several million degrees
of freedom. In the past few years, data-driven parameterization of the
Lorenz-96 model using machine learning and deep learning techniques Gagne2020
(Gagne et al., 2020) has shown promising results.

The paper is organized as follows. In Section 2, we will introduce the
40-dimensional Lorenz-96 model, that was originally suggested by Lorenz
(1996), and describe the characterization of complex dynamics including
statistical properties, Lyapunov exponents, marginal probability density
functions (PDFs) and auto-correlation functions (ACFs). Section 3 will present
the details on data-driven reduction methods such as regression and compressed
sensing with sparse optimization. In Section 4, stochastic paramterization
using compressed sensing approach will be illustrated on Lorenz-96 model.
Section 5 will describe a general ensemble filtering algorithm, where
forecasting step is straightforward because the Markov kernel is known, while
the update step is much more sophisticated. This nonlinear filtering method
will be applied to parameterized models. Finally, conclusions and future work
discussion will follow in Section 6. To ensure that our results are
reproducible, we have provided our code in GithubCode (Aydogdu et al., 2021).

##  2 The Lorenz-96 Atmospheric Model

Let us now turn to a simple atmospheric model, which nonetheless exhibits many
of the difficulties arising in realistic models, to gain insight into
predictability and data assimilation. The Lorenz-96 model was originally
introduced in Lorenz1996 (Lorenz, 1995) to mimic multiscale mid-latitude
atmospheric dynamics for an unspecified scalar meteorological quantity at $K$
equidistant grid points along a latitude circle:

$\displaystyle\dot{X}_{k}$
$\displaystyle=-X_{k-1}(X_{k-2}-X_{k+1})-X_{k}+F-\frac{hc}{b}\sum_{j=1}^{J}Z_{j,k}$
(1) $\displaystyle\dot{Z}_{j,k}$
$\displaystyle=-cbZ_{j+1,k}(Z_{j+2,k}-Z_{j-1,k})-cZ_{j,k}+\frac{hc}{b}X_{k}$
(2)

Equation (1) describes the dynamics of some atmospheric quantity $X$, and
$X_{k}$ can represent the value of this variable at time $t$ in the $k^{\rm
th}$ sector defined over a latitude circle in the mid-latitude region. (Note
that we use subscripts $k$ and $j$ to conform with the typical spatial
indexing notation used for the Lorenz-96 model. In sections that follow,
subscripts $k$ and $j$ will be used as discrete time indices, not to be
confused with the spatial indices of the Lorenz model). The latitude circle is
divided into $K$ sectors.
For convenience, the vector field associated with the slow components or the
resolved variables is assumed a priori known with representable physics, which
is denoted by

$\mathcal{G}_{k}(X):=-X_{k-1}(X_{k-2}-X_{k+1})-X_{k}+F$ (3)

Figure 1 illustrates the behavior of a generic slow state $X_{k}$ (shown in
the inner ring with $K=8$ for illustration), the fast states in the $1$st
sector, that is, $Z_{1,1},\ldots,Z_{J,1}$ (shown in the outer ring with $J=7$
for illustration), and the fast scale forcing that enters (1); for the $X_{k}$
component the fast scale forcing is

$U_{k}=-\frac{hc}{b}\sum_{j=1}^{J}Z_{j,k}$ (4)

Each $X_{k}$ is coupled to its neighbors $X_{k+1}$, $X_{k-1}$, and $X_{k-2}$
by (1). $X_{k}$ and $Z_{j,k}$ are periodic on the spatial domain, and (1) and
(2) applies for all values of $k$ and $j$ by letting $X_{k+K}=X_{k}$,
$Z_{j,k+K}=Z_{j,k}$, $Z_{j+J,k}=Z_{j,k+1}$ so, for example, in Figure 1 we
have $X_{9}=X_{1}$, $Z_{1,9}=Z_{1,1}$, and $Z_{8,1}=Z_{1,2}$.

Figure 1: A recreation of a similar image by Wilks (2005), that illustrates
the two timescale Lorenz-96 model.

Even though this model is not a simplification of any atmospheric system it is
designed to satisfy three basic properties: contains a quadratic discrete
advection-like term that conserves the total energy (i.e. it does not
contribute to $dE/dt$) just like many atmospheric models; it has linear
dissipation (the –$X_{k}$ term) that decreases the energy defined as
$E=\sum_{k=1}^{K}X_{k}$, and an external forcing term $F>0$ that can increase
or decrease the total energy.

The original $K$ dimensional system was extended to study the influence of
multiple spatio-temporal scales on the predictability of atmospheric flow by
the dividing each segment $k$ into $J$ subsectors ($J=10$), and introducing a
fast variable, $Z_{j,k}$ given by (2), associated with each subsector. Thus,
each $X_{k}$ represents a slowly-varying, large amplitude atmospheric
quantity, with $J$ fast-varying, low amplitude, similarly coupled quantities,
$Z_{j,k}$, associated with it. In the context of climate modeling, the slow
component is also known as the resolved climate modes while the rapidly-
varying component is known as the unresolved non-climate modes. The coupling
terms between neighbors model advection between sectors and subsectors, while
the coupling between each sector and its subsectors model damping within the
model. The model is subjected to linear external forcing, $F$, on the slow
timescale.

We consider $40$ dimensional slow variable components and $400$ dimensional
fast variable components. With $K=40$, one can think of $\mathbf{X}$ as
representing an unspecified scalar meteorological quantity equally spaced
along a mid-latitude belt of roughly 30,000 km, thereby giving the model the
correct synodic scale to mimic Rossby waves. Our main research interest in
this paper is to develop a dimensional reduction and parametrization method
based on compressed sensing Cande2006 (Candès et al., 2006, Donoho, 2006). We
compare this sparse regression scheme based on compressed sensing with Wilk’s
polynomial regression Wilks2005 (Wilks, 2005).

###  2.1 Characterization of complex dynamics

According to Devaney1980 (Devaney, 1989), a chaotic dynamical system defined
on a metric space $X$ by a continuous map $f$ is said to be chaotic on $X$ if
it exhibits the following three essential features; (a) $f$ is transitive, (b)
the periodic points of $f$ are dense in $X$, and (c) $f$ has sensitive
dependence on initial conditions. These systems are completely deterministic
in that if the current state of the system $X_{0}$ is known, all future states
at time $t$, denoted by $X_{t},t>0$, are known. However, even small numerical
errors in the computation of $X_{0}$ can lead to very large differences in the
values of $X_{t}$, so while these systems are deterministic, they are not
predictable. Hence, chaotic systems exhibit a large departure in the long-term
evolution for very tiny changes in initial conditions, and hence to study the
limiting behavior of these systems, one cannot study the global evolution
“orbit by orbit”. Instead, one studies the statistical properties of these
systems, as we describe in this section.

####  2.1.1 Statistical properties of chaotic systems

Although individual orbits of chaotic dynamical systems are by definition
unpredictable, the average behavior of typical trajectories can often be given
a precise statistical description. Anosov (1969), Pesin (2004) and Sinai
(1977) showed that a wide class of “chaotic” maps exhibit natural ergodic
invariant measures, the maps being “chaotic” in the sense that the orbits of
two points, which may be arbitrarily close to each other with respect to a
metric, may evolve very differently over time (positive Lyapunov exponents as
discussed below). Loosely speaking, ergodicity is the “law of large numbers”
for dynamical systems.While the ergodic theorems prove convergence to the
mean, natural questions arise about the deviation of these processes from the
mean and about other statistical properties such as rates of mixing. For a
given invariant measure, and a class of observables, the correlation functions
tell whether (and how fast) the system “mixes”, i.e. “forgets” its initial
conditions.

Many statistical properties, such as the Central Limit Theorem and the Law of
Iterated Logarithm, displayed by independent identically distributed
stochastic processes, have been shown to hold for stochastic processes
generated by chaotic dynamical systems. One of the better known techniques for
establishing such results was proved by Young (1998), where it was shown that
for such (chaotic) dynamical systems with certain “towers” property (a
powerful tool for coding the dynamics of a system in such a way so as to
enable the computation of its statistical properties), behave as Markov chains
on spaces with infinitely many states. Should such dynamical systems mix
sufficiently quickly, the process satisfies the central limit theorem and the
completely deterministic process behaves, at least in the first two moments,
as an independent stochastic process. Furthermore, the degree of independence
is characterized by the rate at which the process mixes.

####  2.1.2 Lyapunov Exponents (LEs)

Among the features listed above, sensitivity to initial conditions is widely
known to be associated with chaotic dynamical systems, and a quantitative
measure of this sensitivity is given by the ${\it Lyapunov\;exponents}$ (LEs)
which are the most commonly used quantities to characterize chaotic phenomena.
The LEs are asymptotic measures characterizing the average rate of divergence
(or convergence) of small perturbations to the solutions of a dynamical
system. The value of the maximum LE is an indicator of the chaotic or regular
nature of orbits.

Considering the phase space, which is assumed to be $\mathbb{R}^{d}$ or a
Riemannian manifold, and denoting by $M$ throughout, the Lebesgue or the
Riemannian measure on $M$ is denoted by $m$, and the dynamics are generated by
iterating a self-map of $M$, written $f:M\circlearrowleft$. Given a
differentiable map $f:M\circlearrowleft$, a point $x\in M$ and a tangent
vector $v$ at $x$, we define

$\lambda(x,v)=\lim_{n\rightarrow\infty}\frac{1}{n}\log|Df_{x}^{n}(v)|\hskip
14.22636pt\text{if this limit exists}$ (5)

and if it does not, then the limit is replaced by $\liminf$ or $\limsup$ and
we write $\underline{\lambda}(x,v)$ and $\overline{\lambda}(x,v)$,
respectively. Thus, $\lambda(x,v)>0$ implies the exponential growth of
$|Df_{x}^{n}(v)|$ and can be interpreted to mean exponential divergence of
nearby orbits. The computational procedure for Lyapunov exponents using the
standard method can be found in Wolf1985 (Wolf et al., 1985, Skokos, 2010).

When the chaotic system is dissipative, solutions will eventually reside on an
attractor. The dimension of this attractor, i.e., the number of “active”
degrees of freedom, is also known as the Kaplan-Yorke dimension $D_{KY}$
Eckmann1985 (Eckmann et al., 1985) and is given as

$D_{KY}=r+\frac{1}{|\tilde{\lambda}_{r+1}|}\sum_{i=1}^{r}\tilde{\lambda}_{i}$
(6)

where
$\tilde{\lambda}_{1}\geq\tilde{\lambda}_{2}\geq\cdots\geq\tilde{\lambda}_{d}$
are the Lyapunov exponents counted with multiplicity ($d=\mbox{dim}(M)$) and
$r$ is the largest integer such that $\sum_{i=1}^{r}\tilde{\lambda}_{i}>0$.

The Kolmogorov-Sinai entropy $H_{KS}$ is a measure of the diversity of the
trajectories generated by the dynamical system. This entropy measure is
determined through the Pesin formula, which provides an upper bound to
$H_{KS}$ Eckmann1985 (Eckmann et al., 1985) as

$H_{KS}=\sum_{\tilde{\lambda}_{i}>0}\tilde{\lambda}_{i}.$ (7)

Complexity of the dynamics of the Lorenz-96 model varies considerably with
different values of the constant forcing term $F$. It is shown in Majda2005
(Majda et al., 2005), the model is in weakly chaotic regimes for $F=5,\,6$,
strongly chaotic regime for the values of $F$ from $8$ to $10$, and turbulent
regimes for $F=12,\,16,\,24$. We will present our results for a specific value
of forcing $F=10$.

In our numerical simulations, we begin from random initial conditions and use
a fourth-order Runge-Kutta time integration with a time step of 0.01.
Typically, we integrate forward for 500 time units before starting the
calculation of Lyapunov exponents to ensure that all transients have decayed.
At this point we begin integrating the tangent space equations using the Gram-
Schmidt procedure update every $\tau=0.2$ time units. With $K=40$, $F=10$,
$h=0$ (just the slow variables), $c=10$ and $b=10$, the summary of Lyapunov
Exponents obtained from the above simulation of the Lorenz-96 system is given
Table 1. The time evolution of the Lyapunov exponents for the Lorenz-96 model
is shown in Figure 2. As noted before,the calculation of the Lyapunov
exponents begins at 500 time units and Figure 2 depicts the convergence of the
Lyapunov exponents around 1000 time units.

Largest Lyapunov exponent $\lambda_{1}$ |  2.3098
---|---
Error doubling time |  0.3 time units
Number of strictly positive Lyapunov exponents |  14
Number of neutral, $\lambda\in[-1E-2,1E-2]$, exponents |  1
Number of strictly negative Lyapunov exponents |  26
Kaplan-Yorke dimension |  29.4694
Kolmogorov-entropy |  14.8409
Table 1: Summary of Lyapunov Exponents Figure 2: Time evolution of Lyapunov
exponents for the Lorenz-96 system

####  2.1.3 Marginal Probability Density Functions (PDFs) and Auto-correlation
Functions (ACFs)

Typically, the PDFs are computed using bin-counting also known as histogram
methods. For estimating the PDF of a continuous variable $x$, standard
histogram methods simply partition $x$ into distinct bins of width
$\Delta_{i}$, and then the number of observations $n_{i}$ of $x$ falling into
bin $i$ is converted into a normalized probability density for that bin.

More smoother estimates for PDFs can be obtained by using the Gaussian product
kernel function which results in the estimation of the PDF known as the kernel
density estimation (KDE) Bishop2006 (Bishop, 2006):

$p_{X}(x)=\frac{1}{N}\sum_{j=1}^{N}\frac{1}{(2\pi\lambda^{2})^{1/2}}\exp\\{-\frac{\|x_{j}-x\|^{2}}{2\lambda^{2}}\\}$
(8)

where $\lambda$ represents the standard deviation of the Gaussian components
and $N$ is the total number of data points. Figure 3 depicts the marginal PDFs
of some (randomly selected) slow variables $X_{1}$, $X_{8}$, $X_{19}$ and
$X_{31}$ simulated as described above in Section 2. The convergence of the
marginal PDFs shown in Figure 3 demonstrates the existing symmetry among the
slow variables.

Figure 3: PDFs of the slow variables $X_{1}$, $X_{8}$, $X_{19}$ and $X_{31}$

The strongly mixing character of the slow dynamics can be inferred from the
correlation function, which means that the correlation function decays quickly
to zero or near zero values. The discrete version of the auto-correlation
functions (ACFs) for the slow variables is given by

$\displaystyle C_{k,k}(m\Delta t)=$
$\displaystyle\frac{1}{M-m}\sum_{h=1}^{M-m}(X_{k}(h\Delta
t)-\bar{X})(X_{k}((h+m)\Delta t)-\bar{X})$ (9)

where

$\bar{X}=\frac{1}{M}\sum_{m=1}^{M}X_{k}(m\Delta t)$ (10)

and these formulae can be given similarly for the fast variables. Figure 4
shows the ACF for the slow variable $X_{1}$ which decays very quickly and then
oscillates around zero illustrating the strongly mixing character of the
dynamics.

Figure 4: ACF of the slow variable $X_{1}$

##  3 Data-driven Reduction Methods

The objective of this section is to compare model-based and data-driven tools
for dimensional reduction and reducing uncertainties in prediction. It is
impossible to effectively “learn” from high dimensional data unless there is
some kind of implicit or explicit low-dimensional structure – this can be made
mathematically precise in different ways. The methods presented below provide
the low-dimensional structure either because of the presence of multi-scales
in the model or due to the measurement matrix chosen to satisfy with high
probability the so called “Restricted Isometry Property (RIP)” condition
introduced in compressed sensing Cande2006 (Candès et al., 2006).

Compressed sensing theory addresses the problem of recovering high-dimensional
but sparse vectors from a limited number of measurements - the number of
observed measurements $m$ is significantly smaller than the dimension $n$ of
the original vector. We present a compressed sensing approach for sparse
regression problem and using the ideas from concentration of measure, extend
it to a case where the samples are generated from a chaotic time series, which
are sufficiently mixing, thus ergodic, as shown in subsection 2.1. Here, the
sparse recovery is achieved by constructing a matrix $\Theta$ in section 3.2
from ergodic samples generated by the Lorenz-96, such that it is suitable for
reconstructing sparse parameter vectors via $l_{1}$-minimization.

###  3.1 Regression

This parametrization was proposed by Wilks Wilks2005 (Wilks, 2005), and it
uses a polynomial regression to parametrize the effects (4) of unresolved
variables:

$U_{k}(t)=P(X_{k}(t))+e_{k}(t)$ (11)

where $P$ is a 4th degree polynomial in $X_{k}$ and is given as:

$P(X_{k})=a_{0}+a_{1}X_{k}+a_{2}X_{k}^{2}+a_{3}X_{k}^{3}+a_{4}X_{k}^{4}$ (12)

And $e_{k}$ represents noise, which can be written as a first-order
autoregressive model:

$e_{k}(t_{i})=\varphi e_{k}(t_{i-1})+\sigma z_{k,i}$ (13)

where $z_{k,i}\stackrel{{\scriptstyle i.i.d}}{{\sim}}N(0,1)$ for all $i,k$.
The method of calculating $\varphi$ and $\sigma$ is explained in Appendix A.

###  3.2 Compressed Sensing and Sparse Recovery

Over the past two decades, researchers have focused on sparsity as one type of
low-dimensional structure. Given the recent advances in both compressed
sensing Cande2006 (Candès et al., 2006, Donoho, 2006, Candès et al., 2008) and
sparse regression Tibshirani1996 (Tibshirani, 1996) it has become
computationally feasible to extract system dynamics from large, multimodal
datasets. The term sparse in signal processing context refers to the case
where signals (or any type of data, in general) have merely few non-zero
components with respect to the total number of components. Sparsity plays a
key role in optimization and data sciences. In the context of our work, these
techniques rely heavily on the fact that many dynamical systems can be
represented by governing equations that are sparse in the space of all
possible functions of a given algebraic structure.

Compressed sensing (CS) is a technique for sampling and reconstructing sparse
signals, that is, the signals have few large components. More precisely,
$k$-sparse signals are those that can be represented by $k<<n$ significant
coefficients over an $n$\- dimensional basis. The central goal of CS is to
capture attributes of a signal using very few measurements, which
distinguishes CS from other dimensionality reduction techniques - a
methodology for the recovery of sparse vectors from a small number of linear
measurements. Hence, this allows for polynomial-time reconstruction of the
sparse signal Donoho2006 (Donoho, 2006).

In Donoho (2006) and Cande‘s and Tao (2006), the original sparse ($k$-sparse)
signal is projected onto a lower-dimensional subspace via a random projection
scheme, called the sampling matrix. More precisely, this broader objective is
exemplified by the important special case in which one is interested in
finding a vector $X\in\mathbb{R}^{n}$ using the (noisy) observation or the
measurement data

$Y=\Theta X+\eta,\quad\text{where}\quad\Theta\in\mathbb{C}^{m\times
n}\;\;\text{with}\;\;k<m<n,$ (14)

is the known sensing or sampling matrix and $\eta$ is the measurement noise.

In general, the problem cannot be solved uniquely. However, if $X$ is
$k$-sparse i.e., if it has up to $k$ nonzero entries, the theory of CS shows
that it is possible to reconstruct $X$, a $k$-sparse vector in
$X\in\mathbb{R}^{n}$ uniquely, from $m$ linear measurements even when $m<<n$,
by exploiting the sparsity of $X$. This can be achieved by finding the
sparsest signal consistent with the vector of measurements Donoho2006 (Donoho,
2006), i.e.

$\arg\;\min_{x\in\mathbb{R}^{n}}\|x\|_{0}\quad\text{subject
to}\quad\|y-\Theta\,x\|_{{}_{2}}\leq\varepsilon,$ (15)

where $\|x\|_{0}$ denotes the $l_{0}$ norm for $x$ (i.e., the number of non-
zero entries of $x$), while $\varepsilon$ denotes a parameter that depends on
the level of measurement noise $\eta$. It can be shown that the $l_{0}$
minimization method can exactly reconstruct the original signal in the absence
of noise using a properly chosen sensing matrix $\Theta$ whenever $m>2k$.
However, $l_{0}$ minimization problem (15) is non-convex and combinatorial,
which is NP-hard.

Hence, instead of problem (15) we consider its $l_{1}$ convex relaxation which
may be stated as Chen1998 (Chen et al., 1998)

$\arg\;\min_{x\in\mathbb{R}^{n}}\|x\|_{1}\quad\text{subject
to}\quad\|y-\Theta\,x\|_{{}_{2}}\leq\varepsilon,$ (16)

where $\|\cdot\|_{1}$ represents the $l_{1}$ norm, which is a convex function
and problem (16) is a convex optimization problem which can accurately
approximate the sparse signal $x$ (solution to problem (15)) in polynomial
time if measurement matrix $\Theta$ is chosen to satisfy the necessary RIP
condition with high probability Cande2005 (Candès et al., 2005, Candès et al.,
2006). Loosely speaking, if $\Theta$ satisfies the RIP, then the measurement
matrix approximately preserves the Euclidean length of every $k$-sparse
signal. Equivalently, all subsets of $k$ columns taken from $\Theta$ are
nearly orthogonal. One should note that the $l_{1}$ minimization in (16) is
closely related to the Lasso problem Tibshirani1996 (Tibshirani, 1996).

$\arg\;\min_{x\in\mathbb{R}^{n}}\frac{1}{2}\|y-\Theta\,x\|_{{}_{2}}^{2}+\lambda\,\|x\|_{1}$
(17)

where $\lambda\geq 0$ is a regularization parameter. If $\varepsilon$ and
$\lambda$ in (16) and (17) satisfy some special conditions, the two problems
are equivalent; however, characterizing the relationships between
$\varepsilon$ and $\lambda$ is difficult, except for the special case of
orthogonal sensing matrices $\Theta$.

Even though there are several objectives in compressed sensing, we are mainly
interested in efficiently reconstructing the $k$-sparse vector from the
measurements. Hence, the RIP plays an important role in our calculations.
Unfortunately, direct verification of RIP for a given matrix is not feasible.
Nevertheless, it has been shown that there are some matrices, particularly
random matrices with independent and identically distributed (i.i.d) entries,
satisfy RIP with high probability.

We now enforce this sparse structure within learning. One of the related
problems in the context of sparsity is “dictionary learning”. We refer to the
columns of the ${m\times n}$ matrix $\Theta$, as a dictionary, and they are
assumed to be a set of vectors capable of providing a highly succinct
representation for the signal. One has freedom to determine the particular
dictionary elements; however, the choice of monomial terms provides a model
for the interactions between each of the components of the time-series and is
used for model identification of dynamical systems in the general setting. In
dictionary learning, which we shall use to determine $\Theta$, the aim is to
find a dictionary that can sparsely represent a given set of training data for
further analysis –“sparsifying matrix” from a given set of training data.
Regularized optimization with polynomial dictionaries is used to approximate
the generating function of some unknown dynamic process. When the dictionary
is large enough so that the “true” function is in the span of the candidate
space, the solutions produced by sparse optimization are guaranteed to be
accurate. This means that any of the training data should be presentable using
linear combinations of few columns from the dictionary. In Lorenz-96, the
sampling matrix or the dictionary $\Theta$ is formed as the collection of
matrix of all monomials (stored column-wise) of resolved modes and due to the
ergodic nature of the simulated time series the entries satisfy RIP with high
probability.

In this paper, the effects of unresolved processes (2) on the resolved modes
(1) in the Lorenz-96 model are denoted by
$u_{k}=-\frac{hc}{b}\sum_{j=1}^{J}Z_{j,k}$. The CS method takes snapshot data
${\bf x}(t)\in\mathbb{R}^{n}$ and attempts to discover the structure of a
potentially highly-nonlinear mapping of the form

${\bf u}(t)={\bf f}({\bf x}(t))$ (18)

where the $n$-dimensional vector ${\bf x}(t)=[x_{1}(t)\ x_{2}(t)\ \ldots\
x_{n}]^{T}$ represents the state of the system at time $t$. Basis of this
construction lies on the fact that in the multi scale Lorenz-96 model, the
unresolved dynamics depends on the value of the resolved variables $\bf x$.

In particular, this method searches for a sparse approximation of the function
mapping ${\bf f}({\bf x})=[f_{1}({\bf x})\ f_{2}({\bf x})\ \ldots\ f_{n}({\bf
x})]^{T}\in\mathbb{R}^{n}$ using a dictionary of basis functions yielding

$f_{k}({\bf x})\approx\sum_{j=1}^{p}\theta_{j}({\bf x})s_{jk}={\bf\Theta}({\bf
x}){\bf s}_{k}$ (19)

where ${\bf\Theta}({\bf x})=[\theta_{1}({\bf x}),\theta_{2}({\bf
x}),...,\theta_{p}({\bf x})]\in\mathbb{R}^{1\times p}$ form the dictionary of
basis functions, and ${\bf s}_{k}=[s_{1k}\ s_{2k}\ \ldots\
s_{pk}]^{T}\in\mathbb{R}^{p}$ is the vector of coefficients. Sparsity is
expressed “not” in terms of an orthonormal basis but in terms of an over-
complete dictionary. The majority of the coefficients $s_{jk}$ are zero while
the remaining nonzero entries identify the active terms contributing to the
sparse representation of ${\bf f}({\bf x})$.

To learn the function mapping $\bf{f}$ from data, time-histories of vectors
${\bf x}(t)$ and ${\bf u}(t)$ are collected by performing simulations for the
Lorenz-96 system. The input data matrix ${\bf X}\in\mathbb{R}^{m\times n}$ is
made of $m$ snapshots of the state vector sampled at several times
$t_{1},t_{2},...,t_{m}$ where subscripts $1,2,\ldots,n$ indicates the elements
of the state vector. The output data matrix ${\bf U}\in\mathbb{R}^{m\times n}$
follows the same arrangement, and these two data matrices are given as
follows:

$\displaystyle\bf{X}=\left[\begin{array}[]{c}{\bf x}^{T}(t_{1})\\\ {\bf
x}^{T}(t_{2})\\\ \vdots\\\ {\bf
x}^{T}(t_{m})\end{array}\right]=\left[\begin{array}[]{cccc}x_{1}(t_{1})&x_{2}(t_{1})&\ldots&x_{n}(t_{1})\\\
x_{1}(t_{2})&x_{2}(t_{2})&\ldots&x_{n}(t_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
x_{1}(t_{m})&x_{2}(t_{m})&\ldots&x_{n}(t_{m})\end{array}\right],$
$\displaystyle{\bf U}=\left[\begin{array}[]{c}{\bf u}^{T}(t_{1})\\\ {\bf
u}^{T}(t_{2})\\\ \vdots\\\ {\bf
u}^{T}(t_{m})\end{array}\right]=\left[\begin{array}[]{cccc}u_{1}(t_{1})&u_{2}(t_{1})&\ldots&u_{n}(t_{1})\\\
u_{1}(t_{2})&u_{2}(t_{2})&\ldots&u_{n}(t_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
u_{1}(t_{m})&u_{2}(t_{m})&\ldots&u_{n}(t_{m})\end{array}\right].$

####  3.2.1 Dictionary of Basis Functions

First, we construct an augmented dictionary matrix ${\bf\Theta}({\bf
X})\in\mathbb{R}^{m\times p}$ where each column represents a basis function
that can be a potential candidate for the terms in ${\bf f}$ to be discovered
from data. $p$ is the maximal number of $n$-multivariate monomials of degree
at most $d$. As stated in Brunton2016 (Brunton et al., 2016), there is huge
flexibility in choosing the basis functions to populate the dictionary matrix
${\bf\Theta}({\bf X})$. For example, the dictionary matrix ${\bf\Theta}({\bf
X})$ may consist of constant, polynomial and trigonometric terms as shown in
(25).

$\displaystyle{\bf\Theta}({\bf
X})=\left[\begin{array}[]{ccccccccc}|&|&|&\ldots&|\\\
\bf{1}&\bf{X}&\bf{X}^{P_{2}}&\ldots&\bf{X}^{P_{d}}\\\
|&|&|&\ldots&|\end{array}\right]_{m\times p}$ (25)

The dictionary matrix ${\bf\Theta}({\bf X})$ is constructed by stacking
together, column wise, candidate nonlinear functions of ${\bf X}$. Here,
higher order polynomials are denoted as ${\bf X}^{P_{d}}$ where $d$ is the
order of the polynomial considered. For example, element 1 is a column-vector
of ones, element ${\bf X}$ is as defined above, element ${\bf X}^{P_{2}}$ is
the matrix containing the set of all quadratic polynomial functions of the
state vector ${\bf x}$, and is constructed as shown in (30).

$\displaystyle\bf{X}^{P_{2}}=\left[\begin{array}[]{cccccc}x_{1}^{2}(t_{1})&x_{1}(t_{1})x_{2}(t_{1})&\ldots&x_{2}^{2}(t_{1})&\ldots&x_{n}^{2}(t_{1})\\\
x_{1}^{2}(t_{2})&x_{1}(t_{2})x_{2}(t_{2})&\ldots&x_{2}^{2}(t_{2})&\ldots&x_{n}^{2}(t_{2})\\\
\vdots&\vdots&\ddots&\vdots&\ddots&\vdots\\\
x_{1}^{2}(t_{m})&x_{1}(t_{m})x_{2}(t_{m})&\ldots&x_{2}^{2}(t_{m})&\ldots&x_{n}^{2}(t_{m})\end{array}\right]$
(30)

####  3.2.2 Sparse Optimization

Sparse regression is performed to approximately solve

${\bf U}\approx{\bf\Theta}({\bf X}){\bf S}$ (31)

to determine the vector coefficients ${\bf S}=[{\bf s}_{1}\ {\bf s}_{2}\
\ldots\ {\bf s}_{n}]\in\mathbb{R}^{p\times n}$ that determines the active
terms in ${\bf f}$. Since only a few of the candidate functions are expected
to have an effect on the function mapping ${\bf f}$, all coefficient vectors
${\bf s}_{i}$ are expected to be sparse. In this work, we use the Lasso
algorithm Tibshirani1996 (Tibshirani, 1996) for solving the above sparse
regression problem.

The goal of the sparse optimization to solve the regression problem (31) can
be defined using the Lasso form (17) for $k=1,...,n$

${\bf s}_{k}=\underset{{\bf s}_{k}^{{}^{\prime}}}{argmin}\;||{\bf
u}_{k}-{\bf\Theta}({\bf X}){\bf s}_{k}^{{}^{\prime}}||_{2}^{2}+\lambda||{\bf
s}_{k}^{{}^{\prime}}||_{1}$ (32)

which uses $l_{1}$ norm constraint on the coefficient vector ${\bf s}_{k}$.
The term $\lambda||{\bf s}_{k}^{{}^{\prime}}||_{1}$ is the regularization
term, which penalizes coefficients different from 0 in a linear fashion. It is
the actual promoter of sparsity in the minimization problem.

In the Lasso method, given a collection of $m$ time snapshot pairs $({\bf
x}(t_{i}),{\bf u}(t_{i}))_{i=1}^{m}$, we estimate the coefficients in (31) by
solving the following optimization problem for $k=1,\ldots,n$

$\underset{{\bf
s}_{k}}{minimize}\;\sum_{i=1}^{m}(u_{k}(t_{i})-\sum_{j=1}^{p}\theta_{j}({\bf
x}(t_{i}))s_{jk})^{2}+\lambda\sum_{j=1}^{p}|s_{jk}|$ (33)

for some Lagrange multiplier $\lambda\geq 0$. The Lagrangian multiplier
$\lambda$ can be specified using a cross-validation procedure.

##  4 Stochastic Parametrization: Application to Lorenz-96

Having introduced the data-driven methods for model reduction in the presence
of time-scale separation, let us now turn to a simple atmospheric model, which
nonetheless exhibits many of the difficulties arising in realistic models, to
gain insight into predictability and data assimilation. The dynamics of the
Lorenz-96 equation are shown in (1) and (2). Let $X_{k}$ and $Z_{j,k}$ be the
numerical solutions of the Lorenz-96 equations, generated using the 4th order
Runge-Kutta method, for $k=1,...,K$, $j=1,...,J$. We define $U_{k}$ as in
Equation (4):

$U_{k}=-\frac{hc}{b}\sum_{j=1}^{J}Z_{j,k}$

The resolved modes in (1) is now modified as:

$\dot{X}_{k}=-X_{k-1}(X_{k-2}-X_{k+1})-X_{k}+F+U_{k}=\mathcal{G}_{k}(X)+U_{k}$
(34)

Our goal is to use the methods discussed in Section 3 to form
parametrizations, in other words, express $U_{k}$ as a function of
$X_{1},...,X_{K}$ and compare each of them. Based on the ergodic properties of
the Lorenz-96 system shown in Section 2.1, the Lyapunov exponents stabilize at
$t=500$. The training set is set to $t_{train}=[500,500+h,500+2h,...,1000]$
where $h=0.01$. The test set is set to
$t_{test}=[1000,1000+h,1000+2h,...,2000]$.

Using polynomial regression to solve stochastic parametrizations in the
Lorenz-96 system was proposed by Wilks (2005). The goal of this method is to
find a degree-4 polynomial $P$ such that $U_{k}$ can be approximated by
$P(X_{k})$ for $k=1,...,K$. This requires solving a least squares problem:

$\sum_{k=1}^{K}\sum_{i\in t_{train}}(U_{k}(t_{i})-P(X_{k}(t_{i})))^{2}$ (35)

The polynomial obtained from this method is:

$\displaystyle U_{k}$ $\displaystyle=-0.0003142X_{k}^{4}+0.004882X_{k}^{3}$
$\displaystyle+0.005197X_{k}^{2}-0.47362X_{k}-0.03779$ (36)

The goal of compressed sensing in this context is to find a sparse matrix $S$
that solves for $Y=\Theta S$. Let $Y$ be the matrix with columns
$U_{1},...,U_{K}$, each divided by its $L_{2}$ norm. Every column of $\Theta$
is a basis function of $X_{1},...,X_{k}$ divided by its $L_{2}$ norm. The goal
of compressed sensing is to find unique equations for every $U_{k}$ for
$k=1,...,K$. Since regression gave us a polynomial of degree 4, our goal is to
express $U_{k}$ as a function of $X$ values of order 4. However, with $K=40$,
this means the $\Theta$ matrix will have $135750$ columns, and thus,
compressed sensing will be computationally expensive. Each of the
parametrizations we get are shown in the ”CompressedSensingEquations” folder
in GithubCode (Aydogdu et al., 2021).

A computationally less expensive approach will be to find a solution of order
2. If we average the coefficients of the derived parametrization, we get:

$U_{k}=-0.376X_{k}+0.0166X_{k}^{2}$ (37)

This parametrization shows that every $U_{k}$ is expressed as a function of
$X_{k}$ and $X_{k}^{2}$, thus showing that every basis function $X_{i}X_{j}$
where $i\neq j$ is not relevant, and the $\Theta$ matrix can be reduced to a
matrix with $X_{k}$ and $X_{k}^{2}$ as the columns.

We then explore solutions of order 3. Given $K=40$, $X^{P_{3}}$ has $11480$
columns, thus finding a solution of order 3 is computationally expensive.
However, creating multiple parametrizations, each with a subset of columns of
$X^{P_{3}}$ is computationally less expensive. Thus, in this section, we
create 100 parametrizations, each of which the $\Theta$ matrix contains
$X_{k}^{3}$ and 100 random columns from $X^{P_{3}}$. If we average the
coefficients of the model, we get:

$U_{k}=-0.0338X_{k}+0.00452X_{k}^{2}+0.00142X_{k}^{3}$ (38)

This shows that every $U_{k}$ is expressed as a function of $X_{k}$,
$X_{k}^{2}$ and $X_{k}^{3}$, thus showing that all remaining columns in
$X^{P_{3}}$ are not relevant. Based on the results of the previous models, we
expect $X_{1}^{4},...,X_{K}^{4}$ to be the only relevant columns from
$X^{P_{4}}$. Thus, our $\Theta$ matrix contains the basis functions
$1,X_{k},X_{k}^{2},X_{k}^{3},X_{k}^{4}$ for $k=1,...,K$. We will call this
parametrization the raw degree-4 compressed sensing model. Since this
algorithm will be used for testing purposes, the compressed sensing algorithm
also calculates biases for each $U_{k}$. If we average the coefficients, we
get:

$\displaystyle U_{k}$ $\displaystyle=-0.416X_{k}+0.00277X_{k}^{2}$
$\displaystyle+0.00509X_{k}^{3}-0.000306X_{k}^{4}-0.0243$ (39)

The biases are calculated by the compressed sensing algorithm as
$[2.173746,1.393404,...,-2.973427]$ and their mean is $-0.0243$. The problem
with this parametrization is that there is too much variation in the biases.

We investigate other ways to modify the biases. In the following model, the
biases are strictly set to zero while keeping other coefficients constant. If
we average the coefficients, we get:

$\displaystyle U_{k}$ $\displaystyle=-0.474X_{k}+0.00277X_{k}^{2}$
$\displaystyle+0.00509X_{k}^{3}-0.000306X_{k}^{4}$ (40)

In the next model, the biases are strictly set to the average of the biases
while keeping other coefficients constant. If we average the coefficients, we
get the same equation as (4). Both of these models perform equally well in
predicting the trajectory of the system and significantly better than the raw
4-dimensional model. The model where biases are strictly set to the average of
the biases performs better in the Kullback-Leibler divergence measure, but
does not outperform linear regression.

A small variability was introduced in the biases with the goal of
outperforming linear regression in terms of the Kullback-Leibler divergence
measure. The biases were selected by sampling from $N(\mu,\sigma^{2})$, where
$\mu$ is the average of the biases ($\mu=-0.0243$). Through multiple samples,
$\sigma=0.07$ gave the best results. This leads to a new model, which
outperforms all other models. The model is shown in Appendix B. If we average
the coefficients, we get:

$\displaystyle U_{k}$ $\displaystyle=-0.474X_{k}+0.00277X_{k}^{2}$
$\displaystyle+0.00509X_{k}^{3}-0.000306X_{k}^{4}-0.0241$ (41)

We will use (41) as our compressed sensing model in the later sections.

The equations of the resolved variables associated with the vector field
denoted by $\mathcal{G}_{k}(X)$ are equivariant with respect to a cyclic
permutation of the variables. Therefore, the system with dimension $K$ is
completely determined by the equation for the $k^{\rm th}$ variable. Even
though the unresolved dynamics depend on the value of the resolved variables
$\bf X$, while computing the parametrization of $U_{k}$ by CS, there were no
restrictions placed on the form of the polynomial on the sparse regression.
The consequences of ${\bf\Theta}({\bf X})$ satisfying RIP (hence, sparsity)
are quite profound as shown in results provided in Appendix B. Compared to
regression, compressed sensing also has the advantage that each unresolved
variable $U_{k}$ is given a unique equation that best represents its dynamics
as a function of resolved variables.

###  4.1 Comparison of Trajectories

Section 2.1.2 shows that the largest Lyapynov exponent is 2.3. This shows that
the error-doubling time is $\frac{\ln(2)}{2*2.3}\approx 0.3$. To show the
practicality of each method applied to Lorenz-96, we predict the trajectory of
each $X_{k}$ for a time interval of three times the error doubling time (0.9)
after the end of the training set. The measure of inaccuracy used in this
paper is the Mean Squared Prediction Error (MSPE) MSPE1991 (Pindyck et al.,
1991).

$MSPE=\frac{1}{K}\sum_{k=1}^{K}\frac{1}{I}\sum_{i\in
I}(\hat{X}_{k}(t_{i})-X_{k}(t_{i}))^{2}$ (42)

where $X_{k}$ is derived from (1) and (2) and $\hat{X}_{k}$ is the
approximation to $X_{k}$ derived from a parametrized model (34) through the
4th order Runge-Kutta method, $I$ is the time interval where we compare the
trajectories. In this case, $I=[1000,1000+h,1000+2h,...,1000+0.9]$ where
$h=0.01$. The initial conditions are set as:
$\hat{X}_{k}(1000)=X_{k}(1000),k=1,...,K$ The MSPEs derived for each method
are shown in Table 2.

Method | MSPE | Average K-L Divergence
---|---|---
Regression | 0.03606 | 0.06729
Compressed sensing (Raw) | 5.9993 | 1.9891
Compressed sensing [1] (Setting biases to zero) | 0.03398 | 0.10099
Compressed sensing [2] (Setting biases to average bias) | 0.03440 | 0.07351
Compressed sensing [3] (Adding noise to average bias) | 0.03387 | 0.05919
Table 2: Evaluation of each method

###  4.2 Comparison of Probability Density Functions

This section will compare the probability density functions (PDFs) of
$\hat{X}_{k}$ and $X_{k}$ for $k=1,...,K$ during the time interval $t_{test}$.
The PDFs are computed using the kernel density estimation, given in (8). PDFs
are computed for every $X_{k}$ and $\hat{X_{k}}$ using their data in the time
interval $[500,2000]$.

The Kullback-Leibler divergence was introduced by Kullback and Leibler (1951)
and measures the discrepancy between two PDFs. The Kullback-Leibler divergence
is calculated as

$D_{KL}(P||Q)=\sum_{x\in X}P(x)\log(\frac{P(x)}{Q(x)})$ (43)

where $P$ is the PDF calculated from the data of $X_{k}$ and $Q$ is the PDF
calculated from $\hat{X}_{k}$. We average the Kullback-Leibler divergence in
the PDFs of $\hat{X}_{k}$ and $X_{k}$ for all $k$ to compare each model. The
average Kullback-Leibler Divergence derived for each method is shown in Table
2.

###  4.3 Auto-regressive Models in Residuals

So far, we have investigated the deterministic part in stochastic
parametrization. The goal of this subsection is to investigate the stochastic
part in stochastic parametrization to further improve the regression and
compressed sensing models explored so far. This is done by obtaining
$\varphi,\sigma$ and $\sigma_{e}$ from the first-order AR model in (13). The
method of calculating these values is explained in Appendix A. Throughout this
subsection, we will refer to the compressed sensing model with noise added to
the average bias as the compressed sensing model.

The values for $\varphi,\sigma$ and $\sigma_{e}$ that we obtain are shown in
Table 3.

Method | $\varphi$ | $\sigma$ | $\sigma_{e}$ | MSPE | Ave. K-L Divergence
---|---|---|---|---|---
Regression | 0.9453 | 0.2265 | 0.6945 | 0.02624 | 0.07692
Compressed sensing | 0.9981 | 0.1710 | 2.775 | 0.02066 | 0.07142
Table 3: AR evaluation of residuals

Given the values of $\varphi$ and $\sigma$, we can run a numerical simulation
of the reduced dynamical system with AR residuals. The equation is shown
below:

$\dot{X}_{k}=-X_{k-1}(X_{k-2}-X_{k+1})-X_{k}+F-f_{k}(\textbf{X}(t))+e_{k}(t_{i})$
(44)

where $f_{k}(\textbf{X}(t))$ represents either the Wilk’s parametrization (4)
or the CS parametrization (41) and $e_{k}(t_{i})$ is generated using the first
order AR model below:

$\displaystyle e_{k}(t_{i})$ $\displaystyle=\varphi e_{k}(t_{i-1})+\sigma
z_{k}(t_{i})$ $\displaystyle=\varphi
e_{k}(t_{i-1})+\sigma_{e}(1-\varphi^{2})^{1/2}z_{k}(t_{i})$ (45)

$z_{k}(t_{i})$ is sampled randomly from $N(0,1)$.

Numerical simulations are done using the 4th order Runge-Kutta method.
$e_{k}(t_{0})$ is set to $0$ for all $k$. Similar to section 4.2, we evaluate
our model by assessing its trajectory as well as its PDF. The results are
summarized in Table 3.

A comparison of Table 3 with Table 2 shows that using an AR model in
compressed sensing and regression has significantly improved the prediction of
the trajectory of resolved variables as the MSPE is lower. Compressed sensing
shows the lowest MSPE and Kullback-Leibler divergence.

The goal of this section was to create a deterministic and AR model to improve
the prediction of resolved variables. Based on the results shown, compressed
sensing improves the prediction significantly compared to regression. The
coefficients of the deterministic model as well as the values of
$\varphi,\sigma$ and $\sigma_{e}$ will be useful in the following section
where filtering methods will be used to further improve the prediction of
resolved variables.

##  5 Nonlinear Filtering Applied to Parameterized Models

The second component of this paper deals with nonlinear filtering or data
assimilation (DA). The main focus of filtering is to combine computational
models with sensor data (observations) to predict the dynamics of large-scale
evolving systems. Filtering provides a recursive algorithm for estimating a
signal or state of a random dynamical system based on noisy measurements. The
signal that is represented by a Markov process cannot be accessed or observed
directly and is to be “filtered” from the trajectory of the observation
process which is statistically related to the signal. Suppose we make a
forecast about the behavior at a future time of a complex system with some
uncertainties (randomness) and there is near-continuous data available from
remote sensing instrumentation networks. Then, as new information becomes
available through observations, it is natural to ask how to best incorporate
this new information into the prediction. Data assimilation can be simply
defined as the combination of observations with the model output to produce an
optimal estimate of dynamical model.

The standard notation for the discrete DA problems from the time $t_{k}$ to
$t_{k+1}$ can be formulated as follows:

$x^{f}(t_{k+1})=M_{k+1}[x^{f}(t_{k})]+w_{k}$ (46)

where $x$ is the model state, $M$ is the model dynamics and $w$ is the model
error. The operator $M$ is represented with matrix in discrete case and
differential operator in continuous case. Observations at time $t_{k}$ can be
defined by:

$y_{k}=H_{k}[x^{t}(t_{k})]+v_{k}$ (47)

where $y$ is the observation state and $H$ is the observation operator. $v$ is
the white noise with associated covariance matrix describing the observation
error.

Data assimilation methods consist of analysis (or correction/filtering) and
forecast (or prediction) steps. At time $t_{k}$, we have the output state of
the background forecast $x_{k}^{f}$ and observation state $y_{k}$. The (linear
or nonlinear) analysis step is based on the outputs of $x_{k}^{f}$ and $y_{k}$
and produces analysis $x^{a}_{k}$. Then, applying the model dynamics on the
analysis step generates next forecast output as $x_{k+1}^{f}$.

###  5.1 Ensemble Kalman Filter

It is well documented that particle filters Gordon1993 (Gordon et al., 1993)
suffer from the well-known problem of sample degeneracy Bengtsson2008
(Bengtsson et al., 2008, Snyder et al., 2008). In contrast, the ensemble
Kalman filter (EnKF) Evensen1994 (Evensen, 1994, Evensen, 2003) can handle
some of these difficulties, where the dimensions of states and observations
are large and the number of replicates is small. However, EnKF incorrectly
treats the non-Gaussian features of the forecast distribution that arise in
certain nonlinear systems. Nevertheless, we apply ensemble Kalman filter for
data assimilation, because they are easy to implement in large systems.

For the Kalman filter Kalman60 (Kalman, 1960), the modelling and observation
errors are assumed to be independent with Gaussian/normal probability
distributions:

$w_{k}\sim N(0,Q_{k})\quad\text{and}\quad v_{k}\sim N(0,R_{k})$ (48)

where $Q$ and $R$ are the covariance matrices of the model error and
observation error, respectively. The Kalman filter is of the form:

$x_{k}^{a}=x_{k}^{f}+K_{k}(y_{k}-H_{k}x_{k}^{f})$ (49)

where $K_{k}$ is the Kalman gain and $(y_{k}-H_{k}x_{k}^{f})$ is called the
innovation. After some manipulations, the Kalman gain can be obtained as:

$K_{k}=P_{k}^{f}H^{T}(HP_{k}^{f}H^{T}+R)^{-1}$ (50)

where $P_{k}^{f}$ is the forecast error covariance matrix.

The idea behind EnKF Evensen2003 (Evensen, 2003) is the propagation and
analysis of ensemble members instead of full-rank covariance matrices. Unlike
Kalman filter, the EnKF is applicable for nonlinear dynamical models and
computationally efficient for high dimensional systems. The forecast step of
EnKF can be formulated as follows:

$\displaystyle\hat{v}_{j+1}^{n}$
$\displaystyle=\Psi(v_{j}^{n})+\varepsilon_{j}^{n}\hskip 28.45274ptn=1,...,N$
$\displaystyle\hat{m}_{j+1}$
$\displaystyle=\dfrac{1}{N}\sum_{n=1}^{N}\hat{v}_{j+1}^{n}$ (51)
$\displaystyle\hat{C}_{j+1}$
$\displaystyle=\dfrac{1}{N-1}\sum_{n=1}^{N}(\hat{v}_{j+1}^{(n)}-\hat{m}_{j+1})(\hat{v}_{j+1}^{(n)}-\hat{m}_{j+1})^{T}$

where $v$ is the model state that is propagated by $\Psi$ and
$\varepsilon_{j}^{n}$ are i.i.d’s draws from $N(0,\Sigma)$ or auto-regressive
modeling of the noise. $m$ is the model state mean and $C$ is the model
covariance matrix. Then, the analysis step can be performed as follows:

$\displaystyle S_{j+1}$ $\displaystyle=H\hat{C}_{j+1}H^{T}+\Gamma$
$\displaystyle K_{j+1}$ $\displaystyle=\hat{C}_{j+1}H^{T}S^{-1}_{j+1}$ (52)
$\displaystyle v_{j+1}^{n}$
$\displaystyle=(I-K_{j+1}H)\hat{v}_{j+1}^{(n)}+K_{j+1}y_{j+1}^{(n)}\hskip
14.22636ptn=1,...,N$ $\displaystyle y_{j+1}^{(n)}$
$\displaystyle=y_{j+1}+\eta_{j+1}^{n},\hskip 28.45274ptn=1,...,N$

where $\Gamma$ are i.i.d’s draws from $N(0,\Gamma)$ and $y$ are the
observation states.The ensemble is propagated under the fully nonlinear
dynamics in the forecast step Bocquet2016 (Asch et al., 2016, Stuart et al.,
2015).

###  5.2 Application to Lorenz-96

In this section, we compare the filtering results of Wilk’s (polynomial)
parametrization and compressed sensed parametrization. For each filtering
implementation, 40 ensembles are used and observations are included in every
20 time steps. We observe that compressed sensed dynamics provides better
prediction results compared to other parametrization methods. There are
several softwares that can be used for the filtering problems. Here, we
benefit from the open source provided in PyDA2020 (Pawar et al., 2020)
regarding our parameterized Lorenz-96 models. Observations are generated by
adding normal noise non-parameterized Lorenz-96 model (truth) and filtering
trajectories are obtained by using the set of equations given in (5.1) and
(5.1) for all cases.

####  5.2.1 Autoregressive Wilk’s parametrization

Here, we analyze filtering results of Wilk’s parametrization with auto-
regressive models. Instead of using random noise for the prediction, we use
the auto-regressive modelling of noise obtained in Section 4. While keeping
the parametrized model same, we include modelled noise to the ensemble
dynamics as the following:

$\displaystyle\dot{X}_{k}$
$\displaystyle=\mathcal{G}_{k}(X)-0.03779-0.47362X_{k}+0.005197X_{k}^{2}$
$\displaystyle+0.004882X_{k}^{3}-0.0003142X_{k}^{4}+e_{k}(t)$ (53)

where $\mathcal{G}_{k}(X)$ is defined in (3) and

$e_{k}(t_{i})=0.9453e_{k}(t_{i-1})+0.2265z_{k}(t_{i})$ (54)

$z_{k}(t_{i})$ is sampled randomly from $N(0,1)$.

Figure 5: True trajectories, observations and EnKF with Wilk’s parametrized
model.

Figure 5 shows the trajectories of truth and observations generated from the
truth for the first 3 components of Lorenz-96 model. Figure 7 depicts the
truth and prediction (filtering) results of 40 components and MSPE is 0.00940
which is slightly better than deterministic parametrization.

####  5.2.2 Autoregressive Compressed Sensing

In addition to compressed sensing method in Section 4, we consider auto-
regressive modelling of noise. While keeping the parametrized model same, we
include modelled noise to the ensemble dynamics as the following :

$\displaystyle\dot{X}_{k}$
$\displaystyle=\mathcal{G}_{k}(X)-0.0241-0.474X_{k}+0.00277X_{k}^{2}+0.00509X_{k}^{3}$
$\displaystyle-0.000306X_{k}^{4}+e_{k}(t)$ (55)

where $\mathcal{G}_{k}(X)$ is defined in (3) and

$e_{k}(t_{i})=0.9981e_{k}(t_{i-1})+0.1710z_{k}(t_{i})$ (56)

$z_{k}(t_{i})$ is sampled randomly from $N(0,1)$.

Figure 6: True trajectories, observations and EnKF with compressed sensed
model.

Figure 6 shows the trajectories of truth and observations generated from the
truth for the first 3 components of Lorenz-96 model. Figure 8 depicts the
truth and prediction (filtering) results of 40 components and MSPE is 0.00940
which is slightly better than deterministic parametrization.

In Figures 7 to 8, we see the depictions of filtering results for 40 slow
components of Lorenz96 model. The first graphs of each figure represent the
true trajectories of components where the observations are generated. The
second graphs of each figure show the filtering (prediction) results for the
system. Finally, the third graphs of each figure depict the L1 errors between
the truth and prediction results. All the results about graphs which are shown
in this section are discussed in the previous sections related to the
individual parameterization methods.

Figure 7: Depiction of Lorenz-96 true and EnKF prediction trajectories for 40
components (Auto-regressive Wilk’s parametrization) .

Figure 8: Depiction of Lorenz-96 true and EnKF prediction trajectories for 40
components (Auto-regressive compressed sensing).

Method | MSPE
---|---
Wilk’s parametrization | 0.01006
Compressed sensing | 0.01005
Autoregressive Wilk’s parametrization | 0.00940
Autoregressive compressed sensing | 0.00940
Table 4: Parametrization methods and MSPE results including random noise and
modelled noise with autoregression

Table 4 shows the parameterization methods with the MSPE errors. As a
conclusion, stochastic parameterization methods provide slightly better
prediction results comparing to deterministic parametrizations.

##  6 Conclusion

The sheer number of calculations required in directly solving large-scale
global weather and climate models becomes computationally overwhelming. Hence,
in this paper, using the multi-scale Lorenz-96 system as a testbed, we
introduced a data-driven stochastic-parameterization framework in which the
equations of small-scales are integrated and incorporated using compressed
sensing to reduce the cost of data assimilation. We presented a data-driven,
SPR technique, where we learnt to describe the corresponding unknown or
unresolved physics from the data due to the fact the governing equations of
the unresolved physics are too computationally expensive to solve. Since the
atmospheric model considered in this paper is inherently chaotic, it has
several positive Lyapunov exponents and any small error in the initial
conditions grow in finite time.

Random matrices with independent and identically distributed (i.i.d) entries,
satisfy RIP with high probability, and the RIP plays an important role in our
calculations. We provided a characterization of the dynamical system based on
the Lyapunov exponents for the $F=10$ case. In this case, the model is
strongly chaotic and decorrelation times also suggested the model to be
ergodic. Since the dictionary $\Theta$ was formed as the collection of matrix
of all monomials (stored column-wise) of resolved modes and due to the ergodic
nature of the simulated time series, the entries satisfied the RIP with high
probability. The construction of the optimization problem and computational
method were detailed in section 4.

For parameterization of $U_{k}$, while previous works mostly focused on
suggesting pre-defined polynomial or other model form (or structure) and then
estimating parameters, this work is the first attempt to use CS for
automatically and simultaneously discovering polynomial or model form (or
structure) and estimating parameters for $U_{k}$ as illustrated by the results
in section 4 and Appendix B.

Specifically, a reduced-order ensemble Kalman filter (EnKF) was applied to the
stochastically parametrized model. For the EnKF, particle locations were
corrected based on their respective distances from the observation
(innovation). Magnitude of the correction was proportional to the error
covariance and inverse to the observation noise covariance (Kalman gain). In
between observations, particles in the EnKF sample evolve according to the
original signal dynamics, so particle trajectories are expected to deviate
from the truth as time moves further away from the last observation
correction.

Our results show that the compressed sensing model produces the best
deterministic parametrization results. After incorporating auto-regression to
model residuals, the stochastic parametrizations outperform deterministic
parametrizations in terms of prediction accuracy. Compressed sensing with
auto-regression produces the best stochastic parametrization results.

Solving for small-scale processes with help of Large eddy simulation (LES)
models, on a high-resolution grid embedded (with grid spacing of a few tens of
meters and can explicitly simulate the clouds and turbulence represented by
parameterization schemes) within the low-resolution grid is very useful, but
it’s applicability remains limited to its high computational cost. A long term
goal is to create a data-driven scheme, in which LES models embedded in
selected grid columns of a global model explicitly simulate subgrid-scale
processes which are represented by parameterization schemes in the other
columns. The result is a seamless algorithm which only uses a small scale grid
on a few selected columns that simulates interactively within a running global
simulation, and is far cheaper than eddy-resolving simulation throughout the
global simulation.

## Acknowledgments

The authors acknowledge partial support for this work from Natural Sciences
and Engineering Research Council (NSERC) Discovery grant 50503-10802, TECSIS
/Fields-CQAM Laboratory for Inference and Prediction and NSERC-CRD grant
543433-19. The authors would like to thank Professor Peter Baxendale,
Department of Mathematics, University of Southern California, for stimulating
discussions, contributions and suggestions on many topics of this paper.

## References

  * (1) Anosov, D. V. (1969) Geodesic flows on closed Riemann manifold with negative curvature. Volume 90 of Trundy ordena Lenina Mathematicheskogo Instituta imeni V.A. Steklova. American Mathematical Society.
  * (2) Aydogdu, Y., Mukherjee, A. and Ravichandran, T. (2021). Stochastic Parametrization using Compressed Sensing: Application to the Lorenz-96 Atmospheric Model. GitHub Repository https://github.com/amartyamukherjee/
StochasticParametrization-ApplicationToLorenz96.

  * (3) Asch, M., Bocquet, M., and Nodet, M. (2016). Data assimilation: methods, algorithms, and applications. Society for Industrial and Applied Mathematics.
  * (4) Bengtsson, T., Bickle, P., and Li, B. (2008). Curse-of-dimensionality revisited: Collapse of the particle filter in very large scale systems. IMS Collections: Probability and Statistics Essays in Honor of David A. Freedman, 2, 316–334.
  * (5) Bishop, C. M. (2006). Pattern Recognition and Machine Learning, Springer.
  * (6) Brunton, S. L., Proctor, J. L., and Kutz, J. N. (2016). Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113(15), 3932–3937.
  * (7) Candès, E. J., Romberg, J. K. and Tao, T. (2005). Decoding by linear programming. IEEE Transactions on information theory 51(12), 4203–4215.
  * (8) Candès, E. J. and Tao, T. (2006). Near-optimal signal recovery from random projections: Universal encoding strategies. IEEE Transactions on Information Theory 52(12), 5406–5425.
  * (9) Candès, E. J. and Wakin, M. B. (2006). An introduction to compressed sampling. IEEE signal processing magazine 25.2, 21–30.
  * (10) Chen, S. S., Donoho, D. L. and Saunders, M. A. (2003). Atomic decomposition by basis pursuit. SIAM Journal on Scientific Computing 20(1), 33–61.
  * (11) Devaney, R. L. (1989). An Introduction to Chaotic Dynamical Systems, Addison-Wesley.
  * (12) Donoho, D. L., (2006). Compressed sensing. IEEE Transactions on information theory 52.4, 1289–1306.
  * (13) Eckmann, P., and Ruelle, D. (1985). Ergodic theory of chaos and strange attractors. Rev. Mod. Phys., 57(3), 617-656.
  * (14) Evensen, G.,(1994) Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research, 99,10143–10162.
  * (15) Evensen, G.,(2003) The ensemble Kalman filter: Theoretical formulation and practical implementation. Ocean dynamics, 53(4):343–367, .
  * (16) Gagne, D. J., Christensen, H., Subramanian, A., and Monahan, A. H. (2020). Machine learning for stochastic parameterization: Generative adversarial networks in the lorenz ’96 model. Journal of Advances in Modeling Earth Systems, 12(3).
  * (17) Gordon, N.J., Salmond, D.J. and Smith, A. F. M. (1993). Novel approach to nonlinear/non-gaussian Bayesian state estimation. IEEE Proceedings, F(2), 107-113.
  * (18) Kalman, R. E. (1960). A new approach to linear filtering and prediction problems.
  * (19) Kullback, S. and Leibler, R. A. (1951). On Information and Sufficiency. The Annals of Mathematical Statistics 22, 79–86.
  * (20) Law, K., Stuart, A., and Zygalakis, K. (2015). Data assimilation. Cham, Switzerland: Springer, 214.
  * (21) Lorenz, E. N. (1995). Predictability: A problem partly solved. ECMWF Seminar Proceedings on Predictability, ECMWF.
  * (22) Majda, A., Abramov, R. and Grote, M. (2005). Information Theory and Stochastic for Multiscale Nonlinear Systems (CRM Monograph Series), 25, American Mathematical Society, New York.
  * (23) Pawar, A. S. E., and San, 0. (2020). A Hands-On Introduction to Dynamical Data Assimilation with Python. Fluids 5(4), 225\.
  * (24) Pesin Y B, (2004), Lectures on partial hyperbolicity and stable ergodicity. European Mathematical Society
  * (25) Pindyck, R. S. and Rubinfeld, D. L. (1991). Forecasting with Time-Series Models. Econometric Models & Economic Forecasts 3rd ed., 516–535.
  * (26) Sinai, Y. G. (1977). Introduction to ergodic theory. volume 18 of Mathematical Notes, Princeton University Press.
  * (27) Skokos, Ch. (2010). The Lyapunov Characteristic Exponents and Their Computation. Lect. Notes Phys., 790, 63–135.
  * (28) Snyder, C., Bengtsson, T., Bickel, P., and Anderson, J. (2008). Obstacles to high-dimensional particle filtering. Monthly Weather Review, 136(12), 4629–4640.
  * (29) Tibshirani, R., (1996). Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical Society. Series B (Methodological) 58.1, 267–288.
  * (30) Wilks, D. S., (2005). Effects of stochastic parametrizations in the Lorenz ’96 system. Royal Meteorological Society 131, 389–407.
  * (31) Wolf, A. S., Swinney, J. B. and Vastano, H. L. (1985). Determining Lyapunov exponents from a time series. Physica 16D, 285–317.
  * (32) Young, L. S., (1998). Statistical properties of dynamical systems with some hyperbolicity. The Annals of Mathematics 147(3), 585–650.

## Appendix A: Auto-Regression

Autoregressive (AR) processes are used to predict residuals given a time
series of previous residual values.

The first-order AR model is given by (13). This can also be written as:

$e_{k}(t_{i})=\sigma z_{k,i}+\varphi\sigma z_{k,i-1}+\varphi^{2}\sigma
z_{k,i-2}+...$ (57)

Given the time series $e_{k}(t)$, $\varphi$ is used to approximate
$e_{k}(t_{i})$ given $e_{k}(t_{i-1})$. $\varphi$ is estimated using a least
squares minimization process shown below:

$C(\varphi)=\sum_{k=1}^{K}\sum_{i=2}^{I}(e_{k}(t_{i})-\varphi
e_{k}(t_{i-1}))^{2}$ (58)

In order for the process to be stable, $|\varphi|<1$ must be satisfied so that
the time series $e_{k}(t_{i})$ shown in (57) does not reach large values. In
(57), $\sigma$ is used to measure the variability in $e_{k}(t_{i})$ given
$\varphi e_{k}(t_{i-1})$. From calculating $\varphi$, we can estimate $\sigma$
by dividing the cost function by the degrees of freedom. This is shown by the
following equation:

$\sigma^{2}=\frac{C(\varphi)}{K(I-1)-1}$ (59)

We then calculate the variance of the process
$\sigma_{e}^{2}=Var(e_{k}(t_{i}))$ using (57).

$\displaystyle\sigma_{e}$
$\displaystyle=\sigma^{2}Var(z_{k,i})+\varphi^{2}\sigma^{2}Var(z_{k,i-1})$
$\displaystyle+\varphi^{4}\sigma^{2}Var(z_{k,i-2})+...$
$\displaystyle=\sigma^{2}+\varphi^{2}\sigma^{2}+\varphi^{4}\sigma^{2}+...$
$\displaystyle\approx\frac{\sigma^{2}}{1-\varphi^{2}}$ (60)

The AR model is used in addition to regression and compressed sensing in order
to further improve the prediction of resolved variables. The estimates for
$\varphi,\sigma,\sigma_{e}$ are all important to forming and evaluating the AR
model in regression and compressed sensing presented in section 5.2.

Appendix B: Degree-4 Compressed Sensing Model
---
$U_{1}=-0.467216X_{1}+0.003429X_{1}^{2}+0.004481X_{1}^{3}-0.000262X_{1}^{4}-0.024885$
$U_{2}=-0.47689X_{2}+0.0002X_{2}^{2}+0.005351X_{2}^{3}-0.000299X_{2}^{4}-0.021728$
$U_{3}=-0.478369X_{3}+0.002902X_{3}^{2}+0.005584X_{3}^{3}-0.000358X_{3}^{4}-0.021464$
$U_{4}=-0.472069X_{4}+0.002041X_{4}^{2}+0.005191X_{4}^{3}-0.000307X_{4}^{4}-0.021703$
$U_{5}=-0.481931X_{5}+0.001638X_{5}^{2}+0.005625X_{5}^{3}-0.000345X_{5}^{4}-0.03192$
$U_{6}=-0.486479X_{6}+0.001797X_{6}^{2}+0.006035X_{6}^{3}-0.000385X_{6}^{4}-0.029337$
$U_{7}=-0.46724X_{7}+0.001871X_{7}^{2}+0.004872X_{7}^{3}-0.000281X_{7}^{4}-0.025702$
$U_{8}=-0.487255X_{8}+0.002038X_{8}^{2}+0.00592X_{8}^{3}-0.000368X_{8}^{4}-0.030243$
$U_{9}=-0.466759X_{9}+0.003636X_{9}^{2}+0.004494X_{9}^{3}-0.000258X_{9}^{4}-0.027036$
$U_{10}=-0.463411X_{10}+0.002444X_{10}^{2}+0.004413X_{10}^{3}-0.00024X_{10}^{4}-0.013236$
$U_{11}=-0.478315X_{11}+0.004852X_{11}^{2}+0.004976X_{11}^{3}-0.000318X_{11}^{4}-0.016471$
$U_{12}=-0.467606X_{12}+0.004247X_{12}^{2}+0.004642X_{12}^{3}-0.00028X_{12}^{4}-0.024363$
$U_{13}=-0.46995X_{13}+0.002823X_{13}^{2}+0.004957X_{13}^{3}-0.000297X_{13}^{4}-0.018118$
$U_{14}=-0.470231X_{14}+0.001179X_{14}^{2}+0.004991X_{14}^{3}-0.00028X_{14}^{4}-0.019603$
$U_{15}=-0.48344X_{15}+0.003052X_{15}^{2}+0.005955X_{15}^{3}-0.00039X_{15}^{4}-0.023006$
$U_{16}=-0.470097X_{16}+0.003608X_{16}^{2}+0.004726X_{16}^{3}-0.000275X_{16}^{4}-0.030639$
$U_{17}=-0.473192X_{17}+0.001691X_{17}^{2}+0.005294X_{17}^{3}-0.000313X_{17}^{4}-0.029741$
$U_{18}=-0.47504X_{18}+0.004303X_{18}^{2}+0.005322X_{18}^{3}-0.000345X_{18}^{4}-0.022838$
$U_{19}=-0.473415X_{19}+0.002968X_{19}^{2}+0.004922X_{19}^{3}-0.000289X_{19}^{4}-0.011408$
$U_{20}=-0.474502X_{20}+0.003374X_{20}^{2}+0.005137X_{20}^{3}-0.000306X_{20}^{4}-0.032649$
$U_{21}=-0.472923X_{21}+0.003193X_{21}^{2}+0.005178X_{21}^{3}-0.000327X_{21}^{4}-0.018312$
$U_{22}=-0.466437X_{22}+0.001616X_{22}^{2}+0.005185X_{22}^{3}-0.000318X_{22}^{4}-0.017406$
$U_{23}=-0.474267X_{23}+0.003114X_{23}^{2}+0.005215X_{23}^{3}-0.000323X_{23}^{4}-0.020402$
$U_{24}=-0.465558X_{24}+0.002984X_{24}^{2}+0.004562X_{24}^{3}-0.000262X_{24}^{4}-0.03494$
$U_{25}=-0.474495X_{25}+0.002997X_{25}^{2}+0.005042X_{25}^{3}-0.000302X_{25}^{4}-0.018418$
$U_{26}=-0.468641X_{26}+0.003332X_{26}^{2}+0.004428X_{26}^{3}-0.000249X_{26}^{4}-0.025618$
$U_{27}=-0.479017X_{27}+0.002635X_{27}^{2}+0.005268X_{27}^{3}-0.000322X_{27}^{4}-0.028429$
$U_{28}=-0.466174X_{28}+0.002651X_{28}^{2}+0.004866X_{28}^{3}-0.000289X_{28}^{4}-0.034918$
$U_{29}=-0.476855X_{29}+0.002663X_{29}^{2}+0.005199X_{29}^{3}-0.000311X_{29}^{4}-0.022651$
$U_{30}=-0.478937X_{30}+0.001475X_{30}^{2}+0.00543X_{30}^{3}-0.00032X_{30}^{4}-0.017885$
$U_{31}=-0.45929X_{31}+0.004295X_{31}^{2}+0.00434X_{31}^{3}-0.00026X_{31}^{4}-0.024998$
$U_{32}=-0.477136X_{32}+0.003377X_{32}^{2}+0.005254X_{32}^{3}-0.000323X_{32}^{4}-0.029709$
$U_{33}=-0.472517X_{33}+0.001649X_{33}^{2}+0.004832X_{33}^{3}-0.000264X_{33}^{4}-0.027171$
$U_{34}=-0.470095X_{34}+0.004248X_{34}^{2}+0.004854X_{34}^{3}-0.000301X_{34}^{4}-0.018149$
$U_{35}=-0.476457X_{35}+0.005394X_{35}^{2}+0.005025X_{35}^{3}-0.000326X_{35}^{4}-0.022066$
$U_{36}=-0.4798X_{36}+0.000396X_{36}^{2}+0.005354X_{36}^{3}-0.000286X_{36}^{4}-0.02986$
$U_{37}=-0.468947X_{37}+0.00297X_{37}^{2}+0.004916X_{37}^{3}-0.000289X_{37}^{4}-0.031595$
$U_{38}=-0.476784X_{38}+0.002019X_{38}^{2}+0.005193X_{38}^{3}-0.000302X_{38}^{4}-0.032644$
$U_{39}=-0.480576X_{39}+0.004569X_{39}^{2}+0.005093X_{39}^{3}-0.000322X_{39}^{4}-0.017789$
$U_{40}=-0.482054X_{40}+0.001322X_{40}^{2}+0.005556X_{40}^{3}-0.000332X_{40}^{4}-0.020528$