Original Article Daniel Ayers<EMAIL_ADDRESS>UK Engineering and
Physical Sciences Research Council, Grant Number: EP/N509723/1; UK National
Centre for Earth Observation Grant Number: NCEO02004; Schmidt Futures.

# Supervised machine learning to estimate instabilities in chaotic systems:
estimation of local Lyapunov exponents

Daniel Ayers  Department of Meteorology, University of Reading, Reading, UK
UK National Centre for Earth Observation, Reading, UK  Jack Lau  Department of
Computer Science, University of Reading, Reading, UK  Javier Amezcua
Department of Meteorology, University of Reading, Reading, UK  Tecnologico de
Monterrey, Campus Ciudad de Mexico, Mexico City, Mexico  Alberto Carrassi
Department of Meteorology, University of Reading, Reading, UK  Department of
Physics and Astronomy “Augusto Righi”, University of Bologna, Italy  Varun
Ojha  Department of Computer Science, University of Reading, Reading, UK
School of Computing, Newcastle University, UK

###### Abstract

In chaotic dynamical systems such as the weather, prediction errors grow
faster in some situations than in others. Real-time knowledge about the error
growth could enable strategies to adjust the modelling and forecasting
infrastructure on-the-fly to increase accuracy and/or reduce computation time.
For example one could change the ensemble size, or the distribution and type
of target observations, etc. Local Lyapunov exponents are known indicators of
the rate at which very small prediction errors grow over a finite time
interval. However, their computation is very expensive: it requires
maintaining and evolving a tangent linear model, orthogonalisation algorithms
and storing large matrices.

In this feasibility study, we investigate the accuracy of supervised machine
learning in estimating the current local Lyapunov exponents, from input of
current and recent time steps of the system trajectory, as an alternative to
the classical method. Thus machine learning is not used here to emulate a
physical model or some of its components, but “non intrusively” as a
complementary tool. We test four popular supervised learning algorithms:
regression trees, multilayer perceptrons, convolutional neural networks and
long short-term memory networks. Experiments are conducted on two low-
dimensional chaotic systems of ordinary differential equations, the Rössler
and the Lorenz 63 models. We find that on average the machine learning
algorithms predict the stable local Lyapunov exponent accurately, the unstable
exponent reasonably accurately, and the neutral exponent only somewhat
accurately. We show that greater prediction accuracy is associated with local
homogeneity of the local Lyapunov exponents on the system attractor.
Importantly, the situations in which (forecast) errors grow fastest are not
necessarily the same as those where it is more difficult to predict local
Lyapunov exponents with machine learning.

###### keywords:

supervised machine learning, local Lyapunov exponents, numerical modelling,
chaos

##  1 Introduction

Weather and climate are well known exemplars of chaotic dynamical systems.
These systems exhibit extreme sensitivity to initial conditions, meaning that
initial condition errors are subject to (on average) exponential growth until
they reach saturation Lorenz1963 (47, 39). The rate and the characteristic of
such growth, however, is highly state dependent Lighthill1986 (46, 85). As a
consequence, although chaotic systems have a finite predictability horizon
(about $2$ weeks for the atmosphere, see e.g. Holton2013 (36)), the best
estimate of prediction error growth fluctuates in size along with the system’s
evolution, as the system goes through periods of lower or higher
predictability. For example, the short-term predictability of the atmosphere
depends on the weather regime present at a given time (see, e.g. Palmer1996,
59). Understanding the nature of error growth is essential to characterising a
system, and to enable better prediction. The present work is motivated by the
idea that if the degree of predictability of the system is known in real-time,
it may be possible and beneficial to take adaptive measures. For instance, we
speculate that a local decrease of predictability might be counteracted by
increasing the ensemble size in the context of ensemble-based data
assimilation or probabilistic forecasting, or the distribution and type of
target observations. Conversely, in areas of high predictability, one might
save computational resources (and thus energy consumption) via the opposite
actions. Understanding the impact of such actions would require
experimentation. In this study we investigate the potential of machine
learning (ML) methods Bishop1995 (7, 32) to provide a real-time estimation of
the system’s local predictability.

The mathematical theory of dynamical systems has long been the backbone to
understand and quantify predictability in deterministic chaotic systems. This
is commonly done by studying the instability properties of the solution, i.e.
by analysing the linearised dynamics of small perturbations: the tangent space
evolution of these “small” perturbations are taken as proxies of the dynamics
of unknown initial condition errors ott_2002 (57). In this context, Lyapunov
exponents (LEs) are well-established quantities that measure the asymptotic
rates of error growth for a set of infinitesimally small errors that capture
all directions of phase space Pikovsky2016 (64). In practice, they measure the
average growth of small finite errors over long periods of time. The spectrum
of LEs is characteristic of each given dynamical system. Lyapunov exponents,
and their corresponding Lyapunov vectors (LVs), have been exploited in
geosciences for more efficient uncertainty quantification in data assimilation
(e.g., Palatella2013lyapunov, 58, 65, 1, 17), or for initialising
probabilistic predictions (e.g., Toth1997, 83, 15, 86). The LE spectrum can
also be used to calculate other characteristic properties of a system, such as
the Kolmogorov-Sinai entropy, measuring the rate of information loss Sinai2009
(76), or the Kaplan-Yorke attractor dimension Kaplan1979 (40).

We note that LEs are associated with directions known as covariant Lyapunov
vectors (CLVs). CLVs also provide useful information and can be calculated
numerically (see Ginelli2007 (25, 92, 24)). However, in this work, we focus on
the exponents only.

The LEs are calculated as an average of finite-time Lyapunov exponents which
are here referred to as local Lyapunov exponents, LLEs Benettin1980theory (5,
6, 41). Whereas LEs provide “global” information about the average growth of
small perturbations in the system, the LLEs describe “local” growth rates
along a finite-time section of the trajectory. Notably, the LLEs show the
heterogeneity of the instabilities in phase space: the fluctuation of the
local dynamical stability around the asymptotic value as the system state
varies Sandri1996 (74, 64). This makes the LLEs ideal quantities to measure
the local degree of predictability, yet a bottleneck for their real-time use
in operational scenarios is the huge computational cost. Computational cost
grows quickly with the system’s dimension making it prohibitive even for
moderate size models, let alone for models as large as those currently used in
numerical weather predictions ($\mathcal{O}(10^{9})$ dimensions). Using the
standard method Benettin1980theory (5, 6, 41, 64), calculating the LLEs and
LEs involves computing a long trajectory of the system (including a spin-up
needed to ensure the solution has reached the model attractor), propagating
perturbations (as many as the number of desired LLEs) with the tangent linear
model (i.e., the resolvent of the model Jacobian), and then repeatedly
performing a process of orthogonalisation (e.g. using a QR decomposition
algorithm).

Despite the computational bottleneck, Lyapunov methods (i.e. computing the
local and global LE spectrum and aforementioned associated properties, or the
Lyapunov vectors) have been used for dynamical analysis of geophysical models
of intermediate order ($\mathcal{O}(10^{3})$ to $\mathcal{O}(10^{5})$
variables), for example see Vannitsem2016 (87, 85, 21). Additionally, Lyapunov
methods have been applied to weather reanalysis data to analyse the dynamics
of the North Atlantic Oscillation Quinn2021 (66) and of persistent states of
atmospheric pressure over the European and western Asian continents Quinn2022
(67). In these works, the bottleneck was overcome by reducing the data
dimension (using empirical orthogonal functions) and constructing a reduced
model. Whilst these works demonstrate the utility of Lyapunov methods, they do
not provide a means of calculating LLEs that is cheap enough to be carried out
regularly during a forecasting cycle.

Avoiding the need for such a costly computation whilst attaining an estimate
of the LEs or LLEs has thus great relevance. In their recent work,
Chen_et_al21 (19) show how the outcomes of properly tuned data assimilation
experiments can reveal the first LE as well as the Kolmogorov-Sinai entropy of
the underlying dynamical model. The present work also seeks to avoid the cost
of classical calculation methods, albeit only when time is critical. We
investigate the feasibility of using certain ML methods Bishop1995 (7, 32) to
estimate the LLEs based only on information from the system’s solution. Our
focus is on supervised learning, which uses a data set of input-output pairs.
The targets, i.e. the desired outputs, are LLEs calculated using the classical
method of evolving perturbations via the tangent linear model and
orthogonalising. In this way, the cost of such methods is paid during the
training phase of the ML method, and is avoided when making predictions.

In the area of weather and climate forecasting, supervised learning has been
used for various purposes (see e.g., Reichstein2019, 71, 69, 18, 22, and
references therein). These include: (i) to emulate the full dynamics of a
system Pathak2017 (62, 61, 23, 53, 12, 60, 75, 78), (ii) to improve a physics-
based model with data-driven correction or parameterisation OGorman2018 (55,
70, 9, 68, 10, 54, 28, 13). Both approaches imply an intervention on the
original model: the first approach yields surrogate data-driven models of the
full original system, while the second approach builds hybrid models with
data-driven and physics-based components. In either case the spectrum of the
LEs of these new models can be computed using the standard approach
Benettin1980theory (5, 6, 41), and can be compared to that of the original
model as a way to quantify the goodness of the ML-reconstructed dynamics
Pathak2017 (62, 12).

In contrast to these two families of methods, this study aims towards
improving prediction skills by equipping the model with an external tool to
quantify in real-time the local degree of predictability, and thus guide “non
intrusive” adaptations whereby the model equations are left unaltered. More
specifically, the goal is to use ML to predict the current LLE spectrum given
input of the system state at the current and (possibly) most recent time
steps. We envisage that the trained ML algorithm could then be interrogated
for information about the local dynamical instability whilst performing the
numerical model forward integration. We speculate that such information could
drive a decision process for adaptive modelling, for example adjusting the
ensemble size when doing ensemble-based data assimilation or probabilistic
predictions, changing the distribution and type of target observations, or
adapting the numerical integration scheme. Such adaptations could mitigate
error, improve uncertainty quantification, or reduction computational cost.

In this feasibility study, we test the accuracy of some popular supervised ML
algorithms in this task in two prototypical low-dimensional chaotic dynamical
systems. This study is concerned solely with the predictive capability of ML
methods: the task of optimising the computational cost of making predictions
is left for future work. We anticipate that the latter task will be largely
dependent on the specific use-case and computing hardware. The ML algorithms
we test are regression trees (RTs) breiman1984 (14), multilayer perceptrons
(MLPs) (e.g. see Goodfellow-et-al-2016, 27, Chapter 6), convolutional neural
networks (CNNs) LeCun1990 (43), and long short-term memory networks (LSTMs)
Hochreiter1997 (35, 29). These algorithms encompass three approaches to
exploiting the temporal structure of the input. We evaluate both their point-
wise accuracy and their statistical performance, measured in this case by the
closeness of the distribution of predictions to the distribution of the target
values. We find that on average the machine learning algorithms predict the
stable local Lyapunov exponent accurately, the unstable exponent reasonably
accurately, and the neutral exponent only somewhat accurately. Each exponent
is predicted more accurately in the Lorenz 63 system than in the Rössler
system. We show that greater prediction accuracy is associated with local
homogeneity of the local Lyapunov exponents on the system attractor.
Importantly, the situations in which (forecast) errors grow fastest are not
necessarily the same as those where it is more difficult to predict local
Lyapunov exponents with machine learning.

The rest of this paper is organised as follows. In Section 2, we briefly
review the theory of LEs and detail the method used to compute them. In
Section 3 we pose and conceptualise the ML problem we intend to solve,
motivate the choice of the algorithms, and detail the input and target data
and the evaluation metrics. In Section 4 we present the two systems under
consideration: the Rössler and the Lorenz 63 models and discuss the
characteristics of their Lyapunov spectra. In Section 5 we present the results
and Section 6 concludes with a discussion.

##  2 Lyapunov exponents

###  2.1 Overview of the theory

We review briefly the theory of LEs, with the aim of providing an intuitive
explanation of what they are. The section follows Legras1996 (44, 5, 6, 41,
64) and (Strogatz2018, 81, Chapter 9, Section 3), to which we refer the reader
for a more rigorous and comprehensive treatment.

Consider a deterministic autonomous dynamical system

$\displaystyle\dot{\mathbf{x}}=\mathbf{g}(\mathbf{x}),$ (1)

where $\mathbf{x}\in\mathbb{R}^{n}$ is the state of the system,
$\mathbf{g}\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ is the evolution function,
and $\dot{\mathbf{x}}$ denotes the derivative of $\mathbf{x}$ with respect to
time. A trajectory of the dynamical system starting at initial condition
$\mathbf{x}(0)$ is a set $\\{\mathbf{x}(t)\colon t\in A\\}$, where $A$ is a
connected subset of $\mathbb{R}_{\geq 0}$ containing $0$. Consider the
difference $\mathbf{v}(t)$ between a trajectory started from the “true”
initial condition $\mathbf{x}(0)$, and a trajectory started from the perturbed
initial condition $\mathbf{x}(0)+\mathbf{v}(0)$, where $\mathbf{v}(0)$ is
infinitesimally small. The idea behind LEs is to find $\lambda(t)$, where

$\displaystyle e^{t\lambda(t)}$
$\displaystyle=\frac{\|\mathbf{v}(t)\|}{\|\mathbf{v}(0)\|}$ (2a)
$\displaystyle\Leftrightarrow\lambda(t)$
$\displaystyle=t^{-1}\ln\left(\frac{\|\mathbf{v}(t)\|}{\|\mathbf{v}(0)\|}\right).$
(2b)

In other words, $\lambda(t)$ is the exponential growth rate of the initial
error. In this setting, $\lambda(t)$ is specific to the initial condition
$\mathbf{x}(0)$ and perturbation $\mathbf{v}(0)$, and $\lambda(t)$ varies with
time.

Lyapunov exponents generalise this notion to describe a) average exponential
growth rates of the system, regardless of the initial condition, and b) the
exponential growth rates for infinitessimal perturbations in all directions.
To account for all directions, we consider perturbations contained in an
$n$-sphere of infinitesimal radius. As time progresses, perturbations within
the sphere are mapped into an ellipsoid. The LEs are the time-average of the
exponential growth rate of the ratios between the axes of the sphere and of
the ellipsoid Legras1996 (44).

Fix an initial condition $\mathbf{x}(0)$ and let $\mathbf{v}(0)$ be an
infinitesimally small perturbation, as above. Then the dynamics of the
perturbation are given by

$\displaystyle\dot{\mathbf{v}}=\mathbf{J_{g}}\mathbf{v},$ (3)

where $\mathbf{J_{g}}$ is the Jacobian of $\mathbf{g}$ evaluated at
$\mathbf{x}(t)$, i.e. the linearisation of the evolution function at
$\mathbf{x}(t)$. The solutions to Equation 3 can be found using a fundamental
matrix, i.e. any matrix-valued function $\mathbf{M}(t)$ satisfying

$\displaystyle\dot{\mathbf{M}}=\mathbf{J_{g}}\mathbf{M},$ (4)

and such that $\mathbf{M}(t)$ is non-singular for all $t$. Focusing on a
single perturbation, it follows from Equation 2a that

$\displaystyle
e^{\lambda(t)}=\left(\frac{\|\mathbf{v}(t)\|}{\|\mathbf{v}(0)\|}\right)^{\frac{1}{t}}=\left(\frac{\sqrt{\mathbf{v}(0)^{\rm{T}}\mathbf{M}(t)^{\rm{T}}\mathbf{M}(t)\mathbf{v}(0)}}{\sqrt{\mathbf{v}(0)^{\rm{T}}\mathbf{v}(0)}}\right)^{\frac{1}{t}}.$
(5)

We are interested in the value of $\lambda(t)$ as $t\to\infty$. Rearranging
Equation 5 and taking the limit, we have:

$\displaystyle\lambda=\lim_{t\to\infty}\ln\left[\left(\mathbf{v}(0)^{\rm{T}}\mathbf{M}(t)^{\rm{T}}\mathbf{M}(t)\mathbf{v}(0)\right)^{\frac{1}{2t}}\right].$
(6)

For almost all choices of $v(0)$, the $\lambda$ given by Equation 6 is the
largest LE.

We consider now the full spectra of LEs $\lambda_{i}(t),\ i=1,\dots,n$ that
arise when one considers a sphere of perturbations. The growth of a sphere of
perturbations depends only on $\mathbf{M}(t)^{\rm{T}}\mathbf{M}(t)$. Thus we
consider the limit $\mathbf{W}(\mathbf{x}(0))$ defined by

$\displaystyle\mathbf{W}(\mathbf{x}(0))=\lim_{t\to\infty}\left[\mathbf{M}(t)^{\rm{T}}\mathbf{M}(t)\right]^{\frac{1}{2t}}.$
(7)

By the multiplicative ergodic theorem Oseledets1968 (56, 73) the limit exists,
depends on the initial condition $\mathbf{x}(0)$, and importantly, the
eigendecomposition of $\mathbf{M}(t)^{\rm{T}}\mathbf{M}(t)$ in the limit
exists, which gives

$\displaystyle\mathbf{W}(\mathbf{x}(0))=\mathbf{P}(\mathbf{x}(0))\mathbf{D}\mathbf{P}^{\rm{T}}\\!(\mathbf{x}(0)),$
(8)

where the eigenvector matrix $\mathbf{P}(\mathbf{x}(0))$ is orthonormal. The
matrix of eigenvalues $\mathbf{D}$ is unique and depends neither on
$\mathbf{x}(0)$ nor on the norm of the vector space containing the
perturbations Kuptsov2012 (41). The LEs,
$\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{n}$, are the natural
logarithm of the diagonal elements of $\mathbf{D}$.

We finish with some remarks on the significance of the LEs. A chaotic system
is a system with at least one positive LE. The LEs, above defined in terms of
the growth of axes of a sphere of perturbations, are linked to the growth of
the volume of the $n$-parallelepiped defined by the principle axes of the
resulting ellipsoid (see Wolf1985, 91). Also, the sum of the LEs is equal to
the average divergence of the flow (see Pikovsky2016, 64, Section 2.5.4).
Thus, in dissipative systems, the sum of the LEs is negative. Finally,
continuous chaotic systems have at least one LE equal to zero. This is due to
there being zero growth of an infinitesimal perturbation in the direction of
the flow.

###  2.2 Computation of local and global Lyapunov exponents

The theory does not translate directly to a method for calculation of the LEs,
since in order to approach the limit in Equation 7 one must integrate Equation
4 to find $\mathbf{M}(t)$ for very large $t$. This both accumulates numerical
errors and results in a range eigenvalues of
$\mathbf{M}^{\rm{T}}\\!(t)\mathbf{M}(t)$ that is too large for accurate
numerical calculation Pikovsky2016 (64). Instead, we measure the growth of
perturbations over (finitely) many small time intervals, and compute the
average. Specifically, for each time interval we calculate the LLEs: the
natural logarithm of the growth ratios, divided by the length of the time
interval, as shown in Equation 2b. If the system is ergodic, the arithmetic
mean of the LLEs converges to the LEs as the number of time intervals
increases. Crucially, the perturbation vectors are orthogonalised and resized
between each time interval. Orthogonalising the propagated perturbations is
necessary to keep the perturbations distinct, since perturbations will tend to
be attracted towards the direction of largest growth. The resizing is
necessary to prevent perturbations becoming too small or large to be
represented by floating-point numbers. We now present the algorithm used to
calculate LEs and LLEs in this work, which is based on methods presented in
Benettin1980theory (5, 6, 41).

  1. 1.

Calculate and store a long trajectory $\\{\mathbf{x}(t)\colon
t\in[0,T_{\rm{end}}]\\}$. Discard an initial transient period to ensure the
trajectory is in the attractor. One can alternatively calculate the trajectory
at the same time as integrating Equation 4 in Step 33a below, which avoids the
need to store a long trajectory.

  2. 2.

Initialise a matrix of perturbations
$\mathbf{Q}_{0}=[\mathbf{q}_{0}^{1},\mathbf{q}_{0}^{2},\dots,\mathbf{q}_{0}^{n}]$,
such that the $\mathbf{q}_{0}^{i}\in{\mathbb{R}}^{n}$ are orthogonal and of
unit length, that is, orthonormal.

  3. 3.

Repeat the following iteration $m$ times, where $m$ is large enough to achieve
convergence of the LEs. In iteration $j$, starting with $j=1$, perturbations
are propagated along the trajectory $\\{\mathbf{x}(t)\colon
t\in[(j-1)\tau,j\tau]\\}$, where $\tau$ is typically small. Each iteration
results in $n$ LLEs: $\lambda_{j}^{i},\ i=1,\dots,n$. Henceforth we notate LEs
with hatted lambdas to distinguish the asymptotic LE $\hat{\lambda}_{i}$ from
the LLE $\lambda_{j}^{i}$.

    1. (a)

Propagate the perturbations:
$\mathbf{V}_{j}=\mathbf{M}(j\tau)\mathbf{Q}_{j-1}$, where $\mathbf{Q}_{j-1}$
is from iteration $j-1$, and $\mathbf{M}(j\tau)$ is computed by integrating
Equation 4.

    2. (b)

Orthonormalise the propagated perturbations $\mathbf{V}_{j}$ to get
$\mathbf{Q}_{j}$ using QR-decomposition Golub2013 (26, 80):

$\displaystyle\mathbf{Q}_{j}\mathbf{R}_{j}=\mathbf{V}_{j}.$ (9)

    3. (c)

The diagonal elements $r_{j}^{i}$ of $\mathbf{R}_{j}$ are the desired ratios.
The LLEs at time $j\tau$ are calculated as:

$\displaystyle\lambda_{j}^{i}=\tau^{-1}\ln(r_{j}^{\alpha(i)}),\ \
i=1,\dots,n.$ (10)

In Equation 10 the diagonal element $r_{j}^{\alpha(i)}$ is indexed by
labelling function $\alpha(i)$, where $\alpha(i)$ is determined in Step 4.

  4. 4.

Calculate the LEs:

$\displaystyle\hat{\lambda}_{i}=(m-k)^{-1}\sum_{j=k+1}^{m}\lambda_{j}^{i}\quad\quad\quad
i=1,\dots,n$ (11)

where the LLEs from the first $k$ iterations are discarded. The bijective
function $\alpha$ (in Equation 10) takes inputs and values $1,\dots,n$, and is
chosen such that the global LEs are numbered in descending order, i.e. such
that $\hat{\lambda}_{i}\geq\hat{\lambda}_{i+1}$. We say $i^{\rm{th}}$ LLE to
refer to any set of LLEs $\\{\lambda_{j}^{i}\colon j=k,\dots,m\\}$ that are
associated to the $i^{\rm{th}}$ LE.

In Step 4, a transient period of length $k$ iterations is required to let the
initial perturbations $\mathbf{Q}_{0}$ converge to the dynamics of the
trajectory so that the leading perturbation $\mathbf{q}_{j}^{1}$ is oriented
in the direction of largest growth. As discussed above, the LLEs are defined
in terms of ratios of the axes of the $n$-sphere and the ellipsoid. In
practice, it is unlikely that the chosen initial perturbation
$\mathbf{q}_{0}^{1}$ will be mapped by $\mathbf{M}(\tau)$ onto the leading
axis of the ellipsoid, which will lead to a poor estimation of the LLE.
However, with sufficiently many iterations $k$, $\mathbf{q}_{k}^{1}$ will be
attracted to the direction of largest growth, leading to more accurate
estimates.

#### Computational cost

The algorithm for computing the LLEs and LEs does not scale well. The
computational cost of the steps of the algorithm are as follows. The length of
the required transient period, i.e. the number of "spin-up" iterations,
depends on the system dynamics and, in the worst case, grows proportionally to
the system dimension $n$ (Kuptsov2012, 41, p. 754). For computing LEs, the
total number of iterations $m$ depends on the complexity of the attractor and
the precision required. Each iteration (Step 3) requires the calculation of a
trajectory of length $\tau$, which involves at least $\mathcal{O}(n)$ floating
point operations (flops). The cost of integrating the matrix differential
Equation 4 involves at least one evaluation of the Jacobian matrix, and at
least one multiplication of the Jacobian by another matrix, per timestep. For
a dense, non-trivial Jacobian, it is reasonable to assume that Step 3a
involves $\mathcal{O}(n^{2})$ flops in the simplest case, namely where
$\tau=\Delta t$ and a minimal numerical integration scheme is used. In a less
simple case, Step 3a will require at least $\mathcal{O}(n^{2.3})$ flops: the
cost of multiplying two dense $n\times n$ matrices Alman2021 (3). Step 3b is
by far the most expensive step, since computing eigenvectors and eigenvalues
via QR decomposition requires $25n^{3}$ flops Golub2013 (26, 4). The overall
theoretical time complexity of the LE algorithm is thus
$25n^{3}+\mathcal{O}(n^{2})$ flops. In practice, $n\times n$ matrix
multiplication can be much slower (at least $\mathcal{O}(n^{3})$) due to
memory access latency Albrecht2010 (2). Consequently, computing the full LLE
spectrum for a modern weather prediction system (where $n\approx 10^{9}$) is
too expensive to be done during the forecast cycle.

Obviously, generating a subset of the LLE spectrum costs less. When using the
tangent linear model, one must compute consecutive leading LLEs: it is not
possible to calculate LLE $\lambda^{i}$ without also calculating
$\lambda^{1},\dots,\lambda^{i-1}$. The cost of QR decomposition of an $n\times
i$ matrix scales at $\mathcal{O}(ni^{2})$ Boyd2018 (11). The cost of
multiplying an $n\times n$ matrix by an $n\times i$ matrix is
$\mathcal{O}(n^{2})$ when $i$ is sufficiently smaller than $n$ (in particular,
at least when $\log_{n}(i)<0.31389$, see Huang1998 (38, 20)). In such cases,
where $i\ll n$, the cost of computing the LLEs scales as $\mathcal{O}(n^{2})$
as it is dominated by matrix multiplication rather than QR-decomposition. Use
cases where it suffices to know a subset of the LLE spectrum include
assimilation in the unstable subspace (e.g. see Carrassi2022 (17)) and
computing the local Kaplan-Yorke dimension, which can also be exploited for
better data assimilation Quinn2020 (65).

##  3 Using supervised machine learning to estimate
Lyapunov exponents

###  3.1 Problem statement and evaluation metrics

Supervised learning refers to ML algorithms that use data sets formed of
input-target pairs, whereby the goal is to construct a statistical model that
emulates the idealised function that maps from the input to the target. The
input and target are multidimensional arrays of data, not necessarily of the
same dimensions. A single input-target pair is known as an example; the size
of a ML data set refers to the number of examples it contains.

Supervised learning algorithms construct statistical models by optimising the
model’s parameters using the data. In the problem of this study, the input is
the system state at a set of consecutive recent time steps including the
current time $t$. The target is the vector containing the full spectrum of $n$
LLEs calculated using the method described in Section 2.2 by integrating
perturbations from time $t-\tau$ to $t$. We choose to estimate the full LLE
spectrum, however we note that the ML approaches we use can be easily adapted
to the sub-problem of estimating a subset of the LLE spectrum (henceforth,
“the sub-problem”), such as the unstable and near-neutral LLEs. See also
remarks in Section 2.2 and Section 5.3. Generally, we have

$\displaystyle(\text{input},\text{target})=\left(\
(\mathbf{x}_{k-r},\dots,\mathbf{x}_{k-1},\mathbf{x}_{k})\,,\quad(\lambda_{k}^{1},\dots,\lambda_{k}^{n})\
\right),$ (12)

where we recall that $\mathbf{x}_{k}\in{\mathbb{R}}^{n}$ denotes the system
state at time $k\Delta t$, and $\lambda_{k}^{i}$ is the $i^{\text{th}}$ LLE
computed from the interval $[k\Delta t-\tau,k\Delta t]$.

#### Point-wise accuracy

By point-wise accuracy we refer to the ability of a ML algorithm to predict a
specific LLE at an arbitrary time $t$. For its evaluation we calculate a
separate $R^{2}$ score for each LLE in the spectrum, from a set of $d$
predictions and targets. The $R^{2}$ score, also known as the coefficient of
determination, is given by

$\displaystyle R^{2}\Big{(}\,\left\\{\,(\ y_{j},\hat{y}_{j}\ )\mid
j=1,\dots,d\,\right\\}\,\Big{)}=1-\frac{\sum\limits_{j=1}^{d}(y_{j}-\hat{y}_{j})^{2}}{\sum\limits_{j=1}^{d}(y_{j}-\bar{y})^{2}}\in(-\infty,1],$
(13)

where, for each $j$, $y_{j}\in\mathbb{R}$ is the target output (e.g. the $i$th
LLE), $\hat{y}_{j}\in\mathbb{R}$ is the model’s prediction, and $\bar{y}$ is
mean of the target outputs. In Equation 13, the numerator is known as the sum
of squares of residuals, and the denominator is known as the total sum of
squares. An $R^{2}$ score of 1 is optimal, and an $R^{2}$ score of $0$ is as
good as guessing the mean of the target values every time.

#### Similarity of prediction and target distributions

In addition to the point-wise accuracy, we evaluate the statistical accuracy
of the ML models with quantile-quantile (QQ) plots. QQ plots provide a simple
non-parametric tool to compare the empirical probability distributions
generated by two samples Wilk1968 (90). In our case these are the predicted
and the target values. To generate the plot, a set of quantiles (the
$1000$-quantiles in our experiments) is computed for both samples. These
quantiles are then plotted against each other in a scatter plot. If the two
samples have the same empirical distributions, the scatter plot renders a 45
degree diagonal line (of course this is subject to sampling error, which
diminishes as the sample size grows). Departures from this ideal result show
differences in the location and scale parameters of the empirical
distributions, as well as possible linear and non-linear relationships between
the variables, see e.g. NISThandbook (52) for a more detailed discussion. In
our case, the QQ plots are useful to show which parts of the target
distribution are well represented by the predictions.

###  3.2 Supervised learning algorithms

We test four algorithms, summarised in Table 1, all well known in the ML
community. They are chosen to represent commonly used, proven-successful
supervised learning algorithms. In this section we detail the algorithms,
their structure and their relative characteristics. The final details of the
algorithms, including the number of parameters, are determined by
hyperparameter tuning and described in Section 5.2. We note that superior
performance in supervised learning tasks has been achieved by conducting
neural architecture search (NAS): an extensive (and costly) optimisation of
neural network (NN) architecture from a vast and highly flexible search space
Zoph2018 (94). Here we stop short of conducting such an NAS. Instead, we
choose established architectures for four distinct algorithms, and carry out
hyperparameter optimisation for each, where the hyperparameters include key
architectural choices such as the number of layers and the number of neurons
in each layer. We expect that the results from our selection give a good
indication of the possible performance of supervised learning in this task.

Algorithm | Architecture
---|---
Multilayer perceptron (MLP) | One or more dense layers and one dense output layer
Regression tree (RT) | One tree per target LLE
Convolutional neural network (CNN) | One 1D-convolution layer, max pooling layer, flatten layer, one or more dense layers
Long-short term memory network (LSTM) | One or more LSTM layers, one dense output layer

Table 1:  The four supervised learning algorithms used in this study. The
hyperparameter values (e.g. number of dense layers) used in the experiment are
described in Section 5.2

As we discuss in the following paragraphs, the chosen algorithms take
different approaches to using the temporal structure of the input (when there
are multiple time steps in the input). Here by temporal structure we mean the
temporal sequence of elements in each input vector, as opposed to the pairwise
relation between the inputs and outputs while segmenting the time-series for
data set preparation. The relative success of each algorithm gives insight
into the nature of the problem from the ML perspective.

The first algorithm is the regression tree (RT) breiman1984 (14). The RT is
the only non-neural network algorithm we test; RTs function by evaluating a
finite chain of comparisons on the input features, such as “$x_{k}>0$”. The
chain of comparisons form a tree graph; each leaf node corresponds to an
output value. Thus RTs have finitely many possible output values. The key
advantage of a RT is its low computational cost of making predictions. Other
advantages include the implicit feature selection process and potentially
greater explainability of predictions compared with NNs. In fact, the RT may
make use of only a few features from the set available in the input vector
breiman1984 (14). By setting two hyperparameters (the maximum number of leaf
nodes and the maximum depth, i.e. number of consecutive comparisons before a
leaf node), it is possible to greatly constrain the size of the resulting
tree. For simplicity, we opt to train a separate RT for each target LLE: i.e.,
the prediction of the target vector $(\lambda_{k}^{1},\dots,\lambda_{k}^{n})$
is made by $n$ RTs, where each RT predicts a different $\lambda_{k}^{i}$.

The second algorithm is the multilayer perceptron (MLP), the most basic type
of feedforward artificial NN, (as described in, e.g., Goodfellow-et-al-2016,
27, Chapter 6). The MLP is comprised of several hidden layers and an output
layer, each of which is comprised of many neurons, where each neuron receives
as input the outputs from all neurons in the previous layer (i.e., each layer
is densely/fully connected). Each hidden layer has the same number of neurons;
this value is optimised as a hyperparameter (see Section 5.2). When making a
prediction, the entire input vector is passed to every neuron in the first
hidden layer. These neurons compute their outputs according to their weights
and activation function, and their outputs are passed onto the neurons of the
second layer, etc., until the neurons in the final layer produce outputs which
are taken to be the final output of the algorithm. The MLP can theoretically
approximate any continuous function Hornik1991 (37), however in practise it
has been found that in many problems accurate predictions are more easily
achieved with more sophisticated architectures.

With regards to using temporal structure, the RT and the MLP take the same
approach: to treat each element of the input vector (each feature) as an
independent and identically distributed variable. Both algorithms have no
inherent preference with respect to the temporal structure of the data, since
they are invariant to the choice of order of the input features (a choice
which is made in the preprocessing stage before training the algorithm).

The third algorithm is a convolutional neural network (CNN) LeCun1990 (43). In
our experiments, the CNN is comprised of one 1D-convolution layer with a
kernel of size two and a stride of one, a max-pooling layer (if the number of
time steps in the input is greater than 2), a flattening layer and finally a
set of dense layers (the number of which is optimised as described in Section
5.2). The 1D-convolution layer is “1D” in that the layer convolves only over
the time dimension. The kernel size of two means that if the input is
comprised of system state vectors at $r+1$ consecutive time steps, the
convolution layer is only sensitive to patterns that occur in any time window
of length two within the $r+1$ time steps. The max-pooling layer has a pool
size and stride of two, which results in an invariance to translation of
patterns by a single time step. Thus the CNN takes a different approach to
using the temporal structure of the input, by only being sensitive to patterns
that appear in small time windows within the input. In other words, the
convolution layer predisposes the CNN to be aware which system states in the
input are adjacent to each other in time. CNNs have shown to be very
successful in a range of applications, including in image-related tasks
LeCun2015 (42). The CNN architecture in our experiments is equivalent to the
MLP with an additional 1D-convolution layer at the beginning. We choose a 1D
convolution (i.e. a convolution along the time dimension only), because in the
ODE systems in our experiments, there are only three state variables and there
is not a spatial locality which would motivate a focus on two of the variables
at a time.

The fourth algorithm is the long short-term memory network (LSTM)
Hochreiter1997 (35, 29). An LSTM is a form of recurrent neural network (RNN):
the NN processes data in a sequence, and after each term in the sequence is
processed, information is stored in a hidden state. In an LSTM, the flow of
information in and out of the NN’s hidden state is controlled by learned
gates. LSTMs have been successful in various tasks where there are long-term
dependencies, i.e. where the correct output at a later element in the sequence
requires information from elements further back in the sequence. Such tasks
include speech recognition Graves2013 (30) and machine translation
wu2016googletranslate (93). In machine translation, for instance, it is useful
to retain information about words early on in the sentence to best predict how
to translate words at the end of the sentence. However, the sophistication of
the LSTM architecture comes at a computational cost which we noticed
particularly during training. In the experiments of this study, the LSTMs are
composed of The LSTM layers process the input one time step at a time, in
chronological order. The estimation of the LLEs is made after the last time
step has been processed. Thus, the LSTM takes a third approach to using the
temporal structure by using information from past time steps when processing
future time steps.

The NN algorithms (MLP, CNN and LSTM) were implemented using Tensorflow tf
(50). Additionally, all NN algorithms standardise the input and target data by
subtracting the mean and dividing by the standard deviation, where the mean
and standard deviation are calculated from the training set. Each component of
the input and target is standardised independently. This is to make the data
more amenable to learning. The scaling pipeline was implemented using Scikit-
learn scikit-learn (63). For the RT, we use the implementation in the Scikit-
learn Python module scikit-learn (63).

##  4 Rössler, Lorenz 63 and their local Lyapunov spectra

In this section we present the two dynamical systems used in this work. We
discuss the characteristics of their attractors and their LLEs, as this forms
the data for the ML models and will be important in understanding their
performance. The Rössler Rossler1976 (72) and the Lorenz 63 Lorenz1963 (47)
systems are both three-dimensional, continuous time, ODE dynamical systems,
given respectively by Equations 4 and 4. We use the parameters
$(a,b,c)=(0.37,0.2,5.7)$ and $(\sigma,\rho,\beta)=(10,28,8/3)$. These are
commonly chosen values for which the systems exhibit chaotic behaviour (see
e.g. ott_2002, 57). Under these settings, both systems are dissipative, and
they possess strange attractors of fractal dimension.

Rössler $\displaystyle\openup-3.0pt\begin{split}&\dot{x}=-y-z\\\
&\dot{y}=x+ay\\\ &\dot{z}=b+z(x-c)\end{split}$ (14)

Lorenz 63 $\displaystyle\openup-3.0pt\begin{split}&\dot{x}=\sigma(y-x)\\\
&\dot{y}=x(\rho-z)-y\\\ &\dot{z}=xy-\beta z\end{split}$ (15)

The Lorenz 63 system is famous amongst weather and climate scientists:
initially derived as a truncated model of Rayleigh-Bérnard convection in a
two-dimensional fluid flow, it is the archetypal chaotic system and continues
to be used in weather and climate science, for example in data assimilation
experiments Carrassi2018 (16). For the chosen system parameters, the attractor
of the Lorenz 63 system is formed of two wings, each centred around a non-
stable fixed point. There is a third non-stable fixed point at the origin
Sparrow1982 (79). The Rössler system was introduced as a simpler version of
the Lorenz 63 system, having only one non-linear term ($zx$) instead of the
two in the Lorenz 63 system ($xz$ and $xy$). Equations 4 were derived as a
simplification of a system that combined two chemical reactions: a slow, two-
variable oscillator ($x$ and $y$) and a faster “switching-type” reaction ($z$)
Rossler1976 (72). For a historical review of the development of the Rössler
system, see Letellier2010 (45). The resulting system attractor is composed of
a “disc” in the $xy$ plane ($z$ close to $0$), and a “loop” in which
trajectories rise ($z$ increases rapidly) out of the disc, before folding over
and back into the disc. Consequently there is an imbalance: any infinite, non-
periodic trajectory spends more time in the disc than in the loop. The
attractor has two fixed points: an unstable fixed point in the centre of the
disc, around which system trajectories spiral outwards, and a stable fixed
point located outside of the attractor. We shall see that the difference in
the systems’ dynamics provides useful comparison.

It is easy to compute the analytic Jacobian matrices from Equations 4 and 4.
Furthermore, the small size of both systems allow us to perform exhaustive
experiments with long time series. The data used in our results are generated
following Section 2.2. The evolution equations, Equations 4 and 4, as well as
the corresponding fundamental matrix equation, Equation 4, are integrated
using a fourth-order Runge-Kutta scheme with timestep $\Delta t=0.01$ for both
systems. The LLEs are calculated over time windows of length $\tau=0.04$ (i.e.
four timesteps); the LLEs at time $t$ are calculated by integrating
perturbations from time $t$ to $t+\tau$. We choose $\tau=0.04$, rather than a
smaller value, so that the resulting set of LLEs provide a better coverage of
the attractor. Table 2 shows the resulting LEs. Further details of the data
used in the results are given in Section 5.1.

System | $\hat{\lambda}_{1}$ | $\hat{\lambda}_{2}$ | $\hat{\lambda}_{3}$
---|---|---|---
Rössler | $0.19597\,\pm\,0.00011$ | $0.0000075\,\pm\,0.0001275$ | $-4.30097\,\pm\,0.00068$
Lorenz 63 | $0.90495\,\pm\,0.000145$ | $0.001975\,\pm\,0.000055$ | $-14.571345\,\pm\,0.000105$
Table 2: The LEs of the Rössler and Lorenz 63 systems as calculated following
Section 2. These values are computed from 718,800 LLE iterations, following a
transient of $k=1,200$ iterations. As described in Section 2, the arithmetic
mean of the LLEs converges to the LEs as one includes LLEs from more
iterations, i.e., from a longer trajectory. The convergence is not monotonic:
the series of arithmetic means fluctuates as the number of iterations
increases. To give an indication of the precision of the numerical
calculation, we therefore calculate the LE as the mid-point of the range of
the series of arithmetic means acquired from the first $j$ iterations, where
$j=716,801,\dots,\allowbreak 718,800$. The extent of the range above and below
this value are also given. The proximity of the second LE to zero is a test
for the accuracy of the numerical algorithm, since the middle LE is
theoretically known to be zero in chaotic autonomous continuous-time systems
Pikovsky2016 (64).

In Figures 1 and 2, the top row shows the values of the LLEs along the
systems’ trajectories, specifically for the points $\mathbf{x}(\tau j)$ where
$1,200\leq j\leq 26,199$; the values of the LLEs are given in color. We show
only $25,000$ data points to avoid saturating the figures: the local
heterogeneity of the LLE values in phase space is evident. The bottom row
displays the distribution of LLE values via histograms, for the points
$\mathbf{x}(\tau j)$ where $1,200\leq j\leq 719,999$. The $i$th column shows
the LLEs associated with the $i$th LE.

Figure 1: A long-time trajectory of the Rössler system coloured by LLE values
(top row) shows how LLE values tend to be arranged in the system’s attractor.
The bottom row shows the corresponding statistical distribution of the LLEs
via histograms; the mean of the LLE values is shown by the dotted orange line.
The top panels show the same collection of $25,000$ points. On the other hand,
the histograms are generated from the full set of $718,800$ LLEs; note that
the vertical axis is plotted in logarithmic scale. The mean (i.e. the
corresponding LE) and standard deviation of the LLE values in each column are
shown underneath.

Figure 2: As with Figure 1 but for Lorenz 63.

Let us first consider the statistical distribution of the LLEs. In the Rössler
system (Figure 1) there is a marked difference between the first two LLE
distributions and the third. The first two have a single, tall, thin mode with
long tails and are roughly symmetric. The distribution of the third LLE has a
lower-density mode with thicker tails and is negatively skewed. The range of
values of the third distribution is also much smaller than that of the first
two.

The Lorenz 63 system LLE distributions differ from those of the Rössler
system. All three LLE distributions (Figure 2) have thicker tails and are
positively skewed. The first and third are weakly bimodal, whilst the second
is unimodal.

We now introduce a property that will be important to understanding the
performance of the ML. In the context of the distribution of LLEs in the
attractor of a system, we say that a region $U$ of the attractor $A$ is
(locally) heterogeneous if the function $f\colon A\to\mathbb{R}$ mapping from
$A$ to the $i$th LLE value is non-smooth in $U$. The best definition in more
simple terms is that the LLE values in $U$ vary in a non-smooth fashion across
$U$. For example, an alternating lattice (such as a chess board) would be
locally heterogeneous. If $U$ is not heterogeneous, we say it is homogeneous.
In this work we use homogeneous and heterogeneous only in the sense of locally
across a small region of phase space $U$, and not in the sense of across the
entire attractor (i.e. globally), or along a trajectory (i.e. across time), as
is the case in, e.g., Vannitsem2017 (85), and lucarini2020new (48).

For all three Rössler LLEs, in the disc of the attractor that sits in the
$xy$-plane, the LLE values are homogeneous. In the loop that jumps out of the
disc with positive $z$-values, there are bands of similar values for all three
LLEs. For the first two LLEs, the most extreme LLE values are in the loop, and
there is significant local heterogeneity within the bands. For example, there
are some small regions of the attractor where the first LLE is $10$ or above
for the majority of points, yet for some other points it is as low as $-20$.
This high degree of heterogeneity is reflected in the longer tails of the
distributions of the first and second LLEs (cf Figures 1 and 2). In contrast,
the third LLE is locally homogeneous: the colour change in the graph is
smooth. Note that the trajectory of the Rössler system spends more time in the
disc than in the loop. Consequently, data points in the loop are considerably
sparser.

In the Lorenz 63 system, the first two LLEs have distinct regions of local
heterogeneity and homogeneity. Unlike with the Rössler system, the regions of
greater mixing are along the boundaries between regions of greater local
homogeneity. The most extreme values of the first LLE are at the bottom of the
attractor ($z$ close to 0) and on the top edge of the wings: not in the
regions of local heterogeneity. The third LLE, as with the Rössler system, is
everywhere locally homogeneous. As we will see, the characteristics of the
LLEs we have discussed have implications for the performance of the ML methods
trying to estimate them.

##  5 Experiments and results

###  5.1 Experiment setup

Having described the ML algorithms (Section 3.2) and dynamical systems
(Section 4) used in this work, this section details the experiments and
results. The experiments test a total of $16$ configurations, which consist of
each combination of two dynamical systems, four ML algorithms, and two input
types:

Input Type 1: $\displaystyle(\mathbf{x}_{k})$ (16) Input Type 2:
$\displaystyle(\mathbf{x}_{k-5},\,\mathbf{x}_{k-4},\,\mathbf{x}_{k-3},\,\mathbf{x}_{k-2},\,\mathbf{x}_{k-1},\,\mathbf{x}_{k})$

The first input type consists of the current time step only, whereas the
second has the current time step and the five preceding time steps. In the
envisaged operational application of this approach, storing multiple time
steps of the entire model state poses a severe computational challenge. Since
we envisage making LLE predictions every time step, and we have assumed that
the pattern of input time steps remains fixed, the furthest-back input time
step dictates the number of time steps that must be stored. We choose $5$ time
steps into the past as a balance between testing whether previous time steps
can enable more accurate predictions, and not requiring huge amounts of steps
to be stored. Ignoring the constraints of feasibility and our approach, we
expect that delays of more than $5$ time steps might enable more accurate
predictions. Given the choice of a maximum $5$ time steps into the past, we
include all $6$ time steps in the input and rely on the ML algorithms to
extract useful features thereof. Note that if there is only one time step in
the input, the CNN is equivalent to an MLP where the number of units in the
first hidden layer of the MLP is twice the number of filters in the CNN.

The main results were attained with data sets of $10^{5}$ examples. The data
sets were created using the method described in Section 2.2, and the
parameters given in Section 4. For our chosen value of $\tau=0.04$, the
$10^{5}$ examples are generated from a trajectory of $4\times 10^{5}$ time
steps, equating to $784$ and $3604$ Lyapunov times for Rössler and Lorenz 63,
respectively. We note that another study which estimated the global LEs of the
Lorenz 63 system by emulating the dynamics with reservoir computers used a far
smaller data set of $91$ Lyapunov times Pathak2017 (62). By inspecting plots
of the trajectory (not shown here) we anticipated that $10^{5}$ examples
provide sufficient coverage of the system attractors. The Kullback-Leibler
divergences of input and target variables (not included here) show that
$10^{5}$ examples provides a good representation of the variables’ statistical
distributions (when compared to a data set of $10^{6}$). To assess the impact
of the data set size, we repeated experiments with data sets of $5\times
10^{4}$ and $10^{6}$ examples. We found that the much larger data set resulted
in only slightly more accurate predictions, with the exception of the LSTM.
Thus the initial choice of $10^{5}$ does not limit performance substantially.
These results are discussed further in Section 5.3.

From each system’s data set we generate $30$ data set instances, where each
instance is a unique random shuffle of the original. Each data set instance is
partitioned into training, validation and test sets with a ratio of
$0.6:0.2:0.2$. The resulting data setup is summarised in Table 3. The
training, validation and test sets of different data set instances are
therefore distinct. Given that the Rössler system and Lorenz 63 system are
both ergodic, and the size of the training, validation and test data sets are
large, the coverage of the attractors by examples in each data set instance is
similar to that which would be achieved if we had instead generated $30$ sets
of new data. Each configuration is tested on the $30$ data set instances, i.e.
each configuration is tested in $30$ trials. This provides an estimate of the
variability of performance of each algorithm.

Number of data set instances | $30$
---|---
Number of examples: |
in each data set instance | $100,000$
in the training partition | $60,000$
in the validation partition | $20,000$
in the test partition | $20,000$

Table 3: A summary of key values in the data setup for the $10^{5}$ data sets.
Note that each data set instance contains the same examples in a unique,
shuffled order.

The NNs (MLP, CNN, LSTM) are implemented with two methods for preventing
overfitting: activity regularisation on each layer, and early stopping.
Regularisation penalises large weight values; early stopping selects the model
weights that score optimally on the validation data set, rather than on the
training data set.

###  5.2 Hyperparameter optimisation

In ML, hyperparameter refers to any parameter that has to do with the form of,
or method of optimising, the statistical model, as opposed to the trainable
parameters of the model itself (often referred to as weights). Typically,
hyperparameters are fixed before the model is fit to the data, i.e. before the
weights are optimised. For example, one hyperparameter for a NN is the
learning rate: the amount by which model weights are changed at each step in
the optimisation. Selecting the right hyperparameters is essential for
effective use of ML algorithms Goodfellow-et-al-2016 (27).

We use a Bayesian optimisation method, implemented by Scikit-Optimize
scikit_optimize (33), to optimise the hyperparameters for each ML algorithm
used in this study. Conceptually, a Bayesian optimisation method is an
informed hyperparameter search that generates a probabilistic model (e.g.,
using a surrogate Gaussian process regression) of the true ML model (e.g.,
CNN) to select a set of hyperparameters that maximises the true ML model’s
performance Snoek2012 (77).

We perform separate hyperparameter optimisation for each configuration.
Hyperparameter optimisation was carried out on Google’s tensor processing
units (TPUs) using Google Colaboratory. The search domain for each
hyperparameter was chosen based on users’ knowledge of the algorithms, the
nature of the problem, and common practice in the ML community Hastie2009 (32,
27). The search domains are shown in Table 4. The chosen search space permits
the NNs to be reasonably large (up to $200$ neurons per layer and up to $10$
layers for MLP and CNN, and up to 100 LSTM units per layer and up to 3 LSTM
layers for LSTM) given the low-dimensionality of the task: mapping from three
($3$ variables $\times\ 1$ time steps), or eighteen ($3\times 6$) features to
three outputs. On the other hand, we opted to restrict the maximum number of
layers for the LSTM due to the greater complexity of the algorithm. In fact,
the first attempt to perform hyperparameter optimisation with a $6$-layer LSTM
exceeded the $24$ hours runtime limit for Google Colaboratory. In contrast to
the NNs, we forced the RTs to remain computationally light by limiting maximum
leaf nodes (i.e. the maximum number of possible output values) to $100$. The
entire hyperparameter optimisation process had a combined runtime of $108$
hours. Although the hyperparameter optimisation was computationally expensive
it was affordable in the low dimensional problems at hand. It greatly
increases the chance that we attain maximal performance from each ML
algorithm, thus providing useful insights into the problem from a ML
perspective.

Algo | Hyperparameter | Search domain
---|---|---
RT | maximum depth | $[1,\ 100]$
| maximum leaf nodes | $[5,\ 100]$
| maximum features |  $\\{$None, $\log_{2}$, square root$\\}$
| splitter |  $\\{$best, random$\\}$
| min. cost-complexity pruning parameter | $[1\times 10^{-6},\ 100]$
| minimum examples per leaf | $[1,\ 20]$
| minimum weight fraction per leaf | $\\{0,0.5\\}$
MLP | learning rate | $[1\times 10^{-6},\ 0.9]$
| number of layers | $[1,\ 10]$
| number of neurons per layer | $[1,\ 200]$
| activity regularisation on each layer |  $\\{L1(\alpha=0.001)$, $L2(\alpha=0.001)$, None$\\}$
CNN | learning rate | $[1\times 10^{-6},\ 0.9]$
| number of filters | $[1,\ 100]$
| number of dense layers | $[1,\ 10]$
| number of neurons per dense layer | $[10,\ 200]$
| activity regularisation on each layer |  $\\{L1(\alpha=0.001)$, $L2(\alpha=0.001)$, None$\\}$
LSTM | learning rate | $[1\times 10^{-6},\ 0.9]$
| number of LSTM layers | $[1,\ 3]$
| number of LSTM units per layer | $[1,\ 100]$
| activity regularisation on each layer |  $\\{L1(\alpha=0.001)$, $L2(\alpha=0.001)$, None$\\}$

Table 4: The search domain for hyperparameters of the ML algorithms. For the
MLP, CNN and the LSTM, the number of (dense) layers excludes the final
densely-connected output layer with three units.

|  |  | Optimised value:
---|---|---|---
|  |  | Rössler | Lorenz 63
Algo | Target LLEs | Hyperparameter | 1 time step | 6 time steps | 1 time step | 6 time steps
RT | LLE 1 | maximum depth | 78 | 58 | 100 | 100
|  | minimum examples per leaf | 20 | 13 | 15 | 16
| LLE 2 | maximum depth | 100 | 24 | 100 | 78
|  | minimum examples per leaf | 20 | 20 | 1 | 1
| LLE 3 | maximum depth | 17 | 48 | 100 | 72
|  | minimum examples per leaf | 2 | 1 | 14 | 20
MLP | all | learning rate | 7.300e-05 | 1.222e-04 | 7.628e-05 | 1.939e-05
|  | number of dense layers | 7 | 6 | 8 | 10
|  | number of neurons per dense layer | 200 | 182 | 165 | 200
|  | activity regularisation on each layer | L2 | L2 | L2 | L1
CNN | all | learning rate | 1.081e-04 | 1.586e-04 | 3.287e-04 | 1.016e-04
|  | number of filters | 100 | 37 | 29 | 59
|  | number of dense layers | 6 | 6 | 3 | 6
|  | number of neurons per dense layer | 200 | 200 | 122 | 200
|  | activity regularisation on each layer | None | L1 | L1 | L1
LSTM | all | learning rate | 2.857e-03 | 1.874e-03 | 2.722e-03 | 1.105e-04
|  | number of LSTM layers | 3 | 1 | 2 | 2
|  | number of LSTM units per layer | 62 | 100 | 100 | 100
|  | activity regularisation on each layer | L1 | L1 | L1 | None

Table 5: The optimal hyperparameter values selected by 50 iterations of
Bayesian optimisation. For all RTs, the optimal maximum leaf nodes was $100$,
and the optimal minimum weight fraction per leaf was $0$. Some RT
hyperparameters are excluded for brevity.

Table 5 shows the optimised hyperparameter values, resulting from $50$
iterations of the optimisation algorithm. Some optimal values are at the
boundary of the search domain (Table 4). This suggests that the algorithms
might have performed better with hyperparameter values beyond the chosen
search domain. For instance, in every experiment configuration, the optimised
value of the RT hyperparameter maximum leaf nodes is $100$, the maximum value
in the search domain. Other examples of hyperparameters whose optimised values
are at the boundary of the search domain for some configurations are the
number of neurons per layer for the MLP and the CNN, and the number of LSTM
units per layer. Nevertheless, the number of dense layers in the MLP and the
CNN, as well as the number of LSTM layers, are mostly not at the boundary.
Overall, from the results of the hyperparameter optimisation, we can argue
that although greater prediction accuracy might be achieved using NNs with
larger layers (more neurons or LSTM units), the number of layers (NN depth) in
the experiments is adequate.

To give an idea of the complexity of each of the four optimised ML algorithms,
Table 6 shows the size of the ML models, measured by the maximum number of
comparisons for the RT and the number of trainable parameters in the case of
the NNs. The model size is a function of the optimised hyperparameter values.
For the NNs with $6$ time steps as input, the model size of those that predict
the Lorenz 63 system LLEs is larger than those that predict the Rössler
system. This is perhaps reflective of the more chaotic dynamics of the Lorenz
63 system, which has a larger first LE as well as two non-linear terms in the
equations as opposed to only one in the Rössler system (see Equations 4 and
4). The same is not true of the NNs that take $1$ time step as input. This is
unsurprising since the $1$ time step input contains less information on the
dynamics, i.e., how the state variables are changing in time. In this case,
the size of the MLP and the CNN is smaller in the Lorenz 63 system, whereas
the size of the LSTM is larger in the Lorenz 63 system. Due to their
architectural similarity in the $1$ time step case, it is unsurprising that
the optimal model size for MLP and CNN behaves similarly.

|  | Rössler | Lorenz 63
---|---|---|---
Algo |  |  $1$ time step |  $6$ time steps |  $1$ time step |  $6$ time steps
RT | maximum comparisons | 99 | 99 | 99 | 99
MLP | trainable parameters | 242,603 | 170,537 | 192,888 | 366,203
CNN | trainable parameters | 222,503 | 224,262 | 34,244 | 237,616
LSTM | trainable parameters | 78,557 | 41,903 | 122,303 | 122,303

Table 6: The size of the ML models, measured by the number of comparisons (RT)
and number of trainable parameters (NNs). The maximum number of comparisons
for each RT is $99$ (one less than the maximum leaf nodes), since the max
depth is too large to constrain the number of non-leaf nodes in the tree.

The maximum size of the RTs is the same for all configurations. This is due to
the tight restriction placed on the tree size by the search domain of the
maximum leaf nodes hyperparameter (Table 4). However, the maximum depth
hyperparameter varies across configurations and between LLEs. The optimal
value of the maximum depth depends on the number of linear separations of the
input space that improves the predictions during the generation of the tree
from the training data. For the same input type, and for a given LLE, the RT
for the Lorenz 63 system has a greater maximum depth than the RT for the
Rössler system. The greater maximum depth implies that more accurate
predictions can be made by separating the input space at smaller scales (i.e.
making finer-grained partitions of the input space) in the Lorenz 63 system,
compared to the Rössler system. In other words, there is clearer detail at
smaller relative scales (in the input space) in the Lorenz 63 system than in
the Rössler system.

###  5.3 Results

####  Impact of data set size

Figure 3:  The impact of data set size on mean $R^{2}$ scores, across $30$
trials. The solid lines indicate $1$-time step results, the dashed lines
indicate $6$-time step results.

As discussed in Section 5.1, whilst we focus on results using a data set of
$10^{5}$ examples, we ran all experiments with three data set sizes to
determine the impact on prediction accuracy. We use the same setup for all
data set sizes: hyperparameters as described in Section 5.2, $30$ trials, and
a partition of $0.6:0.2:0.2$ training-validation-testing. The mean $R^{2}$
scores from all data set sizes are shown in Figure 3. In both systems, there
are only small (or negligible) gains in accuracy from using $10^{6}$ examples
compared to $10^{5}$ examples. The one exception is the LSTM: in both systems,
the $R^{2}$ scores of the LSTM significantly increased with data set size. In
the Rössler system, the LSTM becomes the most accurate method when trained on
$10^{6}$ examples, whilst in the Lorenz 63 system the LSTM achieves accuracy
similar to the MLP and CNN (indeed, Figure 3 shows the $6$-time step input
LSTM achieves the best mean $R^{2}$ scores). In both systems, the variation of
the LSTM over the 30 trials is substantially reduced with $10^{6}$ examples
compared to $10^{5}$ examples (notably, for the $1$-time step LSTM in Lorenz
63, the variance of $R^{2}$ scores of LLE 2 reduces from $0.0201$ to
$0.0007$). This suggests that choosing $10^{5}$ examples strongly limits the
performance of the LSTM. We suspect that the LSTM requires more data than the
MLP and CNN due to its more complicated architecture, namely the hidden state
and the three parameterised gates which control information flow into and out
of the hidden state.

With regards to the smaller data set of $5\times 10^{4}$, Figure 3 shows that
impact on accuracy is different in the two dynamical systems: in the Rössler
system, the accuracy of MLP and CNN (especially with the $1$-time step input)
is strongly reduced (compared to the $10^{5}$ data set), whereas for Lorenz 63
the equivalent reduction in accuracy is small. This is likely due to the
greater sparsity of data points in the loop of the Rössler system. With these
insights as context, the remainder of Section 5.3 refers to results from the
$10^{5}$ data sets, unless stated otherwise.

#### Comparisons between systems and across the LLE spectrum

We assess accuracy with the $R^{2}$ score of predictions made on the test data
sets, each of which has $20,000$ examples (see Table 3). There are $30$
$R^{2}$ scores for each configuration: one from each data set instance. These
$R^{2}$ scores are shown in boxplots in Figure 4 and summarised by their mean
and standard deviation in Table 7. The immediate observation is that the
$R^{2}$ scores differ consistently among LLEs (for a given system) and between
systems (for a given LLE). The first LLE is predicted at least reasonably well
in both dynamical systems ($0.54$ for Rössler system and $0.93$ for Lorenz 63
system). The third LLE is well predicted in both systems, by all ML
algorithms, and with both types of input. Apart from one case, the mean
$R^{2}$ scores for LLE $3$ are above $0.98$. This is unsurprising given the
local homogeneity of LLE $3$ on the attractors, as discussed in Section 4. In
all cases, the second LLE is the least well predicted ($0.39$ for Rössler
system and $0.76$ for Lorenz 63 system). This result is to be expected since
it is known that the second LLE is calculated least accurately by the
numerical method Kuptsov2012 (41) and has a slower convergence Bocquet2017
(8). Also, this result aligns with some recent attempts to emulate chaotic
dynamics with ML methods, where the emulators have failed to reproduce the
near-neutral LEs accurately Pathak2017 (62, 12). We note that, particularly in
multiscale systems such as ocean-atmosphere systems, the neutral and near-
neutral exponents play an important role in understanding predictability, (see
e.g. DeCruz2018 (21, 65)) and are connected to the coupling mechanisms
Vannitsem2016 (87, 82).

Next we compare prediction accuracy between dynamical systems. For LLEs $1$
and $2$, predictions of the Lorenz 63 system tend to be better than those of
the Rössler system. The highest mean $R^{2}$ score for LLEs $1$ and $2$ are
$0.9304$ and $0.7613$ (respectively) for Lorenz 63, yet only $0.5365$ and
$0.3897$ for Rössler. For LLE $3$, the mean $R^{2}$ scores are similarly high
in both systems. These results indicate that the LLEs can be predicted and the
variability of the prediction accuracy depends on which LLE and which
dynamical system is being predicted.

$R^{2}$ scores of test data sets for all 30 data set instances

Figure 4: The $R^{2}$ scores of test data sets from the $30$ trials, showing
the variation across data set instances, for each combination of system and ML
method. Perfect predictions have an $R^{2}$ score of one. These results use
the $10^{5}$ data sets:

each test data set has $20,000$ examples.

$R^{2}$ on test data: Mean and (standard deviation) of 30 data set instances
---
| Rössler
target | LLE 1 | LLE 2 | LLE 3
input type | 1 time step | 6 time steps | 1 time step | 6 time steps | 1 time step | 6 time steps
RT | 0.5155 (0.0248) | 0.5161 (0.0268) | 0.3506 (0.0299) | 0.3681 (0.0278) | 0.9944 (0.0002) | 0.9946 (0.0002)
MLP | 0.5323 (0.0211) |  0.5363 (0.0243) | 0.3837 (0.0249) |  0.3897 (0.0274) |  0.9978 (0.0005) | 0.9975 (0.0006)
CNN | 0.5279 (0.0333) | 0.5054 (0.0349) | 0.3711 (0.0518) | 0.3530 (0.0419) | 0.9960 (0.0040) | 0.9956 (0.0044)
LSTM | 0.4657 (0.0633) | 0.5319 (0.0462) | 0.2921 (0.0769) | 0.3788 (0.0571) | 0.9869 (0.0074) | 0.9955 (0.0023)
| Lorenz 63
target | LLE 1 | LLE 2 | LLE 3
input type | 1 time step | 6 time steps | 1 time step | 6 time steps | 1 time step | 6 time steps
RT | 0.8672 (0.0025) | 0.8801 (0.0027) | 0.6166 (0.0046) | 0.6540 (0.0053) | 0.9324 (0.0305) | 0.9936 (0.0003)
MLP | 0.9217 (0.0038) | 0.8350 (0.0438) | 0.7325 (0.0123) | 0.4449 (0.1446) |  0.9993 (0.0003) | 0.9925 (0.0033)
CNN | 0.9169 (0.0081) |  0.9304 (0.0047) | 0.7153 (0.0261) |  0.7613 (0.0159) | 0.9992 (0.0005) |  0.9993 (0.0002)
LSTM | 0.8659 (0.0530) | 0.7594 (0.0302) | 0.5702 (0.1419) | 0.2933 (0.0803) | 0.9955 (0.0044) | 0.9838 (0.0028)

Table 7:  The table shows mean $R^{2}$ scores over 30 trials, and the
corresponding standard deviations in parentheses. The $R^{2}$ score measures
the accuracy of predictions: the optimum score is 1. We calculate the $R^{2}$
on the test data set for each of the 30 trials. The highest mean $R^{2}$ score
for each combination of LLE and system (for both input types) is shown in
bold. “1 (6) time step(s)” refers to the number of time steps in the input.

#### Analysis of predictions on ordered test data

Figures 5 and 6 show time series of target values and predictions for a small
set of ordered test data. The predictions are produced by the algorithm that
achieves the best mean $R^{2}$ scores (on the $10^{5}$ data sets): MLP for the
Rössler system, and CNN for the Lorenz 63 system. In both systems, LLE 3 is
almost perfectly predicted throughout the time series. However, the
predictions of LLEs 1 and 2 have error characteristics that are specific to
each system.

Figure 5 illustrates that the first and second LLEs of the Rössler system vary
intermittently: they are stationary and near the mean value for the majority
of the time and then abruptly change and oscillate for a short period before
returning to be close to the mean. This corresponds to the system trajectory
being in the disc in the $xy$-plane, and then jumping into the “loop” with
positive $z$-values, before returning to the disc. Predictions are extremely
good in the stationary periods, and they are satisfactory during the peaks,
which we label “fluctuation events”. This is particularly true for LLE 1 where
we see that the ML-based predictions always catch the fluctuation event and
often its sign. The predictions of LLE 2 follows similar behaviour to LLE 1,
however the $R^{2}$ score suggests that the pointwise accuracy is slightly
worse than for LLE 1.

On the other hand, in the Lorenz 63 system, Figure 6 shows that LLEs 1 and 2
are constantly oscillating. Certain characteristics of the target time series
are well reproduced by the predictions, e.g. the largest peaks of LLE 1. These
large peaks occur when the system trajectory passes close to the origin (cf.
Figure 2), a region in which LLE 1 is locally homogeneous on the attractor.
Nonetheless, small errors occur frequently. Notably, the higher-frequency
features (such as the secondary peaks of LLE 2 between $t=17$ and $t=20$) are
often relatively poorly reproduced for LLEs 1 and 2.

These time series provide further insight into the lower $R^{2}$ scores for
LLEs 1 and 2 in the Rössler system compared to the Lorenz 63 system. Recall
the definition of $R^{2}$ in Equation 13: the distance from the perfect score
of $1$ is the sum of squared residuals divided by the total sum of squares.
The periods of stationarity in the Rössler system LLEs 1 and 2 contribute
little to the total sum of squares. Consequently the larger errors during
fluctuation events strongly reduce the $R^{2}$ score. In contrast, in the
Lorenz 63 system, the $R^{2}$ score is high despite more frequent prediction
errors because the constant variation of the target values results in a larger
total sum of squares.

Figure 5: Time series of targets and predictions of test data from the Rössler
system. Predictions made by an MLP with 6 input time steps. The $R^{2}$ scores
for the period shown are 0.4393, 0.3175 and 0.9981 for LLEs 1, 2 and 3,
respectively. Figure 6: As with Figure 5 for the Lorenz 63 system. The
predictions are made by a CNN with 6 input time steps. The $R^{2}$ scores for
the period shown are 0.7981, 0.4897 and 0.9953 for LLEs 1, 2 and 3,
respectively.

#### Impact of local heterogeneity in phase space

The variability of the LLEs on the strange attractor of chaotic systems is a
known feature, the immediate consequence of which is a highly state-dependent
predictability horizon: two slightly different initial conditions can generate
trajectories with hugely different degree of instability. In a recent work,
lucarini2020new (48) have for the first time shown how this variability is
related to the presence and distribution of unstable periodic orbits, each
with a different degree of instability, densely filling the attractor.
Arbitrary solutions are bounced among these unstable periodic orbits taking
their local instability features when they are in their proximity.

Recall from Section 4 that in both dynamical systems, there are regions of the
system’s attractor where the values of LLEs 1 and 2 are locally heterogeneous
(LLE 3 is everywhere locally homogeneous). The locally heterogeneous regions
in the Rössler system are in the loop with positive $z$-values, and in the
Lorenz 63 system they form a strip that lies halfway between the outside edge
and the centre of each wing. Figure 7 shows where the larger prediction errors
occur on the attractor, for all configurations with a $6$ time step input.
More precisely, it shows the detraction from the perfect $R^{2}$ score of 1
contributed by each point. We see that for all ML algorithms, larger errors
occur in the locally heterogeneous regions. Moreover, the locally
heterogeneous regions are robustly difficult: for the most accurate algorithm
in the Lorenz 63 system (CNN), larger prediction errors only occur in these
regions. This suggests that local heterogeneity plays a key role in
determining where on the attractor it is possible for ML to make reliably
accurate predictions of LLEs.

Figure 8 shows that a similar pattern occurs for the absolute relative error
of predictions. The relative error compares the size of the error to the size
of the target value. Notably, there are some large relative errors in the
locally homogeneous regions of the Rössler system attractor (i.e. in the disc)
since the relative error is especially punitive when the target value is close
to zero.

The difficulty of making accurate predictions in these regions is in line with
ML theory. ML algorithms work by optimising a model (such as a NN) to
approximate the map from the input to the target of the training data. ML
algorithms are successful if the optimised model also approximates the map
from input to target on unseen data, such as test data. This is possible only
if the training data provide enough information about the unseen data. The
fundamental problem in the locally heterogeneous regions is that the training
data cannot provide enough information because the LLE values are noisy. In
other words, the target values are highly variable even as length scales tend
to zero. Thus the target values of unseen data are likely to be quite
different to those of nearby seen data. Consequently, as local heterogeneity
increases, ML models are less able to generalise from training data to unseen
data.

The characteristics of the local heterogeneity explains the differences in
prediction accuracy between the two dynamical systems. In the Rössler system,
local heterogeneity in the loop of the attractor results in poorer predictions
during the aforementioned fluctuation events. As explained above, errors
during fluctuation events strongly reduce the $R^{2}$ score. On the other
hand, in the Lorenz 63 system the values of LLEs 1 and 2 in their respective
heterogeneous regions (see Figure 2) are relatively close to the mean: the
largest deviations from the mean are in locally homogeneous regions.
Therefore, the prediction errors resulting from locally heterogeneous regions
are likely to be small compared to the deviation of the target values from the
mean. Consequently, the differences in local heterogeneity explain the higher
$R^{2}$ scores achieved for LLEs 1 and 2 of the Rössler system, compared to
the Lorenz 63 system.

|
---|---
|
Figure 7: Recall the definition of the $R^{2}$ score (Equation 13): the distance from the perfect score of $1$ is given by the sum of squares of residuals divided by the total sum of squares. Here we give a local (in phase space) description of the contribution towards that distance made by each prediction in a set of $20,000$ test examples. Points are located at the current time step $\mathbf{x}_{k}$ of the input, and both coloured and sized by the square of the residual, divided by the total sum of squares. Darker points contribute a greater reduction to the $R^{2}$ score. If we denote the colour values as $a_{j}$, then $R^{2}=1-\sum a_{j}$. We show all configurations with a $6$ time step input. For each configuration, we use the data set instance for which the $R^{2}$ score was closest to the mean (as shown in Table 7). |
---|---
|
Figure 8: As with Figure 7, but point colour and size show the absolute
relative error, given by $|y-\hat{y}|/(\epsilon+|\hat{y}|)$, where
$\epsilon=10^{-6}$, $y$ is the prediction and $\hat{y}$ is the target. The
darker the point, the larger the absolute relative error of the prediction.

#### Impact of statistical distribution of targets and predictions

The poorer prediction accuracy for the Rössler system is also explained, to a
lesser extent, by the statistical distribution of the LLE values. As described
in Section 4, the Rössler system LLEs 1 and 2 include “extreme events”, i.e.
values of large magnitude that appear infrequently. Predicting these extreme
events is very challenging for any model, and particularly so for ML: one
would need to enlarge the training data set commensurately to the (very long)
return times of the extreme events. As noted above, we found a marked
improvement in the LSTM performance with a $10$-fold increase in data set size
(see Figure 3).

A complementary picture of the prediction accuracy is given in Figure 9 that
shows target values plotted against predicted values. The panels for the first
and second LLEs of the Rössler system show that all four ML algorithms fail to
predict the larger magnitude targets accurately. Similarly, the QQ-plots in
Figure 10 show that the predictions fail to replicate the extremities of the
true values: the minimum and maximum of the predictions are lower in magnitude
than those of the target. For instance, for LLE 1, we see that the NN methods
make few predictions greater than $10$ in magnitude, despite target values
reaching magnitudes of $25$. This behaviour can also be seen in the time
series plot (Figure 5). In the Lorenz 63 system, however, the larger-magnitude
targets occur more frequently and are thus well-represented in the training
data set. Consequently, we see in Figure 10 that the larger magnitudes are
well represented in the predictions distribution, even if the predictions
sometimes fail to capture the amplitude of the targets on a pointwise basis,
as illustrated in Figure 6.

Figure 9: Scatter plots of targets vs. predictions, for a test data set. Note
that the axis scales differ in each panel. For each method-system combination,
the data are from the data set instance with $R^{2}$ closest to the mean for
all data set instances. Only results from setups with 6 input time steps are
shown. Note that all scatter plots show different levels of
heteroskedasticity, i.e. the variance of a predicted value depends on the
value of the target. This is a well-known challenge for regression methods.
Figure 10: The $1000$-quantiles of the predictions are plotted against those
of the targets, revealing how well the two distributions match. The closer the
graph to the $y=x$ line (dotted white), the closer the prediction distribution
to the target distribution. Note that there is not necessarily any
relationship between proximity of distributions and (point-wise) accuracy of
predictions. In each panel, and for each machine learning method, the
quantiles are from the test data of the data set instance with $R^{2}$ closest
to the mean. Only results from setups with 6 input time steps are shown.

#### Impact of ML approaches to exploiting temporal structure

Each of the ML algorithms we test takes one of three approaches to exploiting
the temporal structure of the input, as discussed in Section 3.2. A comparison
of these approaches can only be made when there is non-trivial temporal
structure: this paragraph refers only to results with the $6$-time step input.
We find that there is no single optimal approach across both systems and all
data set sizes. Figure 4 shows that with the $10^{5}$ data set, the MLP and
LSTM perform comparatively well in the Rössler system but far more poorly in
the Lorenz 63 system. However, with the $10^{6}$ data set, clearer patterns
emerge: in both systems, the LSTM is the most accurate, and the MLP and CNN
perform similarly well. With the larger data set, there is not a clear
distinction in performance between the three approaches.

####  Impact of input type

With the $10^{5}$ data set, the input type has a big impact on the MLP and
LSTM algorithms, particularly for LLE $2$ of Lorenz 63: the mean $R^{2}$ score
for MLP is $65\%$ better with $1$-time step than with $6$, and for the LSTM it
is $94\%$ better. However, with the $10^{6}$ data set, clearer patterns
emerge: in the Rössler system, the $1$-time step input achieves the same or
better mean $R^{2}$ for all algorithms, whereas in Lorenz 63, the $6$-time
step input achieves higher mean $R^{2}$ for all algorithms, and lower variance
of $R^{2}$ for all NN algorithms. This may be due to the more complicated
dynamics of the Lorenz 63 system (see Section 4), or due to Lorenz 63 having a
smaller Lyapunov time ($1.1$ compared to $5.1$ in the Rössler system), meaning
that $6$ time steps provides more information in Lorenz 63. As mentioned in
Section 5.1, in an idealised case (without computational constraints) one
would optimise the input to include more distant time steps. The comparison
between the $10^{5}$ and $10^{6}$ data sets shows that some NN algorithms
require more data in order to fully exploit the additional information
provided in the $6$-time step input. With the $10^{6}$ data set, the maximal
difference of mean $R^{2}$ scores between input types (for any given
algorithm) is small: $1$-time step is $10\%$ better in Rössler (LSTM, LLE 2)
and $6$-time steps is $9\%$ in Lorenz 63 (LSTM, LLE 2).

Among two input types we tested, the $R^{2}$ scores in Table 7 (for the
$10^{5}$ data set) suggest that across all algorithms, there is only a small
advantage to be gained by providing $6$ time steps in the input rather than
$1$. However, for a given algorithm, the input type can have a big impact.
Figure 3 shows this is especially true for the smaller data set of $5\times
10^{4}$ in the Rössler system, where the NN’s prediction accuracy were far
lower with the $1$-time step input.

#### Computation time

As stated previously, this is a feasibility study focused on accuracy.
Nevertheless, we discuss here briefly the computational cost, keeping in mind
however that optimising the latter was not our priority. These arguments are
thus included for context and completeness, and not as a proof of viability.
Table 8 shows the mean (and standard deviation) time elapsed per prediction of
the three LLEs, computed over 4000 trials on a CPU processor in a personal
computer. The RT is two orders of magnitude faster than the NNs. This is
expected: the RT has far fewer trainable parameters (see Table 6) and, unlike
the NNs, does not require the evaluation of activation functions. The RT is
the only ML algorithm that is faster than the standard method of computing
LLEs (propagating perturbations and orthogonalising, ignoring the time taken
by the spin up iterations). Despite the comparative simplicity of the RT, it
achieves $R^{2}$ scores that are close to those of the NNs, although it is
less close in the Lorenz 63 system. It is likely that the number of output
values (constrained by the maximum leaf nodes hyperparameter) limited the RTs’
performance more strongly in the Lorenz 63 system, due to the greater
complexity of the Lorenz 63 attractor.

Execution time per set of 3 LLEs
---
|  | Standard Method | RT | MLP | CNN | LSTM
Rössler | mean | 2.32e-4 | 9.96e-5 | 2.37e-2 | 2.38e-2 | 2.34e-2
| $\sigma$ | 2.16e-5 | 3.51e-6 | 2.36e-2 | 2.57e-3 | 2.51e-2
Lorenz 63 | mean | 9.69e-4 | 9.88e-5 | 2.36e-2 | 2.37e-2 | 2.31e-1
| $\sigma$ | 2.53e-3 | 2.26e-6 | 2.85e-2 | 2.83e-2 | 1.60e-3

Table 8: The time elapsed for one set of LLEs to be calculated (by the
standard method) or predicted (by the ML algorithms). The table shows the mean
and standard deviation of elapsed time from 4000 trials.

Finally, we speculate on the cost of making predictions in the sub-problem
(see Section 3.1). The energy cost (measured in flops) of the RT scales
linearly with the number of LLEs, since a separate tree is trained for each
scalar target. However, the time cost remains the same as the RTs can be
executed in parallel. For the NNs, the potential cost saving is not clear
without further experimentation. Less expensive NNs could be attained via
pruning or distilling methods, e.g. see Molchanov2017 (51).

##  6 Discussion and summary

This study discusses the use of supervised machine learning (ML) to support
numerical forecasting of chaotic dynamics. A huge amount of work has appeared
recently at the crossroads between ML and the geosciences, whereby the former
has provided novel data-driven solutions to complement or, in some ideal
scenarios, to substitute the physical models (see e.g. Sonnewald2021, 78). In
this work we took a different approach that we referred to as “non intrusive”.

We did not pursue improving the given physical model with a ML model, but
rather using ML models as a supplementary tool that provides information to
drive adaptive decisions while running the prediction. The range of possible
desirable information is ample, e.g. anticipating a regime change or the onset
of intense convective events, as is that of consequent actions. In this work,
we focused on chaotic systems where real time knowledge of the unstable
properties of the system’s state is of paramount relevance. Local Lyapunov
exponents (LLEs) provide this knowledge in the form of the local (in time)
exponential rates at which errors about the system’s state evolve
Benettin1980theory (5, 64). Nevertheless, they are notably difficult to
compute, require the coding and maintaining of a tangent linear model, and the
computational cost grows fast with the system’s size. This work is a
feasibility study that investigates the accuracy with which supervised ML can
estimate the LLEs of a dynamical system trajectory based only on the system
state at the current time step and a few recent time steps.

We tested four supervised ML algorithms (a regression tree (RT), a multilayer
perceptron (MLP), a convolutional neural network (CNN) and a long short-term
memory network (LSTM)) on two dynamical systems (the Rössler system and Lorenz
63 systems). The dynamical systems are chaotic, dissipative, 3-variable ODE
systems. The algorithms encompass three approaches to exploiting the temporal
structure of the input.

Our results indicate that the best algorithm depends on the dynamical system,
the size of the data set, and on the number of time steps included in the
input. Overall, the results show that in certain conditions, the LLEs can be
predicted well: this depends on the system dynamics, the LLE being predicted,
and importantly on the local heterogeneity of the LLE in the proximity of the
given state. In particular, the average accuracy was lowest for the neutral
LLE. Further work is required to see if this result also holds in ocean-
atmosphere systems (and multiscale systems more generally), where the neutral
and near-neutral exponents are key to determining local predictability
DeCruz2018 (21, 65). Our results suggest that the feasibility of using
supervised ML to drive adaptive actions in an operational setting will depend
on the specific use case: the forecasting model, the desired target
information, and the intended adaptive actions.

Additionally, we investigated the impact of the size of the data set used to
train the ML algorithms. We found that with data sets of $10^{6}$ examples,
compared to $10^{5}$ examples, variance of $R^{2}$ score reduced but there
were only marginal improvements in the mean $R^{2}$ score. With the $10^{6}$
data sets, the LSTM performed best in both systems. However, with the $10^{5}$
data set, the LSTM was limited: the MLP performed best in the Rössler system
whereas the CNN performed best in the Lorenz 63 system. The RT achieves
accuracy that is close to the best-performing algorithm in both systems,
whilst being computationally much cheaper than the NNs. We tested two input
types: one with $1$ time step, and one with $6$ time steps (of the system
state). We found the best input type depends on data set size and dynamical
system: in the Rössler system, $6$ input time steps is better for the smallest
data set ($5\times 10^{4}$), whilst in the Lorenz 63 system, much more data
($10^{6}$) is required for the LSTM and MLP to achieve comparable performance
with the $6$-time step input. We further show that large prediction errors
occur when the current state is in a region of local heterogeneity on the
system attractor. Outside of the locally heterogeneous regions, the best-
performing algorithms make consistently accurate predictions. The differences
in local heterogeneity between the two systems explain the lower $R^{2}$
scores achieved for LLEs 1 and 2 in the Rössler system, compared to the Lorenz
63 system. We explain that local heterogeneity is an insurmountable problem
for deterministic ML predictions. This challenge could be mitigated if the ML
prediction also included a reliable uncertainty quantification. We suspect
that an uncertainty quantification could be made either by using Bayesian NNs
Wang2016 (89) or by including a measurement of the nearby local heterogeneity
in the target of each example.

The low-dimensional setting permitted extensive experimentation in this work,
providing lessons that will be useful should this “non-intrusive” approach be
taken in weather and climate prediction. The next steps will be to apply the
approach of this work to spatially-extended models with more dimensions. There
are several foreseeable challenges on the path from the very low-order models
of this work to the envisioned setting of operational weather prediction
models. The first challenge is to generate suitable data sets, since the
calculation of LLEs does not scale well and requires a tangent linear model
(see Section 2.2). However, the requirement of a tangent linear model can be
avoided by using bred vectors (e.g. Toth1997, 83, 84). Additionally, the
attractor of any numerical weather prediction (NWP) model is complex and high-
dimensional: very large data sets will be required if sufficient attractor
coverage is to be obtained. For a NWP model with $\mathcal{O}(10^{8})$
variables, one would need $\mathcal{O}(10^{12})$ input-target pairs to obtain
a ratio between the degrees of freedom of the NWP model and the size of the ML
data set that is similar to the ratio used in this study. If one were to use
ERA5 reanalysis data era5 (34) as input, a data set of $\mathcal{O}(10^{12})$
single time step inputs would amount to approximately $\mathcal{O}(10^{9})$ TB
of data, which is unfeasibly large. Furthermore, due to the long time scales
involved in teleconnection events, the number of such events can be small even
in long time-series. Given the number of input features, this can lead to a
“small data problem”, see e.g. Vecchi2022 (88) and references therein.

Therefore, it may be necessary to generate training data using a reduced-
dimension version of the operational model. For example, in Quinn2021 (66,
67), LLEs are computed by reducing the data dimension (via empirical
orthogonal functions) and constructing a multi-state vector auto-regressive
model.

Once initial training data has been generated, the cost of making predictions
with the ML model can be reduced by further reducing the dimension of the
training data, i.e. by performing feature-extraction Guyon2006 (31). For
example, ML techniques such as autoencoders can be used for dimension
reduction, e.g. see Mack2020 (49). Also, it may be possible to curate training
data sets strategically to reduce their size. Finally, if the intended use-
case requires only part of the LLE spectrum, then cost savings can be made a)
when generating training data, which scales as $\mathcal{O}(n^{2})$ rather
than $\mathcal{O}(n^{3})$ (see Section 2.2), and b) when making predictions
(see Section 5.3).

The computational benefit of the ML approach investigated here is two-fold.
The ML approach estimates LLEs directly from the current system state, thus
avoiding the cost of the long spin up that is required by the conventional
method for calculating LLEs. Second, the ML approach has the potential to be
cheaper per iteration of LLEs. We found that the lightest algorithm we tested,
the RT, was computationally cheaper (by a factor of 10) than the conventional
method for calculating LLEs (see Table 8).

Although the NNs were comparatively costly in this setting, we expect that in
a higher-dimensional, operational setting, NNs may be competitive. It is
unknown how the required NN size will increase with the system dimension: this
will require experimentation. The time-cost of making predictions with NNs may
be reduced (relative to the size of the NNs) by using purpose-built ML
hardware. On the other hand, the cost of numerically calculating LLEs (by
propagating perturbations and orthogonalising) will scale as
$\mathcal{O}({n^{3}})$ if computing the full spectrum, or
$\mathcal{O}({n^{2}})$ if the number of LLEs computed is much smaller than the
dimension of the NWP model $n$.

## Acknowledgements

DA is funded by a studentship from the Engineering and Physical Sciences
Research Council (EP/N509723/1). JA and AC acknowledge the support of the UK
National Centre for Earth Observation (grant no. NCEO02004). AC is also
supported by the project SASIP funded by Schmidt Futures – a philanthropic
initiative that seeks to improve societal outcomes through the development of
emerging science and technologies.

## References

  * (1) Aishah Albarakati, Marko Budišić, Rose Crocker, Juniper Glass-Klaiber, Sarah Iams, John Maclean, Noah Marshall, Colin Roberts and Erik S. Van Vleck  “Model and data reduction for data assimilation: Particle filters employing projected forecasts and data with application to a shallow water model”  In _Computers & Mathematics with Applications_ Elsevier, 2021  DOI: 10.1016/j.camwa.2021.05.026
  * (2) Martin Albrecht, Gregory Bard and William Hart  “Algorithm 898: Efficient Multiplication of Dense Matrices over GF(2)”  In _ACM Trans. Math. Softw._ 37.1 New York, NY, USA: Association for Computing Machinery, 2010  DOI: 10.1145/1644001.1644010
  * (3) Josh Alman and Virginia Vassilevska Williams  “A Refined Laser Method and Faster Matrix Multiplication”  In _Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA)_ Philadelphia, PA: Society for IndustrialApplied Mathematics, 2021, pp. 522–539  DOI: 10.1137/1.9781611976465.32
  * (4) Peter Arbenz  “The QR Algorithm. In “Lecture notes on Numerical Methods for Solving Large Scale Eigenvalue Problems”, pages 63-90, 2016.”  In _Numerical Methods for Solving Large Scale Eigenvalue Problems_ , Available at http://people.inf.ethz.ch/arbenz/ewp/lnotes.html, accessed 2021/06/20, 2016, pp. 63–90
  * (5) Giancarlo Benettin, Luigi Galgani, Antonio Giorgilli and Jean-Marie Strelcyn  “Lyapunov Characteristic Exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. Part 1: Theory”  In _Meccanica_ 15, 1980, pp. 9–20  DOI: 10.1007/BF02128236
  * (6) Giancarlo Benettin, Luigi Galgani, Antonio Giorgilli and Jean-Marie Strelcyn  “Lyapunov Characteristic Exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. Part 2: Numerical application”  In _Meccanica_ 15, 1980, pp. 21–30  DOI: 10.1007/BF02128237
  * (7) Christopher M. Bishop  “Neural networks for pattern recognition”  Oxford: Clarendon Press, 1995
  * (8) Marc Bocquet, Karthik S. Gurumoorthy, Amit Apte, Alberto Carrassi, Colin Grudzien and Christopher K… Jones  “Degenerate Kalman Filter Error Covariances and Their Convergence onto the Unstable Subspace”  In _SIAM/ASA Journal on Uncertainty Quantification_ 5.1, 2017, pp. 304–333  DOI: 10.1137/16M1068712
  * (9) Thomas Bolton and Laure Zanna  “Applications of Deep Learning to Ocean Data Inference and Subgrid Parameterization”  In _Journal of Advances in Modeling Earth Systems_ 11.1, 2019, pp. 376–399  DOI: 10.1029/2018MS001472
  * (10) Massimo Bonavita and Patrick Laloyaux  “Machine Learning for Model Error Inference and Correction”  In _Journal of Advances in Modeling Earth Systems_ 12.12, 2020  DOI: 10.1029/2020MS002232
  * (11) Stephen Boyd and Lieven Vandenberghe  “Introduction to Applied Linear Algebra: Vectors, Matrices, and Least Squares”  Cambridge University Press, 2018  DOI: 10.1017/9781108583664
  * (12) Julien Brajard, Alberto Carrassi, Marc Bocquet and Laurent Bertino  “Combining data assimilation and machine learning to emulate a dynamical model from sparse and noisy observations: a case study with the Lorenz 96 model”  In _Journal of Computational Science_ 44 Elsevier, 2020, pp. 101171
  * (13) Julien Brajard, Alberto Carrassi, Marc Bocquet and Laurent Bertino  “Combining data assimilation and machine learning to infer unresolved scale parametrization”  In _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 379.2194, 2021, pp. 20200086  DOI: 10.1098/rsta.2020.0086
  * (14) L. Breiman, J. Friedman, C.J. Stone and R.A. Olshen  “Classification and Regression Trees”  Taylor & Francis, 1984  DOI: 10.1201/9781315139470
  * (15) Roberto Buizza  “Introduction to the special issue on “25 years of ensemble forecasting””  In _Quarterly Journal of the Royal Meteorological Society_ 145 Wiley Online Library, 2019, pp. 1–11
  * (16) Alberto Carrassi, Marc Bocquet, Laurent Bertino and Geir Evensen  “Data assimilation in the geosciences: An overview of methods, issues, and perspectives”  In _WIREs Climate Change_ 9.5, 2018, pp. 1–79  DOI: 10.1002/wcc.535
  * (17) Alberto Carrassi, Marc Bocquet, Jonathan Demaeyer, Colin Grudzien, Patrick Raanes and Stéphane Vannitsem  “Data Assimilation for Chaotic Dynamics”  In _Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications (Vol. IV)_ Cham: Springer International Publishing, 2022, pp. 1–42  DOI: 10.1007/978-3-030-77722-7_1
  * (18) Matthew Chantry, Hannah Christensen, Peter Dueben and Tim Palmer  “Opportunities and challenges for machine learning in weather and climate modelling: hard, medium and soft AI”  In _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 379.2194, 2021, pp. 20200083  DOI: 10.1098/rsta.2020.0083
  * (19) Y. Chen, A. Carrassi and V. Lucarini  “Inferring the instability of a dynamical system from the skill of data assimilation exercises”  In _Nonlinear Processes in Geophysics_ 28.4, 2021, pp. 633–649  DOI: 10.5194/npg-28-633-2021
  * (20) Matthias Christandl, François Le Gall, Vladimir Lysikov and Jeroen Zuiddam  “Barriers for rectangular matrix multiplication”  arXiv, 2020  DOI: 10.48550/ARXIV.2003.03019
  * (21) L. De Cruz, S. Schubert, J. Demaeyer, V. Lucarini and S. Vannitsem  “Exploring the Lyapunov instability properties of high-dimensional atmospheric and climate models”  In _Nonlinear Processes in Geophysics_ 25.2, 2018, pp. 387–412  DOI: 10.5194/npg-25-387-2018
  * (22) Peter Düben, Umberto Modigliani, Alan Geer, Stephan Siemen, Florian Pappenberger, Peter Bauer, Andy Brown, Martin Palkovic, Baudouin Raoult, Nils Wedi and Vasileios Baousis  “Machine learning at ECMWF: A roadmap for the next 10 years”, 2021  DOI: 10.21957/ge7ckgm
  * (23) Ronan Fablet, Said Ouala and Cedric Herzet  “Bilinear Residual Neural Network for the Identification and Forecasting of Geophysical Dynamics”  In _2018 26th European Signal Processing Conference (EUSIPCO)_ IEEE, 2018, pp. 1477–1481  DOI: 10.23919/EUSIPCO.2018.8553492
  * (24) Gary Froyland, Thorsten Hüls, Gary P. Morriss and Thomas M. Watson  “Computing covariant Lyapunov vectors, Oseledets vectors, and dichotomy projectors: A comparative numerical study”  In _Physica D: Nonlinear Phenomena_ 247.1, 2013, pp. 18–39  DOI: 10.1016/j.physd.2012.12.005
  * (25) F. Ginelli, P. Poggi, A. Turchi, H. Chaté, R. Livi and A. Politi  “Characterizing Dynamics with Covariant Lyapunov Vectors”  In _Phys. Rev. Lett._ 99 American Physical Society, 2007, pp. 130601  DOI: 10.1103/PhysRevLett.99.130601
  * (26) Gene H. Golub and Charles F. Van Loan  “Matrix Computations”  Baltimore, MD: Johns Hopkins University Press, 2013  URL: https://jhupbooks.press.jhu.edu/title/matrix-computations
  * (27) Ian Goodfellow, Yoshua Bengio and Aaron Courville  “Deep Learning” http://www.deeplearningbook.org MIT Press, 2016
  * (28) Georg A. Gottwald and Sebastian Reich  “Supervised learning from noisy observations: Combining machine-learning techniques with data assimilation”  In _Physica D: Nonlinear Phenomena_ 423, 2021, pp. 132911  DOI: 10.1016/j.physd.2021.132911
  * (29) Alex Graves  “Long Short-Term Memory”  In _Supervised Sequence Labelling with Recurrent Neural Networks. Studies in Computational Intelligence, vol 385._ Springer Berlin Heidelberg, 2012, pp. 37–45  DOI: 10.1007/978-3-642-24797-2_4
  * (30) Alex Graves, Abdel-Rahman Mohamed and Geoffrey Hinton  “Speech recognition with deep recurrent neural networks”  In _2013 IEEE International Conference on Acoustics, Speech and Signal Processing_ , 2013, pp. 6645–6649  DOI: 10.1109/ICASSP.2013.6638947
  * (31) Isabelle Guyon, Steve R. Gunn, Masoud Nikravesh and Lotfi A. Zadeh  “Feature Extraction”  In _Feature Extraction_ Springer-Verlag Berlin Heidelberg, 2006  DOI: 10.1007/978-3-540-35488-8
  * (32) Trevor Hastie, Robert Tibshirani and Jerome Friedman  “The Elements of Statistical Learning”  Springer, New York, NY, 2009  DOI: 10.1007/978-0-387-84858-7
  * (33) Tim Head, Manoj Kumar, Holger Nahrstaedt, Gilles Louppe and Iaroslav Shcherbatyi  “scikit-optimize/scikit-optimize”  Zenodo, 2021  DOI: 10.5281/zenodo.5565057
  * (34) Hans Hersbach, Bill Bell, Paul Berrisford, Shoji Hirahara, András Horányi, Joaquín Muñoz-Sabater, Julien Nicolas, Carole Peubey, Raluca Radu, Dinand Schepers, Adrian Simmons, Cornel Soci, Saleh Abdalla, Xavier Abellan, Gianpaolo Balsamo, Peter Bechtold, Gionata Biavati, Jean Bidlot, Massimo Bonavita, Giovanna De Chiara, Per Dahlgren, Dick Dee, Michail Diamantakis, Rossana Dragani, Johannes Flemming, Richard Forbes, Manuel Fuentes, Alan Geer, Leo Haimberger, Sean Healy, Robin J. Hogan, Elías Hólm, Marta Janisková, Sarah Keeley, Patrick Laloyaux, Philippe Lopez, Cristina Lupu, Gabor Radnoti, Patricia Rosnay, Iryna Rozum, Freja Vamborg, Sebastien Villaume and Jean-Noël Thépaut  “The ERA5 global reanalysis”  In _Quarterly Journal of the Royal Meteorological Society_ 146.730, 2020, pp. 1999–2049  DOI: 10.1002/qj.3803
  * (35) Sepp Hochreiter and Jürgen Schmidhuber  “Long Short-Term Memory”  In _Neural Computation_ 9.8, 1997, pp. 1735–1780  DOI: 10.1162/neco.1997.9.8.1735
  * (36) James R. Holton and Gregory J. Hakim  “An Introduction to Dynamic Meteorology”  Elsevier, 2013  DOI: 10.1016/C2009-0-63394-8
  * (37) Kurt Hornik  “Approximation capabilities of multilayer feedforward networks”  In _Neural Networks_ 4.2, 1991, pp. 251–257  DOI: 10.1016/0893-6080(91)90009-T
  * (38) Xiaohan Huang and Victor Y. Pan  “Fast Rectangular Matrix Multiplication and Applications”  In _Journal of Complexity_ 14.2, 1998, pp. 257–299  DOI: 10.1006/jcom.1998.0476
  * (39) Eugenia Kalnay  “Atmospheric Modeling, Data Assimilation and Predictability”  Cambridge University Press, 2002  DOI: 10.1017/CBO9780511802270
  * (40) James L. Kaplan and James A. Yorke  “Chaotic behavior of multidimensional difference equations”  In _Functional Differential Equations and Approximation of Fixed Points_ Berlin, Heidelberg: Springer Berlin Heidelberg, 1979, pp. 204–227
  * (41) Pavel V. Kuptsov and Ulrich Parlitz  “Theory and Computation of Covariant Lyapunov Vectors”  In _Journal of Nonlinear Science_ 22.5, 2012, pp. 727–762  DOI: 10.1007/s00332-012-9126-5
  * (42) Yann LeCun, Yoshua Bengio and Geoffrey Hinton  “Deep learning”  In _Nature_ 521.7553, 2015, pp. 436–444  DOI: 10.1038/nature14539
  * (43) Yann LeCun, Bernhard Boser, John Denker, Donnie Henderson, R. Howard, Wayne Hubbard and Lawrence Jackel  “Handwritten Digit Recognition with a Back-Propagation Network”  In _Advances in Neural Information Processing Systems_ 2 Morgan-Kaufmann, 1990  URL: https://proceedings.neurips.cc/paper/1989/file/53c3bce66e43be4f209556518c2fcb54-Paper.pdf
  * (44) B Legras and R Vautard  “A Guide To Liapunov Vectors”  In _Seminar on Predictability_ 1 Reading, UK: ECMWF, 1996, pp. 135–146  URL: https://www.ecmwf.int/en/learning/workshops-and-seminars/past-workshops/1995-annual-seminar
  * (45) Christophe Letellier and Valérie Messager  “Influences on Otto E. Rössler’s Earliest Paper on Chaos”  In _International Journal of Bifurcation and Chaos_ 20.11, 2010, pp. 3585–3616  DOI: 10.1142/S0218127410027854
  * (46) James Lighthill, J M T Thompson, A K Sen, A G M Last, D T Tritton and P Mathias  “The Recently Recognized Failure of Predictability in Newtonian Dynamics [and Discussion]”  In _Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences_ 407.1832 The Royal Society, 1986, pp. 35–50  URL: http://www.jstor.org/stable/2397780
  * (47) Edward N. Lorenz  “Deterministic Nonperiodic Flow”  In _Journal of the Atmospheric Sciences_ 20.2, 1963, pp. 130–141  DOI: 10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2
  * (48) Valerio Lucarini and Andrey Gritsun  “A new mathematical framework for atmospheric blocking events”  In _Climate Dynamics_ 54.1 Springer, 2020, pp. 575–598
  * (49) Julian Mack, Rossella Arcucci, Miguel Molina-Solana and Yi-Ke Guo  “Attention-based Convolutional Autoencoders for 3D-Variational Data Assimilation”  In _Computer Methods in Applied Mechanics and Engineering_ 372, 2020, pp. 113291  DOI: 10.1016/j.cma.2020.113291
  * (50) Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S., Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu and Xiaoqiang Zheng  “TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems” Software available from tensorflow.org, 2015  URL: https://www.tensorflow.org/
  * (51) Pavlo Molchanov, Stephen Tyree, Tero Karras, Timo Aila and Jan Kautz  “Pruning Convolutional Neural Networks for Resource Efficient Inference”  In _International Conference on Learning Representations_ , 2017  URL: https://openreview.net/forum?id=SJGCiw5gl
  * (52) National Institute of Standards and Technology (U.S.)  “NIST/SEMATECH e-Handbook of Statistical Methods” [Online; accessed 21-December-2021], 2012  DOI: 10.18434/M32189
  * (53) Duong Nguyen, Said Ouala, Lucas Drumetz and Ronan Fablet  “EM-like Learning Chaotic Dynamics from Noisy and Partial Observations”, 2019  arXiv:1903.10335 [cs.LG]
  * (54) Duong Nguyen, Said Ouala, Lucas Drumetz and Ronan Fablet  “Variational Deep Learning for the Identification and Reconstruction of Chaotic and Stochastic Dynamical Systems from Noisy and Partial Observations”, 2021  arXiv:2009.02296 [cs.LG]
  * (55) Paul A. O’Gorman and John G. Dwyer  “Using Machine Learning to Parameterize Moist Convection: Potential for Modeling of Climate, Climate Change, and Extreme Events”  In _Journal of Advances in Modeling Earth Systems_ 10.10, 2018, pp. 2548–2563  DOI: 10.1029/2018MS001351
  * (56) Valery Oseledets  “A multiplicative ergodic theorem. Lyapunov characteristic numbers for dynamical systems”  In _Trans. Moscow Math. Soc._ 19, 1968, pp. 197–231
  * (57) Edward Ott  “Chaos in Dynamical Systems”  Cambridge University Press, 2002  DOI: 10.1017/CBO9780511803260
  * (58) Luigi Palatella, Alberto Carrassi and Anna Trevisan  “Lyapunov vectors and assimilation in the unstable subspace: theory and applications”  In _Journal of Physics A: Mathematical and Theoretical_ 46.25 IOP Publishing, 2013, pp. 254020
  * (59) T.. Palmer  “Predictability of the Atmosphere and Oceans: From Days to Decades”  In _Decadal Climate Variability_ Berlin, Heidelberg: Springer Berlin Heidelberg, 1996, pp. 83–155
  * (60) Dhruvit Patel, Daniel Canaday, Michelle Girvan, Andrew Pomerance and Edward Ott  “Using machine learning to predict statistical properties of non-stationary dynamical processes: System climate,regime transitions, and the effect of stochasticity”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 31.3 AIP Publishing LLC, 2021, pp. 033149  DOI: 10.1063/5.0042598
  * (61) Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu and Edward Ott  “Model-Free Prediction of Large Spatiotemporally Chaotic Systems from Data: A Reservoir Computing Approach”  In _Physical Review Letters_ 120.2, 2018, pp. 024102  DOI: 10.1103/PhysRevLett.120.024102
  * (62) Jaideep Pathak, Zhixin Lu, Brian R. Hunt, Michelle Girvan and Edward Ott  “Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 27.12 AIP Publishing, 2017, pp. 121102  DOI: 10.1063/1.5010300
  * (63) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot and E. Duchesnay  “Scikit-learn: Machine Learning in Python”  In _Journal of Machine Learning Research_ 12, 2011, pp. 2825–2830
  * (64) Arkady Pikovsky and Antonio Politi  “Lyapunov Exponents”  Cambridge: Cambridge University Press, 2016  DOI: 10.1017/CBO9781139343473
  * (65) C. Quinn, T.. O’Kane and V. Kitsios  “Application of a local attractor dimension to reduced space strongly coupled data assimilation for chaotic multiscale systems”  In _Nonlinear Processes in Geophysics_ 27.1, 2020, pp. 51–74  DOI: 10.5194/npg-27-51-2020
  * (66) Courtney Quinn, Dylan Harries and Terence J. O’Kane  “Dynamical Analysis of a Reduced Model for the North Atlantic Oscillation”  In _Journal of the Atmospheric Sciences_ 78.5 Boston MA, USA: American Meteorological Society, 2021, pp. 1647–1671  DOI: 10.1175/JAS-D-20-0282.1
  * (67) Courtney Quinn, Terence J. O’Kane and Dylan Harries  “Systematic calculation of finite-time mixed singular vectors and characterization of error growth for persistent coherent atmospheric disturbances over Eurasia”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 32.2, 2022, pp. 023126  DOI: 10.1063/5.0066150
  * (68) Stephan Rasp  “Coupled online learning as a way to tackle instabilities and biases in neural network parameterizations: general algorithms and Lorenz 96 case study (v1.0)”  In _Geoscientific Model Development_ 13.5, 2020, pp. 2185–2196  DOI: 10.5194/gmd-13-2185-2020
  * (69) Stephan Rasp, Peter D. Dueben, Sebastian Scher, Jonathan A. Weyn, Soukayna Mouatadid and Nils Thuerey  “WeatherBench: A Benchmark Data Set for Data-Driven Weather Forecasting”  In _Journal of Advances in Modeling Earth Systems_ 12.11, 2020  DOI: 10.1029/2020MS002203
  * (70) Stephan Rasp, Michael S. Pritchard and Pierre Gentine  “Deep learning to represent subgrid processes in climate models”  In _Proceedings of the National Academy of Sciences_ 115.39, 2018, pp. 9684–9689  DOI: 10.1073/pnas.1810286115
  * (71) Markus Reichstein, Gustau Camps-Valls, Bjorn Stevens, Martin Jung, Joachim Denzler, Nuno Carvalhais and Prabhat  “Deep learning and process understanding for data-driven Earth system science”  In _Nature_ 566.7743 Springer US, 2019, pp. 195–204  DOI: 10.1038/s41586-019-0912-1
  * (72) O.E. Rössler  “An equation for continuous chaos”  In _Physics Letters A_ 57.5, 1976, pp. 397–398  DOI: 10.1016/0375-9601(76)90101-8
  * (73) David Ruelle  “Ergodic theory of differentiable dynamical systems”  In _Publications Mathématiques de l’IHÉS_ 50 Institut des Hautes Études Scientifiques, 1979, pp. 27–58  URL: http://www.numdam.org/item/PMIHES_1979__50__27_0/
  * (74) Marco Sandri  “Numerical calculation of Lyapunov exponents”  In _Mathematica Journal_ 6.3, 1996, pp. 78–84
  * (75) M.. Schultz, C. Betancourt, B. Gong, F. Kleinert, M. Langguth, L.. Leufen, A. Mozaffari and S. Stadtler  “Can deep learning beat numerical weather prediction?”  In _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 379.2194, 2021, pp. 20200097  DOI: 10.1098/rsta.2020.0097
  * (76) Yakov Sinai  “Kolmogorov-Sinai entropy”  In _Scholarpedia_ 4.3, 2009, pp. 2034  DOI: 10.4249/scholarpedia.2034
  * (77) Jasper Snoek, Hugo Larochelle and Ryan P. Adams  “Practical Bayesian Optimization of Machine Learning Algorithms”  In _Proceedings of the 25th International Conference on Neural Information Processing Systems - Volume 2_ , NIPS’12, 2012, pp. 2951–2959
  * (78) Maike Sonnewald, Redouane Lguensat, Daniel C Jones, Peter D Dueben, Julien Brajard and V Balaji  “Bridging observations, theory and numerical simulation of the ocean using machine learning”  In _Environmental Research Letters_ 16.7 IOP Publishing, 2021, pp. 073008  DOI: 10.1088/1748-9326/ac0eb0
  * (79) Colin Sparrow  “The Lorenz Equations”  Springer New York, NY, 1982  DOI: 10.1007/978-1-4612-5767-7
  * (80) Gilbert Strang  “Introduction to Linear Algebra”  Wellesley-Cambridge Press, 2016
  * (81) Steven H. Strogatz  “Nonlinear Dynamics and Chaos”  CRC Press, 2018  DOI: 10.1201/9780429492563
  * (82) Maxime Tondeur, Alberto Carrassi, Stephane Vannitsem and Marc Bocquet  “On Temporal Scale Separation in Coupled Data Assimilation with the Ensemble Kalman Filter”  In _Journal of Statistical Physics_ 179.5, 2020, pp. 1161–1185  DOI: 10.1007/s10955-020-02525-z
  * (83) Z. Toth and E. Kalnay  “Ensemble Forecasting at NCEP and the Breeding Method”  In _Monthly Weather Review_ 12.125, 1997, pp. 3297–3319  URL: https://doi.org/10.1175/1520-0493(1997)125%3C3297:EFANAT%3E2.0.CO;2
  * (84) F. Uboldi and A. Trevisan  “Multiple-scale error growth in a convection-resolving model”  In _Nonlinear Processes in Geophysics_ 22.1, 2015, pp. 1–13  DOI: 10.5194/npg-22-1-2015
  * (85) Stéphane Vannitsem  “Predictability of large-scale atmospheric motions: Lyapunov exponents and error dynamics”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 27.3, 2017, pp. 032101  DOI: 10.1063/1.4979042
  * (86) Stéphane Vannitsem and Wansuo Duan  “On the use of near-neutral Backward Lyapunov Vectors to get reliable ensemble forecasts in coupled ocean–atmosphere systems”  In _Climate Dynamics_ 55 Springer, 2020, pp. 1125–1139
  * (87) Stéphane Vannitsem and Valerio Lucarini  “Statistical and dynamical properties of covariant lyapunov vectors in a coupled atmosphere-ocean model—multiscale effects, geometric degeneracy, and error dynamics”  In _Journal of Physics A: Mathematical and Theoretical_ 49.22 IOP Publishing, 2016, pp. 224001  DOI: 10.1088/1751-8113/49/22/224001
  * (88) Edoardo Vecchi, Luká$v$ s Pospí$v$ sil, Steffen Albrecht, Terence J. O’Kane and Illia Horenko  “eSPA+: Scalable Entropy-Optimal Machine Learning Classification for Small Data Problems”  In _Neural Computation_ 34.5, 2022, pp. 1220–1255  DOI: 10.1162/neco_a_01490
  * (89) Hao Wang and Dit-Yan Yeung  “Towards Bayesian Deep Learning: A Framework and Some Existing Methods”  In _IEEE Transactions on Knowledge and Data Engineering_ 28.12, 2016, pp. 3395–3408  DOI: 10.1109/TKDE.2016.2606428
  * (90) M.. Wilk and R. Gnanadesikan  “Probability plotting methods for the analysis for the analysis of data”  In _Biometrika_ 55.1, 1968, pp. 1–17  DOI: 10.1093/biomet/55.1.1
  * (91) Alan Wolf, Jack B. Swift, Harry L. Swinney and John A. Vastano  “Determining Lyapunov exponents from a time series”  In _Physica D: Nonlinear Phenomena_ 16.3, 1985, pp. 285–317  DOI: 10.1016/0167-2789(85)90011-9
  * (92) Christopher L. Wolfe and Roger M. Samelson  “An efficient method for recovering Lyapunov vectors from singular vectors”  In _Tellus A_ 59.3, 2007, pp. 355–366  DOI: 10.1111/j.1600-0870.2007.00234.x
  * (93) Yonghui Wu, Mike Schuster, Zhifeng Chen, Quoc V. Le, Mohammad Norouzi, Wolfgang Macherey, Maxim Krikun, Yuan Cao, Qin Gao, Klaus Macherey, Jeff Klingner, Apurva Shah, Melvin Johnson, Xiaobing Liu, Łukasz Kaiser, Stephan Gouws, Yoshikiyo Kato, Taku Kudo, Hideto Kazawa, Keith Stevens, George Kurian, Nishant Patil, Wei Wang, Cliff Young, Jason Smith, Jason Riesa, Alex Rudnick, Oriol Vinyals, Greg Corrado, Macduff Hughes and Jeffrey Dean  “Google’s Neural Machine Translation System: Bridging the Gap between Human and Machine Translation”, 2016  arXiv:1609.08144 [cs.CL]
  * (94) Barret Zoph, Vijay Vasudevan, Jonathon Shlens and Quoc V. Le  “Learning Transferable Architectures for Scalable Image Recognition”  In _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 8697–8710  DOI: 10.1109/CVPR.2018.00907