∎

11institutetext: G. Margazoglou 22institutetext: Imperial College London,
Aeronautics Department, South Kensington Campus London SW7 2AZ, United Kingdom
22email<EMAIL_ADDRESS>33institutetext: L.
Magri44institutetext: Imperial College London, Aeronautics Department, South
Kensington Campus London SW7 2AZ, United Kingdom 55institutetext: The Alan
Turing Institute, 96 Euston Road, NW1 2DB, London, United Kingdom
55email<EMAIL_ADDRESS>
# Stability analysis of chaotic systems from data

Georgios Margazoglou  Luca Magri

(Received: N/A / Accepted: N/A)

###### Abstract

The prediction of the temporal dynamics of chaotic systems is challenging
because infinitesimal perturbations grow exponentially. The analysis of the
dynamics of infinitesimal perturbations is the subject of stability analysis.
In stability analysis, we linearize the equations of the dynamical system
around a reference point, and compute the properties of the tangent space
(i.e., the Jacobian). The main goal of this paper is to propose a method that
learns the Jacobian, thus, the stability properties, from observables (data).
First, we propose the Echo State Network (ESN) with the Recycle Validation as
a tool to accurately infer the chaotic dynamics from data. Second, we
mathematically derive the Jacobian of the Echo State Network, which provides
the evolution of infinitesimal perturbations. Third, we analyse the stability
properties of the Jacobian inferred from the ESN, and compare them with the
benchmark results obtained by linearizing the equations. The ESN correctly
infers the nonlinear solution, and its tangent space with negligible numerical
errors. In detail, (i) the long-term statistics of the chaotic state; (ii) the
covariant Lyapunov vectors; (ii) the Lyapunov spectrum; (iii) the finite-time
Lyapunov exponents; (iv) and the angles between the stable, neutral, and
unstable splittings of the tangent space (the degree of hyperbolicity of the
attractor). This work opens up new opportunities for the computation of
stability properties of nonlinear systems from data, instead of equations.

###### Keywords:

Data-driven learning, Lyapunov Exponents, Covariant Lyapunov Vectors, Echo
State Network

###### MSC:

68T07 34D08 37D45 37M22

††journal: -

##  1 Introduction

Chaotic behavior has been observed and extensively studied in diverse
scientific fields, initially in meteorologyLorenz63  and later in
physicsOtt2002 , chemistry, biology and engineeringStrogatz2018  to name a
few. Chaos appears from deterministic nonlinear equations in the form of
sensitivity to initial conditions, aperiodic behavior, and short
predictability. A successful mathematical tool for the analysis of chaos is
provided by stability analysis. By applying infinitesimal perturbations to a
system’s trajectory, we can classify its stability along different directions,
and compute the properties of its linear tangent space.

Stability analysis relies on the linearization of the dynamical equations,
which requires the Jacobian of the system. The key quantities that
characterize chaotic dynamics, and many other related physical properties,
such as dynamical entropies and fractal dimensions, are the Lyapunov Exponents
(LEs)Ruelle1980 ; Eckmann_Ruelle1985 , which are the eigenvalues of the
Oseledets matrixOseledets1968 . There are several numerical methods to extract
the LEs based on the Gram-Schmidt orthogonalization procedureBenettin1980 ;
Shimada_Nagashima1979 ; Eckmann_Ruelle1985 . The relevant eigenvectors are the
corresponding Lyapunov vectors that constitute a coordinate dependent
orthogonal basis of the linear tangent space. Instead, an intrinsic and norm-
independent basis, which is also time invariant and covariant with the
dynamics is given by the Covariant Lyapunov vectors (CLVs). Crucially, CLVs
are able to provide information on the local structure of chaotic
attractorsGinelli2013 . This viewpoint allows the study of an attractor’s
topology with the occurrence of critical transitionsSchubert2015 ;
Vannitsem2016 ; Sharafi2017 ; Brugnago2020 , paving the way for CLVs to be
considered as precursors to such phenomena.

The previous exposition is traditionally related with model-based approaches,
as it relies on the knowledge of a system’s dynamical equations. However,
studying the stability properties of observed data, where equations are not
necessarily known, is hard; there are few approaches e.g. Wolf1985 ;
Rosenstein1993 , relying on the delayed coordinates attractor reconstruction
by TakensTakens1981 . The recent breakthrough of data-driven (model-free)
approaches poses the reasonable question: Can we use the rich knowledge of
dynamical systems theory for model-free approaches? Indeed, although at early
steps, the use of advanced Machine Learning (ML) techniques for complex
systems has shown promising potential in applications ranging from weather and
climate prediction and classification Rasp2018 ; Dueben2018 ; Margazoglou2021
to fluid flows prediction and optimization Verma2018 ; Kochkov2021 ;
Huhn_Magri_2022 , among others. The overarching goal of this work is to
propose a machine learning approach to accurately learn and infer the ergodic
properties of prototypical chaotic attractors, and in particular to extract
LEs and CLVs from data.

The Recurrent Neural Networks (RNNs) constitute a promising type of ML to
address chaotic behavior. Thanks to their architecture, the RNNs are suitable
for processing sequential data, typically encountered in speech and language
recognition, or time-series predictionGoodfellow2016 . In particular, they are
proven to be universal approximatorsSchafer2006 ; Grigoryeva2018 , and are
able to capture long-term temporal patterns, i.e. they possess memory. A key
piece of their architecture is that they maintain a hidden state that evolves
dynamically, effectively allowing the RNNs to be treated as dynamical systems,
and in particular as discrete neural differential equationsChen2018 . Thus,
RNNs lend themselves to being analysed with dynamical systems theory, allowing
the study of stability properties from the dynamics they have learned. By
exploiting this here, we derive the RNN’s Jacobian and infer the linear
dynamics from data.

Recently there have been significant advancements in employing RNNs to learn
chaotic dynamics Lu2018 ; Pathakchaos2017 ; PathakPRL2018 ; Vlachas2018 ;
Vlachas2020 ; Borra2020 ; Doan2020 ; Doan2021_prsa ; Racca2021 , where two
core objectives are studied: (1) the time-accurate prediction of chaotic
fluctuations and maximization of the prediction horizon; and (2) accurately
learning the ergodic properties of chaotic attractors. The first objective has
been addressed by one of the co-authors in Doan2020 ; Doan2021_prsa ;
Racca2021  for several prototypical chaotic dynamical systems using the same
RNN architecture as the present work. Here we address the second objective by
extending the recent worksPathakchaos2017 ; PathakPRL2018 ; Vlachas2020 ,
where the LEs of the Lorenz 63Lorenz63  and the one-dimensional Kuramoto-
Sivashinsky equation Kuramoto1978  were retrieved from trained RNNs.

In this work we employ a specific architecture of the RNN, a type of reservoir
computer, the Echo State Network (ESN)Lukosevicius2012  and train it with a
diverse set of four prototypical chaotic attractors. The objective of this
paper is twofold; first the accurate learning and inference of the ergodic
properties of the chaotic attractors by the ESN. This is accomplished by
thoroughly comparing the long-term statistics of (i) degrees of freedom, (ii)
LEs, (iii) finite-time LEs, and (iv) angles of the CLVs. Second, by comparing
the distribution of (i) finite-time LEs, and (ii) angles of CLVs on the
topology of the attractor, providing a strict test of the ESN’s capability to
accurately learn intrinsic chaotic properties.

The paper is organized as follows. Section 2 presents the necessary tools for
our study. In particular, Sec. 2.1 provides a brief introduction to the
relevant concepts and quantities from dynamical systems, such as LEs and CLVs.
Then Sec. 2.2 describes the architecture of the ESN, while Sec. 2.3 its
validation strategies. Section 3 presents our main results, which are divided
into two subsections; Sec. 3.1 devoted in low-dimensional systems, namely the
Lorenz 63 Lorenz63  and RösslerRossler1976  attractors; and Sec. 3.2 showing
results on the Charney-DeVoreCdV1979  and the Lorenz 96Lorenz96  attractors.
Finally, we summarize our results and provide future perspectives in the
conclusions in Sec. 4. The appendix A presents the two algorithms to extract
the LEs and CLVs from the ESN.

##  2 Background

In the following two subsections, we summarise the key theory that underpins
the stability of chaotic systems (Sec. 2.1) and reservoir computers (Sec.
2.2).

###  2.1 Stability of chaotic systems

We consider a state ${\bm{x}}(t)\in\mathbb{R}^{D}$ with $D$ degrees of
freedom, which is governed by a set of nonlinear ordinary differential
equations

$\frac{d{\bm{x}}}{dt}=f({\bm{x}}),$ (1)

where $f({\bm{x}}):\mathbb{R}^{D}\to\mathbb{R}^{D}$ is a sufficiently smooth
nonlinear function. Equation (1) defines an autonomous dynamical system.
Hence, the dynamical system exists in a phase space of dimension $D$, equipped
with a suitable metric, and is associated with a certain measure $\mu$ that we
assume to be preserved (invariant). To investigate the stability of the
dynamical system (1), we perturb the state by first-order perturbations as

$\displaystyle{\bm{x}}+{\bm{u}},$
$\displaystyle{\bm{x}}\sim\mathcal{O}(1),\;\;\;{\bm{u}}\sim\mathcal{O}(\epsilon),\;\;\;\epsilon\to
0.$ (2)

By substituting decomposition (2) into (1) and collecting the first-order
terms $\sim\mathcal{O}(\epsilon)$, we obtain the governing equation for the
first-order perturbations (i.e., linear dynamics)

$\frac{d{\bm{u}}}{dt}={\bm{J}}({\bm{x}}(t)){\bm{u}},$ (3)

where $J_{ij}=\frac{\partial f(x_{i})}{\partial x_{j}}\in\mathbb{R}^{D\times
D}$ are the components of the Jacobian, ${\bm{J}}({\bm{x}}(t))$, which is in
general time-dependent. The perturbations ${\bm{u}}$ evolve on the linear
tangent space at each point ${\bm{x}}(t)$. The goal of stability analysis is
to compute the growth rate of infinitesimal perturbations, which is achieved
by computing the Lyapunov exponents and a basis of the tangent space with the
Covariant Lyapunov Vectors. To do so, we numerically time-march $K\leq D$
tangent vectors, ${\bm{u}}_{i}\in\mathbb{R}^{D}$, as columns of the matrix
${\bm{U}}\in\mathbb{R}^{D\times K}$,
${\bm{U}}=[{\bm{u}}_{1},{\bm{u}}_{2},\dots,{\bm{u}}_{K}]$

$\frac{d{\bm{U}}}{dt}={\bm{J}}({\bm{x}}(t)){\bm{U}}.$ (4)

Geometrically, Eq. (4) describes the tangent space around the state
${\bm{x}}(t)$. Starting from ${\bm{x}}(t=t_{0})={\bm{x}}_{0}$ and
${\bm{U}}(t_{0})=\mathbb{I}$, Eqs. (1) and (4) are numerically solved with a
time integrator. As explained in the subsequent paragraphs, in a chaotic
system, almost all nearby trajectories diverge exponentially fast with an
average rate equal to the leading Lyapunov exponent. Hence, the tangent
vectors align exponentially fast with the leading Lyapunov vector,
${\bm{u}}_{1}$. (‘Almost all’ means that the set of perturbations that do not
grow with the largest Lyapunov exponents has a zero measure.) To circumvent
this numerical issue, it is necessary to periodically orthonormalize the
tangent space basis during time evolution, using a QR-decomposition of
${\bm{U}}$, as ${\bm{U}}(t)={\bm{Q}}(t){\bm{R}}(t,\Delta t)$ (see Benettin1980
; Shimada_Nagashima1979 ) and updating the columns of ${\bm{U}}$ with the
columns of ${\bm{Q}}$, i.e. ${\bm{U}}\leftarrow{\bm{Q}}$. The matrix
${\bm{R}}(t,\Delta t)\in\mathbb{R}^{K\times K}$ is upper-triangular and its
diagonal elements $[{\bm{R}}]_{i,i}$ are the local growth rates over a time
span $\Delta t$ of the (now) orthonormal vectors ${\bm{U}}$, which are also
known as backward Gram-Schmidt vectors (GSVs)Ginelli2007 ; Ginelli2013 . The
Lyapunov spectrum, which is sorted from the largest to the smallest exponent,
contains $K\leq D$ exponents and is computed by taking the time average of the
logarithms of $[{\bm{R}}]_{i,i}$,

$\lambda_{i}=\lim\limits_{T\to\infty}\frac{1}{T}\int_{t_{0}}^{T}\ln[{\bm{R}}(t,\Delta
t)]_{i,i}dt.$ (5)

The algorithm 1 in the appendix A is a pseudocode for the calculation of the
LEs for the ESN following Pathakchaos2017 ; Vlachas2020 . The sign of the
Lyapunov exponents indicates the type of the attractor. If the leading
exponent is negative, $\lambda_{1}<0$, the attractor is a fixed point. If
$\lambda_{1}=0$, and the remaining exponents are negative, the attractor is a
periodic orbit. If at least a Lyapunov exponent is positive, $\lambda_{1}>0$,
the attractor is chaotic. In chaotic systems, the Lyapunov time
$\tau_{\lambda}=\frac{1}{\lambda_{1}}$ defines a characteristic timescale for
two nearby orbits to separate, which gives a scale of the system’s
predictability horizonBoffetta2002 . The existence of the Lyapunov exponents,
under fairly general conditions, is proven by the Oseledets’
theoremOseledets1968 ; Eckmann_Ruelle1985 .

The GSVs, ${\bm{U}}$, constitute a norm-dependent orthonormal basis, which is
not time-reversible, due to the frequent orthogonalizations via the QR
decomposition. Instead, the Covariant Lyapunov Vectors (CLVs)
${\bm{V}}=[{\mathbf{v}}_{1},{\mathbf{v}}_{2},\dots,{\mathbf{v}}_{K}]$ (each
CLV ${\mathbf{v}}_{i}\in\mathbb{R}^{D}$ is a column of ${\bm{V}}$) form a
norm-independent and time-invariant basis of the tangent space, which is
covariant with the dynamics. The latter features of the CLVs, which are not
possessed by the GSVs, allow us to examine individual expanding and
contracting directions of a given dynamical system, thus providing an
intrinsic geometrical interpretation of the attractorEckmann_Ruelle1985 ;
Ginelli2013 , as well as a hierarchical decomposition of spatiotemporal chaos,
thanks to their generic localization in physical spaceGinelli2007 . Each
bounded non-zero CLV, i.e. $0<||{\mathbf{v}}_{i}||<\infty$, satisfies the
following equation

$\frac{d{\mathbf{v}}_{i}}{dt}={\bm{J}}({\bm{x}}(t)){\mathbf{v}}_{i}-\lambda_{i}{\mathbf{v}}_{i},$
(6)

which shows that the CLV is evolved by the tangent dynamics
${\bm{J}}({\bm{x}}(t)){\mathbf{v}}_{i}$, while the extra term
$-\lambda_{i}{\mathbf{v}}_{i}$ guarantees that its norm is
boundedHuhn_Magri_2020 . The name “covariant” means that the $i$th CLV at time
$t_{1}$, ${\mathbf{v}}_{i}({\bm{x}}(t_{1}))$, maps at
${\mathbf{v}}_{i}({\bm{x}}(t_{2}))$ at time $t_{2}$, and vice versa.
Mathematically, if ${\bm{M}}(t,\Delta t)=\exp(\int_{t}^{t+\Delta
t}{\bm{J}}({\bm{x}},\tau)d\tau)$ is the system’s tangent evolution operator
(which contains a path-ordered exponential), covariant means
${\bm{M}}(t,\Delta t){\mathbf{v}}_{i}(t)={\mathbf{v}}_{i}(t+\Delta t)$; time-
invariance of CLVs naturally arises from the previous expression, as
${\bm{M}}(t,-\Delta t){\mathbf{v}}_{i}(t+\Delta t)={\mathbf{v}}_{i}(t)$. If
the Lyapunov spectrum is non-degenerate (such as for the cases considered
here) each CLV ${\mathbf{v}}_{i}$ is associated with the Lyapunov exponent
$\lambda_{i}$, and is uniquely defined (up to a phase).

An important subclass of chaotic systems are uniformly hyperbolic systems,
which have a uniform splitting between expanding and contracting directions,
i.e., there are no tangencies between the unstable, neutral, and stable
subspaces Ruelle1997  that form the tangent space. Because of their simple
geometrical structure many theoretical tools have been developed in recent
years. Hyperbolic systems have structurally stable dynamics, and linear
response, meaning that their statistics vary smoothly with parameter
variationsLucarini2014 . In practice, violations of hyperbolicity are commonly
reported in the literatureBlonigan2017 ; Wormell2022 ; Huhn_Magri_2020 ,
whereas true hyperbolic systems are rareKuznetsov2012 . Thanks to the chaotic
hypothesis Ruelle1980 ; Gallavotti_Coen1995_PRL ; Gallavotti_Coen1995_JSP ,
high-dimensional chaotic systems can be practically treated as hyperbolic
systems, i.e. using techniques developed for hyperbolic systems, regardless of
hyperbolicity violations. This is because many convenient statistical
properties of uniformly hyperbolic systems, such as ergodicity, existence of
physical invariant measures, exponential mixing and well-defined time averages
with large deviation laws LebowitzSpohn1999 ; Lepri1998 , can be found in the
macroscopic scale dynamics of certain large non-uniformly hyperbolic
systemsLucarini2014 .

An application of CLVs is to assess the degree of hyperbolicity of the
underlying chaotic dynamics. The tangent space of hyperbolic systems, at each
point ${\bm{x}}$, can be directly decomposed into three invariant subspaces,
${\bm{E}}^{U}_{\bm{x}}\oplus{\bm{E}}^{N}_{\bm{x}}\oplus{\bm{E}}^{S}_{\bm{x}}$.
Here ${\bm{E}}^{U}_{\bm{x}}$ is the unstable subspace composed by the CLVs
associated with positive LEs, ${\bm{E}}^{N}_{\bm{x}}$ is the neutral subspace
spanned by the CLVs associated with the zero LEs, and ${\bm{E}}^{S}_{\bm{x}}$
is the stable subspace spanned by the CLVs associated with negative LEs. In
hyperbolic systems, the distribution of angles between subspaces is bounded
away from zero. In Sec. 3 we will study in detail the angles $\theta_{U,N}$,
$\theta_{U,S}$, and $\theta_{N,S}$ between pairs of the subspaces, and compare
the ability of the ESN to accurately learn both the long-term statistics, and
the phase space finite-time variability of the angles. Because the GSVs are
mutually orthogonal, they cannot assess the degree of hyperbolicity of the
attractor. Moreover, CLVs are key to the optimisation of chaotic acoustic
oscillations Huhn_Magri_2020 , as well as in reduced-order modellingYang2009 ;
they can reveal two uncoupled subspaces of the tangent space, one that
comprises the physical modes carrying the relevant information of the
trajectory, and another composed of strongly decaying spurious
modesGinelli2013 . Two recent attempts to extract CLVs from data driven
approaches, which do not employ a neural network, can be found inViennet2022 ;
Christoph2022 .

We explain the algorithm we employ to compute the CLVs; for further details we
refer the interested reader to Ginelli2007 ; Ginelli2013 ; Huhn_Magri_2020 .
The GSVs are generated by numerically solving Eqs. (1) and (4) simultaneously,
and performing a QR-decomposition every $m$ timesteps. In this way, after a
time-lapse $\Delta t$, the GSVs at time $t+\Delta t$ are given by

${\bm{M}}(t,\Delta t){\bm{U}}(t)={\bm{U}}(t+\Delta t){\bm{R}}(t,\Delta t).$
(7)

We can define the CLVs ${\bm{V}}(t)$ in terms of the GSVs as

${\bm{V}}(t)={\bm{U}}(t){\bm{C}}(t),$ (8)

where ${\bm{C}}$ is an upper triangular matrix that contains the CLV expansion
coefficients, $[{\bm{C}}(t)]_{ji}=c^{j,i}(t)$, for $j\leq i$. Hence, the
objective is to calculate ${\bm{C}}(t)$. Because the CLVs have by choice a
unit norm, each column of the matrix ${\bm{C}}$ has to be normalized
independently, i.e. $\sum_{j=1}^{i}(c_{j,i}(t))^{2}=1,\forall i$.

We start by writing the evolution equation of the CLVs as

${\bm{M}}(t,\Delta t){\bm{V}}(t)={\bm{V}}(t+\Delta t){\mathbf{D}}(t,\Delta
t).$ (9)

We can re-write Eq. (9) via Eq. (8)

$\displaystyle{\bm{U}}(t+\Delta t){\bm{C}}(t+\Delta t){\mathbf{D}}(t,\Delta
t)$ $\displaystyle=$ $\displaystyle{\bm{M}}(t,\Delta t){\bm{U}}(t){\bm{C}}(t)$
$\displaystyle=$ $\displaystyle{\bm{U}}(t+\Delta t){\bm{R}}(t,\Delta
t){\bm{C}}(t),$

and solve with respect to ${\bm{C}}(t)$

$\displaystyle{\bm{C}}(t)={\bm{R}}^{-1}(t,\Delta t){\bm{C}}(t+\Delta
t){\mathbf{D}}(t,\Delta t).$ (11)

This equation is evolved backwards in time starting from the end of the
forward-in-time simulation. We employ the solve_triangular routine of
scipySciPy2020  to invert ${\bm{R}}(t,\Delta t)$ and solve with respect to
${\bm{C}}(t)$. The ${\bm{C}}$ and ${\mathbf{D}}$ matrices are initialized to
the identity matrix $\mathbb{I}$. We leave a sufficient spin-up and spin-down
transient time at the beginning and end of our total time window, before we
compute the CLVs via Eq. (8), to ensure that they are converged. The algorithm
2 in the appendix A is a pseudocode for the calculation of the CLVs.

To estimate the expansions and contractions of the tangent space on finite-
time intervals of length $\Delta t=t_{2}-t_{1}$, we compute the finite time
Lyapunov exponents (FTLEs) as $\Lambda_{i}=\frac{1}{\Delta
t}\ln[{\bm{R}}]_{i,i}$. Hence, $\lambda_{i}$ is the long-time average of
$\Lambda_{i}$. The FTLE $\Lambda_{1}$ physically quantifies the exponential
growth rate of a vector ${\bm{u}}_{1}$ during the time interval $\Delta t$,
therefore, $\Lambda_{2}$ quantifies the exponential growth rate of the vector
${\bm{u}}_{2}$ that is orthogonal to ${\bm{u}}_{1}$ by construction. Hence, as
the GSVs form an orthogonal basis, looking at individual FTLEs for
$\Lambda_{i}$, $i\geq 2$, lacks a physical meaning. Instead, the sum of the
first $n$ FTLEs is a growth rate in $\Delta t$ for a typical $n$-dimensional
volume ${\rm Vol_{n}}$ in the tangent space Shimada_Nagashima1979 ;
kuptsov2018

$\displaystyle\sum_{i=1}^{n}\Lambda_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{\Delta t}\sum_{i=1}^{n}\ln[{\bm{R}}(t,\Delta
t)]_{i,i}=\frac{1}{\Delta t}\ln\prod_{i=1}^{n}[{\bm{R}}(t,\Delta t)]_{i,i}$
(12) $\displaystyle=$ $\displaystyle\frac{1}{\Delta t}\ln{{\rm Vol_{n}}(\Delta
t)}.$

Accordingly, the diagonal matrix ${\mathbf{D}}(t,\Delta t)$ contains the CLV
local growth factors of $\gamma_{i}(t,\Delta t)=||{\bm{M}}(t,\Delta
t){\mathbf{v}}_{i}(t)||$, i.e. $[{\mathbf{D}}(t,\Delta
t)]_{i,j}=\delta_{i,j}\gamma_{i}(t,\Delta t)$. We can extract the finite-time
Covariant Lyapunov exponents (FTCLEs) from the logarithm of these growth
factors for a time interval $\Delta t$

$\Lambda^{c}_{i}=\frac{1}{\Delta t}\ln[{\mathbf{D}}]_{ii}.$ (13)

Each FTCLE quantifies a finite-time exponential expansion or contraction rate
along a covariant direction given by ${\mathbf{v}}_{i}$. Hence each individual
FTCLE has a physical interpretation, in contrast to the FTLEs, as explained
before. On the other hand, now the sums of FTCLEs lack a physical
meaningkuptsov2018 . The long-time average of the FTCLEs is equal to the
Lyapunov exponents,
$\lambda_{i}=\lim\limits_{T\to\infty}\frac{1}{T}\int_{t_{0}}^{T}\Lambda^{c}_{i}(t)dt$.

###  2.2 Echo State Network

The solution of a dynamical system is a time series. From a data point of
view, a time series is a sequentially ordered set of values, in which the
order is provided by time. In a discrete setting, time can be thought of as an
ordering index. For sequential data, and hence time series, recurrent neural
networks (RNNs) are designed to infer the temporal dynamics through their
internal hidden state. However, training RNNs, such as Long-Short term memory
(LSTM)Hochreiter1997  networks and Gated Recurrent Units (GRUs)Cho2014 ,
requires backpropagation through time, which can be a demanding computational
task due to the long-lasting time dependencies of the hidden states Werbos1990
. This issue is overcome by Echo State Networks (ESNs)Jaeger2004 ;
Lukosevicius2012 , a RNN that is a type of reservoir computer, of which the
recurrent weights of the hidden state (commonly named “reservoir”) are
randomly assigned and possess low connectivity. Therefore, only the hidden-to-
output weights are trained leading to a simple quadratic optimization problem,
which does not require backpropagation (see Fig. 1(a) for a graphical
representation). The reservoir acts as a memory of the observed state history.
ESNs have demonstrated accurate inference of chaotic dynamics, such as
inLu2018 ; Pathakchaos2017 ; PathakPRL2018 ; Vlachas2018 ; Vlachas2020 ;
Borra2020 ; Doan2020 ; Doan2021_prsa ; Racca2021 ; Huhn2020proc .

An Echo state network maps the state from time index ${\rm{t}}_{i}$ to index
${\rm{t}}_{i+1}$ as follows (with a slight abuse of notation, the discrete
time is denoted ${\rm{t}}_{i}$). The evolution equations of the reservoir
state and output are governed, respectively, by Lukosevicius2012 ; Racca2021

$\displaystyle{\bm{r}}({\rm{t}}_{i+1})$ $\displaystyle=$
$\displaystyle\tanh\left([\hat{{\mathbf{y}}}_{\mathrm{in}}({\rm{t}}_{i});b_{\mathrm{in}}]^{T}\mathbf{W}_{\mathrm{in}}+{\bm{r}}({\rm{t}}_{i})^{T}\mathbf{W}\right),$
(14) $\displaystyle{\mathbf{y}}_{\mathrm{p}}({\rm{t}}_{i+1})$ $\displaystyle=$
$\displaystyle[{\bm{r}}({\rm{t}}_{i+1});1]^{T}\mathbf{W}_{\mathrm{out}};$ (15)

where at any discrete time t the input vector,
${\mathbf{y}}_{\mathrm{in}}({\rm{t}}_{i})\in\mathbb{R}^{N_{y}}$, is mapped
into the reservoir state ${\bm{r}}\in\mathbb{R}^{N_{r}}$, by the input matrix,
$\mathbf{W}_{\mathrm{in}}\in\mathbb{R}^{(N_{y}+1)\times N_{r}}$, where
$N_{r}\gg N_{y}$. The updated reservoir state ${\bm{r}}({\rm{t}}_{i+1})$ is
calculated at each time iteration as a function of the current input
$\hat{{\mathbf{y}}}_{\mathrm{in}}({\rm{t}}_{i})$ and its previous value
${\bm{r}}({\rm{t}}_{i})$ via Eq. (14), and then is involved in the calculation
of the predicted output,
${\mathbf{y}}_{\mathrm{p}}({\rm{t}}_{i+1})\in\mathbb{R}^{N_{y}}$ via Eq. (15).
Here, $\hat{(\;\;)}$ indicates normalization by the range component-wise,
$(^{T})$ indicates matrix transposition, (;) indicates array concatenation,
$\mathbf{W}\in\mathbb{R}^{N_{r}\times N_{r}}$ is the state matrix,
$b_{\mathrm{in}}$ is the input bias and
$\mathbf{W}_{\mathrm{out}}\in\mathbb{R}^{(N_{r}+1)\times N_{y}}$ is the output
matrix. In our applications, the dimension of the input and output vectors is
equal to the dimension of the physical system of Eq. (1), i.e. $N_{y}\equiv
D$.

The matrices $\mathbf{W}_{\mathrm{in}}$ and $\mathbf{W}$ are (pseudo)randomly
generated and fixed, whilst the weights of the output matrix,
$\mathbf{W}_{\mathrm{out}}$, are the only trainable elements of the network.
The input matrix, $\mathbf{W}_{\mathrm{in}}$, has only one element different
from zero per row, which is sampled from a uniform distribution in
$[-\sigma_{\mathrm{in}},\sigma_{\mathrm{in}}]$, where $\sigma_{\mathrm{in}}$
is the input scaling. The state matrix, W, is an Erdös-Renyi matrix with
average connectivity $d$, in which each neuron (each row of $\mathbf{W}$) has
on average only $d$ connections (i.e. non-zero elements), which are obtained
by sampling from a uniform distribution in $[-1,1]$. The echo state property
enforces the independence of the reservoir state on the initial conditions,
which is satisfied by rescaling W by a multiplication factor, such that the
absolute value of the largest eigenvalue Lukosevicius2012 , i.e., the spectral
radius, is smaller than unity. Following Pathakchaos2017 ; Huhn2020proc ;
Racca2021 ; Racca2022 , we add a bias in the input and output layers to break
the inherent symmetry of the basic ESN architecture. Specifically, the input
bias, $b_{\mathrm{in}}$ is a hyperparameter, selected in order to have the
same order of magnitude as the normalized inputs,
$\mathbf{\hat{y}}_{\mathrm{in}}$. Differently, the output bias is determined
by training the weights of the output matrix, $\mathbf{W}_{\mathrm{out}}$.

Figure 1: (a) Schematic representation of the echo state network. (b) Open-
loop and (c) closed-loop configurations.

In Figs. 1(b)-(c) we present the two types of configurations with which the
ESN can run, i.e. in open-loop or closed-loop, respectively. Running in open-
loop is necessary for the training stage, as the input data is fed at each
step, allowing for the calculation of the reservoir timeseries
${\bm{r}}({\rm{t}}_{i})$, ${\rm{t}}_{i}\in[0,T_{\mathrm{train}}]$, which need
to be stored. There is an initial transient time window, the “washout
interval”, where the output ${\mathbf{y}}_{\mathrm{p}}({\rm{t}}_{i})$ is not
computed. This allows for the reservoir state to satisfy the echo state
property, i.e. making it independent of the arbitrarily chosen initial
condition, ${\bm{r}}(t_{0})={0}$, while also synchronizing it with respect to
the current state of the system.

The training of the output matrix, $\mathbf{W}_{\mathrm{out}}$, is performed
after the washout interval, and involves the minimization of the mean square
error between the outputs and the data over the training set

$\textrm{MSE}=\frac{1}{N_{\mathrm{tr}}N_{y}}\sum_{i=0}^{N_{\mathrm{tr}}}||{\mathbf{y}}_{\mathrm{p}}({\rm{t}}_{i})-{\mathbf{y}}_{\mathrm{in}}({\rm{t}}_{i})||^{2},$
(16)

where $||\cdot||$ is the $L_{2}$ norm, $N_{\mathrm{tr}}+1$ is the total number
of data in the training set, and ${\mathbf{y}}_{\mathrm{in}}$ the input data
on which the ESN is trained. Training the ESN is performed by solving with
respect to $\mathbf{W}_{\mathrm{out}}$ via ridge regression of

$(\mathbf{R}\mathbf{R}^{T}+\beta\mathbb{I})\mathbf{W}_{\mathrm{out}}=\mathbf{R}\mathbf{Y}_{\mathrm{d}}^{T},$
(17)

where $\mathbf{R}\in\mathbb{R}^{(N_{r}+1)\times N_{\mathrm{tr}}}$ and
$\mathbf{Y}_{\mathrm{d}}\in\mathbb{R}^{N_{y}\times N_{\mathrm{tr}}}$ are the
horizontal concatenation of the reservoir states with bias,
$[{\bm{r}}({\rm{t}}_{i});1]$, ${\rm{t}}\in[0,T_{\mathrm{train}}]$, and of the
output data, respectively; $\mathbb{I}$ is the identity matrix and $\beta$ is
the Tikhonov regularization parameter Tikhonov1995 .

On the other hand, in the closed-loop configuration (Fig. 1(c)) the output
${\mathbf{y}}_{\mathrm{p}}$ at time step ${\rm{t}}_{i}$, is used as an input
at time step ${\rm{t}}_{i+1}$, in a recurrent manner, allowing for the
autonomous temporal evolution of the network. The closed-loop configuration is
used for validation (i.e. hyperparameter tuning, see Sec.2.3) and testing, but
not for training. For our purposes, we independently train
$N_{\mathrm{ESN}}\in[5,10]$ networks, of which we take the ensemble average to
increase the statistical accuracy of the prediction and evaluate its
uncertainty. We start with $N_{\mathrm{ESN}}=10$ trained networks, but during
post-processing we may discard any network that shows spurious temporal
evolution.

####  2.2.1 Jacobian of the ESN

In this subsection, we mathematically derive the Jacobian of the Echo State
Network. Equations (14)-(15) are a discrete mapOtt2002 ; Vlachas2020 ,

$\displaystyle{\bm{r}}({\rm{t}}_{i+1})$ $\displaystyle=$ $\displaystyle
f({\mathbf{y}}({\rm{t}}_{i}),{\bm{r}}({\rm{t}}_{i}))=\tanh\left([\hat{{\mathbf{y}}}({\rm{t}}_{i});b_{\mathrm{in}}]^{T}\mathbf{W}_{\mathrm{in}}+{\bm{r}}({\rm{t}}_{i})^{T}\mathbf{W}\right),$
$\displaystyle{\mathbf{y}}({\rm{t}}_{i+1})$ $\displaystyle=$
$\displaystyle[{\bm{r}}({\rm{t}}_{i+1});1]^{T}\textbf{W}_{\mathrm{out}},$

and the continuous-time formulae derived for the Lyapunov exponents and CLVs
in Sec. 2.1 can be adapted for a discrete-time system. The Jacobian of the ESN
reservoir is the total derivative of the hidden state dynamics at a single
timestepPathakchaos2017

$\displaystyle{\mathbf{J}}({\bm{r}}({\rm{t}}_{i+1}))$ $\displaystyle=$
$\displaystyle\frac{d{\bm{r}}({\rm{t}}_{i+1})}{d{\bm{r}}({\rm{t}}_{i})}=\frac{df({\mathbf{y}}({\rm{t}}_{i}),{\bm{r}}({\rm{t}}_{i}))}{d{\bm{r}}({\rm{t}}_{i})}$
(18) $\displaystyle=$ $\displaystyle\frac{\partial
f({\mathbf{y}}({\rm{t}}_{i}),{\bm{r}}({\rm{t}}_{i}))}{\partial{\mathbf{y}}({\rm{t}}_{i})}\frac{\partial{\mathbf{y}}({\rm{t}}_{i})}{\partial{\bm{r}}({\rm{t}}_{i})}+\frac{\partial
f({\mathbf{y}}({\rm{t}}_{i}),{\bm{r}}({\rm{t}}_{i}))}{\partial{\bm{r}}({\rm{t}}_{i})}$
$\displaystyle=$
$\displaystyle(1-\tanh^{2}[\cdot])\mathbf{W}^{T}_{\mathrm{in}}\mathbf{W}^{T}_{\mathrm{out}}+(1-\tanh^{2}[\cdot])\mathbf{W}^{T}$
$\displaystyle=$
$\displaystyle(1-{\bm{r}}({\rm{t}}_{i+1})^{2})\left(\mathbf{W}^{T}_{\mathrm{in}}\mathbf{W}^{T}_{\mathrm{out}}+\mathbf{W}^{T}\right),$

where from Eq. (14)
${\bm{r}}({\rm{t}}_{i+1})^{2}=\tanh^{2}([\hat{{\mathbf{y}}}_{\mathrm{in}}({\rm{t}}_{i});b_{\mathrm{in}}]^{T}\,\mathbf{W}_{\mathrm{in}}+{\bm{r}}({\rm{t}}_{i})^{T}\mathbf{W})$
is the updated squared hidden state at timestep ${\rm{t}}_{i+1}$. The Jacobian
of the ESN is cheap to calculate as the expression
$\left(\mathbf{W}^{T}_{\mathrm{in}}\mathbf{W}^{T}_{\mathrm{out}}+\mathbf{W}^{T}\right)$
is a constant matrix, which is fixed after the training of
$\mathbf{W}_{\mathrm{out}}$. The only time-varying component is the hidden
state. The Jacobian ${\mathbf{J}}\in\mathbb{R}^{N_{r}\times N_{r}}$ is used
for the extraction of the Lyapunov spectrum and the CLVs of a trained ESN. We
time-march $D$ Lyapunov vectors ${\bm{u}}_{i}\in\mathbb{R}^{N_{r}}$, and
periodically perform QR decompositions, where
${\bm{Q}}\in\mathbb{R}^{N_{r}\times D}$, and ${\bm{R}}\in\mathbb{R}^{D\times
D}$. The same CLV algorithm described in Sec. 2.1 is employed to extract $D$
covariant Lyapunov vectors ${\mathbf{v}}_{i}\in\mathbb{R}^{N_{r}}$ from a
trained ESN. The pseudocode is given in algorithm 2.

###  2.3 Validation

The dataset is split into three subsets, which are the training, validation,
and testing subsets in a time-ordered fashion. During training the ESN runs in
open-loop, while during validation and testing, the ESN runs in closed-loop
and the prediction at each step becomes the input for the next step. After
training the ESN, its validation is necessary for the determination of the
hyperparameters. The objective is to compute the hyperparameters that minimize
the logarithm of the MSE (16). The logarithm of the MSE is preferred because
the error varies by orders of magnitude for different hyperparameters, as
explained inRacca2021 . In general, instead of Eq. (16), other types of error
functions can be used for the hyperparamer tuning, such as the maximization of
the prediction horizon Pathakchaos2017 ; Doan2020 ; Racca2021  or the
minimization of the kinetic energy differencesRacca2022 . Here the input
scaling, $\sigma_{\mathrm{in}}$, the spectral radius, $\rho$, and the Tikhonov
parameter, $\beta$, are the ESN hyperparameters that are being tuned
Lukosevicius2012 ; Racca2022 . In order to select the optimal hyperparameters,
$\sigma_{\mathrm{in}}$ and $\rho$, we employ a Bayesian Optimization, which is
a strategy for finding the extrema of objective functions that are expensive
to evaluateHoffman2011 ; Racca2022 . Within the optimal
$[\sigma_{\mathrm{in}},\rho]$ we perform a grid search to select $\beta$
Racca2022 . In particular, $[\sigma_{\mathrm{in}},\rho]$ are searched in the
hyperparameter space $[0.1,5]\times[0.1,1]$ in logarithmic scale, while for
$\beta$ we test $\\{10^{-6},10^{-8},10^{-10},10^{-12}\\}$. The Bayesian
Optimization starts from a grid of $6\times 6$ points in the
$[\sigma_{\mathrm{in}},\rho]$ domain, and then it selects five additional
points through the gp-hedge algorithm Hoffman2011 . We set
$b_{\mathrm{in}}=1$, $d=3$ and add Gaussian noise with zero mean and standard
deviation, $\sigma_{n}=0.0005\sigma_{y}$, where $\sigma_{y}$ is the standard
deviation of the data component-wise, to the training and validation data.
Adding noise to the data improves the performance of ESNs in chaotic dynamics
by alleviating overfitting Vlachas2020 . A summary of the hyperparameters is
shown in Table 1.

Table 1: Echo state networks’ hyperparameters. Multiple values indicate that the parameter is optimized within the range. Parameter | Name | Value
---|---|---
$\rho$ | Spectral radius | $[0.1,\,1]$
$\sigma_{\mathrm{in}}$ | Input scaling | $[0.1,\,5]$
$\beta$ | Tikhonov parameter | $\\{10^{-6},10^{-8},10^{-10},10^{-12}\\}$
$d$ | Connectivity | $3$
$b_{\mathrm{in}}$ | Input bias | $1$
$\sigma_{n}$ | Noise (training) | $0.0005\sigma_{u}$

One of the most commonly used validation strategy for RNNs is the Single Shot
Validation (SSV) Lukosevicius2019 , in which the data are split into a
training set, followed by a single small validation set; see Fig. 2(a). As the
ESN now runs in closed-loop, the size of the validation set is limited by the
chaotic nature of the signal. In particular, at the beginning of the
validation set, the input ${\mathbf{y}}({\rm{t}}_{0})$ of the ESN is
initialized to the target value. However, chaos causes the predicted
trajectory to quickly diverge from the target trajectory in a few Lyapunov
times $\tau_{\lambda}$. The validation interval is therefore small and not
representative of the full training set, which causes poor performance in the
test setRacca2022 . An improvement to the performance with cheap computations
is achieved by the the Recycle Validation (RV), which was recently proposed
byRacca2022 . In the RV, the network is trained only once on the entire
training dataset (in open-loop), and validation is performed on multiple
intervals already used for training (but now in closed-loop); see Fig. 2(b).
In this work we use the chaotic Recycle Validation (RVC), where the validation
interval simply shifts as a small multiple of the first Lyapunov exponent,
$N_{\mathrm{val}}=3\lambda_{1}$.

Figure 2: Schematic representation of the (a) Single Shot, and (b) Recycling
Validation strategies. Here, ${\mathbf{y}}$ represents the degrees of freedom
of the data. Three sequential validation intervals are shown for the Recycle
Validation Racca2022 .

##  3 Results

In this section, we present the numerical results, which include a thorough
comparison between the statistics produced by the autonomous temporal
evolution of the ESN and the target dynamical system. The selected observables
are the statistics of the degrees of freedom, the Lyapunov exponents, the
angles between the CLVs or subspaces composed of CLVs, and the finite-time
covariant Lyapunov exponents. We separate our analysis into two subsections,
which contain two low-dimensional systems and then two higher dimensional
systems.

###  3.1 Low dimensional chaotic systems

As a first case, we consider two low-dimensional dynamical systems that
exhibit chaotic behavior: Lorenz 63 (L63)Lorenz63  and Rössler Rossler1976
attractors. The Lorenz 63 system is a reduced-order model of atmospheric
convection for a single thin layer of fluid that is heated uniformly from
below and cooled from above, which is defined by

$\displaystyle\frac{dx_{1}}{dt}$ $\displaystyle=$
$\displaystyle\sigma(x_{2}-x_{1})$ $\displaystyle\frac{dx_{2}}{dt}$
$\displaystyle=$ $\displaystyle x_{1}(\rho-x_{3})-x_{2}$ (19)
$\displaystyle\frac{dx_{3}}{dt}$ $\displaystyle=$ $\displaystyle
x_{1}x_{2}-\beta x_{3}.$

We chose the parameters $[\sigma,\beta,\rho]=[10,8/3,28]$ to ensure a chaotic
behavior. The Rössler attractor, which models equilibrium in chemical
reactions, is governed by

$\displaystyle\frac{dx_{1}}{dt}$ $\displaystyle=$
$\displaystyle-(x_{2}+x_{3})$ $\displaystyle\frac{dx_{2}}{dt}$
$\displaystyle=$ $\displaystyle x_{1}+ax_{2}$ (20)
$\displaystyle\frac{dx_{3}}{dt}$ $\displaystyle=$ $\displaystyle
b+x_{3}(x_{1}-c),$

We choose the parameters $[a,b,c]=[0.1,0.1,14]$ to ensure a chaotic behavior.

To generate the target set we evolve the dynamical systems forward in time
with a 4${}^{\text{th}}$ order Runge-Kutta (RK4) integrator and a timestep
$dt=0.005$ for both L63 and Rössler, which is sufficiently small for a good
temporal resolution. (We tested slightly larger/smaller timesteps with no
significant differences. Results not shown.) We perform a QR decomposition
every $m=1$ timesteps for L63 and every $m=5$ timesteps for Rössler. For all
systems, we generate a training set of size $1000\tau_{\lambda}$, and a test
set of size $4000\tau_{\lambda}$, for the CLV statistics to converge, where
$\tau_{\lambda}=1/\lambda_{1}$ is the Lyapunov time, which is the inverse of
the maximal Lyapunov exponent $\lambda_{1}$.

Figure 3: Comparison of the Target (straight black line) and ESN (red dashed
line) Probability Density Functions (PDF) of the three degrees of freedom,
$x_{1}$, $x_{2}$, and $x_{3}$ of the Lorenz 63 system (3.1) (a-c) and the
Rössler system (3.1) (d-f).

First, we test whether the ESN correctly learns the chaotic attractor from a
statistical point of view, i.e., whether the ESN correctly learns the long-
term statistics of the degrees of freedom when it evolves in the closed-loop
(autonomous) mode. By estimating the probability density function (PDF) of the
degrees of freedom of the ESNs, as a normalized histogram, and comparing it
with the corresponding PDF of the target set, we extract information on the
invariant measure of the considered chaotic system. This is shown in Fig. 3
for L63 and Rössler attractors, in which the black lines show the target
statistics and the red dashed lines show the ESN statistics. In Figs. 3, 5, 7,
and Table 2 we have used $N_{\mathrm{ESN}}$ ESNs trained on $N_{\mathrm{ESN}}$
independent target systems, starting from different initial conditions, and
averaged among the estimated observables, where $N_{\mathrm{ESN}}=6$, and
$N_{\mathrm{ESN}}=8$ for Rössler and L63, respectively. We perform the
ensemble calculation to quantify the uncertainty of the predictions and the
robustness of the ESN for different initializations.

Second, we test whether the ESN correctly learns the Lyapunov spectrum. Table
2 shows the ESN predictions on the Lyapunov exponents for the L63 and Rössler
attractors, which are compared with the target exponents. The leading exponent
is accurately predicted with a 0.2% error in the L63, and 1.5% error in the
Rössler system. In chaotic systems, there exists a neutral Lyapunov exponent,
which is associated with the direction of $\frac{d{\bm{x}}}{dt}$. In these
cases the neutral Lyapunov exponents are $\lambda_{2}=0$ for both systems,
which are correctly inferred by the ESN within a $\mathcal{O}(10^{-5})$ error,
or less. For the smallest, and negative exponent, which is generally harder to
extract because it is highly damped, the relative error is about 0.6% for L63
and 2.1% for Rössler. Therefore, the ESNs can accurately capture the tangent
dynamics of a low dimensional chaotic attractor.

Table 2: Estimates of Lyapunov exponents $\lambda_{i}$ for the two low dimensional systems, the Lorenz 63 and Rössler attractors. Comparison between the target and echo state network. | Lorenz 63 | Rössler
---|---|---
$\lambda_{i}$ | target | ESN | target | ESN
1 | 0.9050 | 0.9067 | 0.071 | 0.070
2 | 9$\times 10^{-5}$ | -8.$\times 10^{-5}$ | $2\times 10^{-6}$ | $1\times 10^{-6}$
3 | -14.572 | -14.664 | -13.88 | -14.17

Third, we investigate the angles between the CLVs. We assess whether the ESNs
learn the long-term statistics of these quantities, but also whether, they
correctly infer the distribution and fluctuations of those observables in the
phase space. In other words, whether the ESNs learn the geometrical structure
of the attractor and its tangent space.

Figure 4: Comparison of Target (left column) and ESN (right column) for a
$300\tau_{\lambda}$ trajectory of the Lorenz 63 system (3.1) in the test set,
coloured by the CLV principal angles (in $\deg$). First row: $\theta_{U,N}$,
second row: $\theta_{U,S}$, and third row: $\theta_{N,S}$.

Figure 5: Comparison of the Target (straight black line) and ESN (red dashed
line) Probability Density Functions (PDF) of the three principal angles
between the Covariant Lyapunov Vectors, where $U$ refers to unstable, $N$ to
Neutral and $S$ to stable CLVs. The top row (a-c) is for Lorenz 63 (3.1) and
the bottom row (d-f) for Rössler (3.1). All $y$-axes are in logarithmic scale
and the $x$-axis is in degrees. The shaded region indicates the error bars
derived by the standard deviation.

In Fig. 4, we present an analysis of the distribution of principal angles
between the CLVs,

$\theta_{a,b}=\frac{180^{\circ}}{\pi}\cos^{-1}(|{\mathbf{v}}_{a}\cdot{\mathbf{v}}_{b}|),$
(21)

$\theta_{a,b}\in[0^{\circ},90^{\circ}]$, on the topology of the L63 attractor.
The attractor is well reproduced by a selected ESN (right column), compared to
the target (left column). The size of both trajectories is
$300\tau_{\lambda}$. In this case, there are three principal angles between
the CLVs; $\theta_{U,N}$ is the angle between the unstable and neutral CLV;
$\theta_{U,S}$ is the angle between the unstable and stable CLV;
$\theta_{N,S}$ is the angle between the neutral and stable CLV. The colouring
of the attractor is associated with the measured $\theta_{a,b}$. The black and
dark red colours identify small angles, i.e., regions of the attractor where
near-tangencies between the CLVs occur. Possible tangencies between CLVs or
invariant manifolds composed of CLVs (as will be discussed later for higher
dimensional chaotic systems) are of significant importance, as they signify
that the attractor is non-hyperbolicGinelli2013  (see Sec. 2.1). Figure 4
shows that the ESN is able to accurately learn the dynamics of the tangent
linear space of the attractor.

In Fig. 5, we show the PDF of the principal angles between the three CLVs, for
which there is agreement between target and ESN results in all cases for both
L63 and Rössler, even for smaller angles. The non-zero count of events close
to $\theta\to 0$ indicates that the two considered systems are non-hyperbolic,
which is consistent with the literaturekuptsov2018 .

Fourth, we analyse the distribution on the attractor, as well as the
statistics, of the Finite-time Covariant Lyapunov Exponents, for a time-lapse
of $\Delta t=m\,dt$ timestep , and assess the accuracy of the trained ESNs.
For the considered low-dimensional systems there are three FTCLEs with each
showing the finite-time growth rate of the corresponding Covariant Lyapunov
Vectors.

In Fig. 6 we visualize the distribution of the single timestep FTCLEs, in the
case of the Rössler attractor, which is well reproduced by a selected ESN
(right column), compared to the target (left column). The size of both
trajectories is $300\tau_{\lambda}$. FTCLE 1 is the finite-time exponent for
the unstable CLV, FTCLE 2 is for the neutral CLV, and FTCLE 3 is for the
stable CLV. The colouring is associated with the values of the FTCLEs. Large
positive FTCLEs correspond to high finite-time growth rates and, thus, reduced
predictability. The distribution of the leading FTCLE on the attractor is
similar between the target and ESN. The second FTCLE, which corresponds to the
neutral CLV also has agreement between the two, while the third FTCLE has some
minor, and hard to detect, disagreement in the region of the attractor where
the large variation in the $z$-axis occur.

Finally, Fig. 7 shows the PDF of the three FTCLEs. There is agreement between
the ESN-inferred quantities and the target in all cases, in particular in the
Rössler attractor for the most-probable statistics. The small deviation in
Fig. 7(a) for L63 corresponds to the statistics around the peak of the first
FTCLE, $\Lambda_{1}^{c}$, but the tails of the distributions are well
reproduced. The mean of the $\Lambda_{i}^{c}$ distributions coincides with the
LEs $\lambda_{i}$, which holds true for all our results. A behavior as in Fig.
7(a) implies that in this case the finite-time values $\Lambda_{1}^{c}$ are
less peaked around the mean value, even though their long-time average
coincides with the Lyapunov exponent $\lambda_{1}$. Nevertheless, in Figs.
7(d-f) for Rössler the statistics around the peak (and beyond) are well
captured.

Figure 6: Comparison of Target (left column) and ESN (right column) for a
$300\tau_{\lambda}$ trajectory of the Rössler system (3.1) in the test set,
coloured by the three FTCLEs. First row: FTCLE 1, second row: FTCLE 2, and
third row: FTCLE 3.

Figure 7: Comparison of the Target (straight black line) and ESN (red dashed
line) Probability Density Functions (PDF) of the three finite-time Covariant
Lyapunov Exponents. The top row (a-c) is for Lorenz 63 (3.1) and the bottom
row (d-f) for Rössler (3.1). All $y$-axes are in logarithmic scale.

We refer the interested reader to our supplementary material where the
corresponding results of Figs. 4,6 for both attractors are shown. Also, the
statistics of FTLEs, as well as their distribution on the chaotic attractors
are presented in the supplementary material.

Figure 8: Comparison of the Target (straight black line) and ESN (red dashed
line) PDF of the first six degrees of freedom, $x_{1}$, …, $x_{6}$ of the (a)
Charney-DeVore (3.2) and (b) for Lorenz 96 (24) for $D=20$.

###  3.2 Higher dimensional chaotic systems

We follow the same analysis and approach as in Sec. 3.1 for two higher
dimensional chaotic systems, both of which are related to atmospheric physics
and meteorology. The first is a reduced-order model of atmospheric blocking
events by Charney and DeVoreCdV1979  (CdV), which is a six-dimensional
truncation of the equations for barotropic flow with orography. We employ the
formulation of DeSwart1989 ; Crommelin2004 , which is forced by a zonal flow
profile that can be barotropically unstable. The governing equations are

$\displaystyle\frac{dx_{1}}{dt}$
$\displaystyle=\gamma_{1}^{*}x_{3}-C(x_{1}-x_{1}^{*}),$
$\displaystyle\frac{dx_{2}}{dt}$
$\displaystyle=-(\alpha_{1}x_{1}-\beta_{1})x_{3}-Cx_{2}-\delta_{1}x_{4}x_{6},$
$\displaystyle\frac{dx_{3}}{dt}$
$\displaystyle=(\alpha_{1}x_{1}-\beta_{1})x_{2}-\gamma_{1}x_{1}-Cx_{3}+\delta_{1}x_{4}x_{5},$
(22) $\displaystyle\frac{dx_{4}}{dt}$
$\displaystyle=\gamma_{2}^{*}x_{6}-C(c_{4}-x_{4}^{*})+\epsilon(x_{2}x_{6}-x_{3}x_{5}),$
$\displaystyle\frac{dx_{5}}{dt}$
$\displaystyle=-(\alpha_{2}x_{1}-\beta_{2})x_{6}-Cx_{5}-\delta_{2}x_{4}x_{3},$
$\displaystyle\frac{dx_{6}}{dt}$
$\displaystyle=(\alpha_{2}x_{1}-\beta_{2})x_{2}-\gamma_{2}x_{4}-Cx_{6}+\delta_{2}x_{4}x_{2},$

where the model coefficients are

$\displaystyle\alpha_{m}$
$\displaystyle=\frac{8\sqrt{2}}{\pi}\frac{m^{2}}{4m^{2}-1}\frac{b^{2}+m^{2}-1}{b^{2}+m^{2}},\quad\quad\beta_{m}=\frac{\beta
b^{2}}{b^{2}+m^{2}},$ $\displaystyle\delta_{m}$
$\displaystyle=\frac{64\sqrt{2}}{15\pi}\frac{b^{2}-m^{2}+1}{b^{2}+m^{2}},\quad\quad\quad\gamma^{*}_{m}=\gamma\frac{4m}{4m^{2}-1}\frac{\sqrt{2}b}{\pi},$
(23) $\displaystyle\epsilon$
$\displaystyle=\frac{16\sqrt{2}}{5\pi},\quad\quad\quad\quad\quad\quad\gamma_{m}=\gamma\frac{4m^{3}}{4m^{2}-1}\frac{\sqrt{2}b}{\pi(b^{2}+m^{2})}.$

Eq. (3.2) is integrated with RK4 and $dt=0.1$. The constants are set to
$(x_{1}^{*},x_{4}^{*},C,\beta,\gamma,b)=(0.95,-0.76095,0.1,1.25,0.2,0.5)$, for
which the CdV model generates regime transitionsDeSwart1989 ; Crommelin2004 .
In particular, the CdV model allows for two metastable states, the so-called
“zonal” state, which represent the approximately zonally symmetric jet stream
in the mid-latitude atmosphere, and the “blocked” state, which refers to a
diverse class of weather patterns that are a persistent deviation from the
zonal state. Blocking events are known to be associated with regional extreme
weather, from heatwaves in summer to cold spells in winterWoollings2018 . The
dynamical properties of CLVs in connection to blocking events were recently
investigated for a series of more complex atmospheric models than CdV
Vannitsem2016 ; Schubert2015 ; Schubert2016 , which demonstrated that CLVs are
good candidates for blockings precursors, as wells as a good basis for model
reduction. In previous workDoan2020 , the CdV system was used as a training
model for the ESN, with the purpose of studying short-term accurate prediction
of chaos, and quantifying the benefit of Physics informed Echo State
NetworksDoan2020 .

The second higher-dimensional system that we consider is the Lorenz 96 (L96)
modelLorenz96 , which is a system of coupled ordinary differential equations
that describes the large-scale behavior of the mid-latitude atmosphere, and
the transfer of a scalar atmospheric quantity. Three characteristic processes
of atmospheric systems (advection, dissipation, and external forcing) are
included in the model, whose equations are

$\frac{d{\bm{x}}_{i}}{dt}=({\bm{x}}_{i+1}-{\bm{x}}_{i-2}){\bm{x}}_{i-1}-{\bm{x}}_{i}+F$
(24)

where ${\bm{x}}=[x_{1},x_{2},\dots,x_{D}]\in\mathbb{R}^{D}$. We set periodic
boundary conditions, i.e. $x_{1}=x_{D+1}$. In our analysis we chose $D=20$
degrees of freedom. The external forcing is set to $F=8$, which ensures a
chaotic evolutionVlachas2020 . We integrate the system with RK4 and $dt=0.01$.
We perform a QR decomposition every $m=5$ timesteps for CdV and every $m=10$
timesteps for L96. Similar to the previous section, we generate a training set
of size $1000\tau_{\lambda}$ and a test set of size $4000\tau_{\lambda}$.

First, Fig. 8 shows the PDF of the six degrees of freedom of CdV, and the
first six from L96 (the PDFs of the rest 13 dofs have similar shape and
agreement between ESN and target). We use a semilogarithmic scale to emphasize
that the agreement between target (black line) and ESN (red dashed line) is
accurate for the tails of the distributions, which effectively correspond to
the edges of each attractor. As in Sec. 3.1 in order to evaluate uncertainty
and robustness, we start with $N_{\mathrm{ESN}}=10$ trained networks, but
during post-processing we discard any network that shows spurious temporal
evolution, and perform a further averaging of the PDFs of each network’s
observable. Therefore, the PDFs of Fig. 8 are the outcome of averaging
$N_{\mathrm{ESN}}=5$ and $N_{\mathrm{ESN}}=9$ PDFs with the same binning, for
CdV and L96, respectively.

Second, Fig. 9 shows the Lyapunov Exponents spectrum of (a) CdV and (b) L96
for $D=20$, and compares the target (black squares) with the ESN prediction
(red circles). The CdV model has a single positive Lyapunov exponent, with the
average value of 5 ESNs resulting in $\lambda_{1}=0.0214$, and for the 5
independent target sets, $\lambda^{\text{targ}}_{1}=0.0232$ with an $8\%$
absolute error. The second Lyapunov exponent is zero (to numerical error), and
corresponds to the neutral direction, with $\lambda_{2}=-3\times 10^{-5}$ for
ESN, and $\lambda^{\text{targ}}_{2}=-7\times 10^{-6}$ for the target. The low
order of magnitude achieved by the ESN assures its ability to capture the
neutral exponent. Finally, the four remaining negative exponents are well
learned by the ESN, i.e., $\lambda_{3-6}=[-0.077$, $-0.103$, $-0.224$,
$-0.234]$ and $\lambda^{\text{targ}}_{3-6}=[-0.079$, $-0.101$, $-0.218$,
$-0.226]$ for the target. Overall, excluding $\lambda_{2}$, the mean absolute
error of the CdV Lyapunov spectrum here is 3.7%, which is negligibly small.

Figure 9: Comparison of the Target (black squares) and ESN (red circles)
Lyapunov spectrum for (a) Charney-DeVore (3.2), and (b) Lorenz 96 (24) at
$D=20$.

With respect to the L96 Lyapunov spectra in Fig. 9(b), the agreement between
target and ESN across all 20 exponents. In particular, there are 6 positive, 1
zero and 13 negative exponents. The maximal exponent predicted from the
ensemble of $N_{\mathrm{ESN}}=9$ ESNs is equal to $\lambda_{1}=1.551$, and for
the 9 independent target sets, $\lambda^{\text{targ}}_{1}=1.557$, meaning a
$0.4\%$ absolute error. The rest of the positive exponents are well captured
by the ESN, with $\lambda_{2-6}=[1.221$, $0.936$, $0.668$, $0.416$, $0.151]$
and $\lambda^{\text{targ}}_{2-6}=[1.217$, $0.937$, $0.673$, $0.413$, $0.152]$
for the target. The zero exponent is sufficiently small with
$\lambda_{7}=-10^{-4}$ for ESN, and $\lambda^{\text{targ}}_{7}=4\times
10^{-4}$ for target. Albeit more difficult to predict because of large
numerical dissipation, the negative Lyapunov exponents are accurately learned
by the ESN, with the smallest ones reading $\lambda_{15-20}=[-1.84$, $-2.22$,
$-2.71$, $-3.45$, $-4.24$, $-4.73]$ and accordingly
$\lambda^{\text{targ}}_{15-20}=[-1.85$, $-2.21$, $-2.71$, $-3.45$, $-4.25$,
$-4.75]$ for the target. Those directions in tangent space decay exponentially
fast and the accuracy that the ESN achieves is consistent. For L96 the mean
absolute error of the Lyapunov spectrum is approximately 0.5%.

Figure 10: Comparison of the Target (straight black line) and ESN (red dashed
line) PDF of the three minimum principal angles between the three subspaces
composed by the CLVs, where $U$ refers to unstable, $N$ to neutral and $S$ to
stable CLVs. The top row (a-c) is for Charney-DeVore (3.2) and the bottom row
(d-f) for Lorenz 96 (24) at $D=20$. Both $x$ and $y$ axes are in logarithmic
scale and the $x$-axis is in degrees. Only in (a) a logarithmic binning was
used being denser close to $\theta_{U,N}\to 0$, while PDFs in (b-f) are
linearly binned in $x$-axis.

Figure 11: Comparison of the Target (straight black line) and ESN (red dashed
line) PDF of six Finite time Covariant Lyapunov Exponents ($\Lambda^{c}_{i}$)
for (a) Charney-DeVore (3.2) for which $\lambda_{1}>0$, $\lambda_{2}=0$ and
the rest are $\lambda_{i}<0$, and (b) Lorenz 96 (24) at $D=20$, where for
$i=1,2,3$ $\lambda_{i}>0$, for $i=7$ $\lambda_{i}=0$, and for $i=11,17$
$\lambda_{i}<0$. All $y$-axes are in logarithmic scale.

To further elaborate, the L96 is known to be an extensive system
Ruelle1982large ; Karimi2010 , which means that quantities such as the surface
width, the entropy and the attractor dimension scale linearly with its
dimensionality $D$. For the Lyapunov spectrum this means that the proportion
of positive to negative exponents is roughly the same ($\approx 1/2$) as $D$
changes. For this reason, our chosen $D=20$ is sufficient for our purposes.

Third, we investigate the statistics of the principal angles,
$\theta\in[0^{\circ},90^{\circ}]$, between the three subspaces that partition
the invariant manifolds, which are the unstable ${\bm{E}}^{U}_{\bm{x}}$,
neutral ${\bm{E}}^{N}_{\bm{x}}$ and stable ${\bm{E}}^{S}_{\bm{x}}$, spanned by
the corresponding CLVs. The extraction of the principal angles between two
linear subspaces requires a singular value decomposition of their matrix
product ${\bm{\Gamma}}^{a,b}={\bm{E}}^{a}_{\bm{x}}{\bm{E}}^{b}_{\bm{x}}$
(assuming the CLVs are ordered as stacked columns, according to their Lyapunov
exponent order), because all paired products between the CLVs spanning the
subspaces do not provide all the anglesKuptsov2009 ; Ginelli2013 . The angles
are given by

$\theta_{a,b}=\frac{180^{\circ}}{\pi}\cos^{-1}\left({\text{svd}[{\bm{\Gamma}}^{a,b}]}\right),$
(25)

and we analyse the smallest singular value. Here, we use the implemented
routine scipy.linalg.subspace_angles of the scipy packageSciPy2020  in python,
and analyse the minimum angle in order to track homoclinic tangencies between
the subspaces. This implementation is based on the algorithm presented in
Knyazev2002 , which is has improved accuracy with respect to Eq. (25) in the
estimation of small angles.

In Fig. 10, we study the PDFs of the three principal angles between the linear
subspaces for CdV and L96. In CdV, the unstable and neutral subspaces are
spanned only by the corresponding CLVs, while the stable subspace is spanned
by the remaining four CLVs, of which $\lambda_{i}<0$. In L96 with $D=20$, the
unstable subspace is spanned by the first six CLVs, the neutral subspace is
spanned only by the $7^{\text{th}}$ CLV, and the stable subspace is spanned by
the remaining 13 CLVs. Focusing on Figs. 10(a-c) for CdV, we notice that this
system is non-hyperbolic because the PDFs are populated close to $\theta\to
0$. Specifically for Fig. 10(a), the binning is geometrically spaced and
denser close to $\theta\to 0$. Interestingly, the PDF of $\theta_{U,N}$ of CdV
for small angles follows a power-law $PDF(\theta)\sim\theta^{-\alpha}$ for
$\theta\to 0$ and until $\approx 10^{\circ}$, before it saturates. A different
shape that is still highly non-hyperbolic is shown for the PDFs of
$\theta_{U,S}$ and $\theta_{N,S}$ in Figs. 10(b-c), in which the binning is
linear and both axes in logarithmic scale. Figures 10(d-f) show the same
statistics in the case of L96, which is also non-hyperbolic, as there is
strong frequency of tangencies, $\theta\to 0$. In all plots of Fig. 10 the
agreement of the subspace angle statistics between target and ESN is good,
which demonstrates that the ESN has achieved a robust and accurate learning of
the ergodic properties from higher-dimensional data.

Fourth, the statistics of FTCLEs ($\Lambda^{c}_{i}$), for a time-lapse of
$\Delta t=m\,dt$ timesteps, in the cases of CdV and L96 are shown in Fig. 11.
All six $\Lambda^{c}_{i}$ are shown for CdV, while a representative set of six
$\Lambda^{c}_{i}$ are shown for L96, such that $\lambda_{i}>0$ for $k=1,2,3$,
$\lambda_{i}=0$ for $k=7$, and $\lambda_{i}<0$ for $k=11,17$. For CdV, the
most probable statistics are well captured by the ESN, which is in agreement
with the target data. There are slight deviations at the tails of the
distributions, which are still in agreement within error bars (shaded region).
In the case of L96, the agreement is good for both the most probable
statistics and the tails, for all FTCLEs (also those not shown). The first
moment of the distributions, i.e. the mean of the FTCLEs timeseries, must be
equal to the Lyapunov exponents,
$\lambda_{i}=\frac{1}{T}\int_{0}^{T}\Lambda^{c}_{i}$, which indeed holds for
all the cases considered here. The agreement between ESN and target sets in
Fig. 11 shows that the ESN is able to accurately learn the finite-time
variability of the CLV growth rates also for higher dimensional systems that
are characterized by many Lyapunov exponents.

Finally, in Table 3 we show the estimated Kaplan-Yorke dimension
KaplanYorke1979  for all the considered systems and compare the outcomes of
the ESN and target. This dimension is an upper bound of the attractor’s
fractal dimensionOtt2002 , which is defined as

$D_{\text{KY}}=k+\frac{\sum_{i=1}^{k}\lambda_{i}}{|\lambda_{i+1}|},$ (26)

where $k$ is such that the sum of the first $k$ LEs is positive and the sum of
the first $k+1$ LEs is negative. We observe a good agreement in all cases with
$\leq 1\%$ error. This observation further confirms the ability of the ESN to
accurately learn the properties of the chaotic attractor.

Table 3: Estimates of the Kaplan-Yorke dimension for all attractors, comparing between the target and echo state network. The error is the quantity $\frac{\rm target-ESN}{\rm target}\times 100\%$. | Target | ESN | % error
---|---|---|---
Lorenz 63 | 2.0621 | 2.0618 | 0.014
Rossler | 2.0051 | 2.0049 | 0.009
CdV | 2.294 | 2.277 | 0.76
Lorenz 96 | 13.4697 | 13.4721 | 0.017

##  4 Conclusions

Stability analysis is a principled mathematical tool to quantitatively answer
key questions on the behaviour of nonlinear systems: Will infinitesimal
perturbations grow in time (i.e., is the system linearly unstable)? If so,
what are the perturbations’ growth rates (i.e., how linearly unstable is the
system)? What are the directions of growth? To answer these questions,
traditionally, we linearize the equations of the dynamical system around a
reference point, and compute the properties of the tangent space, the dynamics
of which is governed by the Jacobian. The overarching goal of this paper is to
propose a method that infers the stability properties directly from data,
which does not rely on the knowledge of the dynamical differential equations.
We tackle chaotic systems, which have a linearized behaviour that is more
general and intricate than periodic or quasi-periodic oscillations. First, we
propose the Echo State Network with the Recycle Validation as a tool to
accurately learn the chaotic dynamics from data. The data is provided by the
integration of low- and higher- dimensional prototypical chaotic dynamical
systems. These systems are qualitatively different from each other, and are
toy models that describe diverse physical settings, ranging from climatology
and meteorology to chemistry. Second, we mathematically derive the Jacobian of
the Echo State Network (Eq. (18)). In contrast to other recurrent neural
networks, such as long short-term memory networks or gated recurrent units,
the Jacobian of the ESN is mathematically simple and computationally
straightforward. Third, we analyse the stability properties inferred from the
ESN, and compare them with the target properties (ground truth) obtained by
linearizing the equations. The ESN correctly infers quantities that
characterize the chaotic dynamics and its tangent space (i) the long-term
statistics of the solution, for which we compute the probability density
function of each state variable; (ii) the covariant Lyapunov vectors, which
are a physical basis for the tangent space that is covariant with the
dynamics; (ii) the Lyapunov spectrum, which is the set of eigenvalues of the
Oseledets matrix that are the perturbations’ average exponential growths;
(iii) the finite-time Lyapunov exponents, which are the finite-time growth
along the covariant Lyapunov vectors; and (iv) the angles between the stable,
neutral, and unstable splittings of the tangent space, which informs about the
degree of hyperbolicity of the attractor. We show that these quantities can be
accurately learned from data by the ESN, with negligible numerical errors.

As mathematically and numerically shown in Huhn_Magri_2020 , the stability
properties of fixed points (with eigenvalue analysis) and periodic solutions
(with Floquet analysis) can be inferred from covariant Lyapunov analysis.
Therefore, this work opens up new opportunities for the inference of stability
properties from data in nonlinear systems, from simple fixed points, through
periodic oscillations, to chaos.

###### Acknowledgements.

This research has received financial support from the ERC Starting Grant No.
PhyCo 949388. LM gratefully acknowledges financial support TUM Institute for
Advanced Study (German Excellence Initiative and the EU 7th Framework
Programme No. 291763). We are grateful to Alberto Racca for insightful
discussions regarding the ESN. GM is also grateful to Valerio Lucarini for
insightful discussions regarding dynamical systems theory.

## Data availability

The implementation of the ESN followsRacca2021  and the code can be found in
the github repository https://github.com/gmargazo/ESN-CLVs.git.

## Conflict of interest

The authors declare that they have no conflict of interest.

##  Appendix A Algorithms to compute LEs and CLVs

Here we present the two main algorithms for the computation of LEs and CLVs,
provided as pseudocodes. Algorithm 1 follows Pathakchaos2017 ; Vlachas2020
and contains the main steps for the calculation of $D$ Lyapunov exponents of
the echo state network. The total number of Lyapunov exponents of the ESN is
equal to the dimensionality of the hidden state, $N_{r}$, but we only compute
the first $D$, which are the physical ones. Algorithm 2 follows Ginelli2007
and computes the first $D$ Covariant Lyapunov vectors, and is the same for the
ESN and the target chaotic systems.

${\bm{U}}\leftarrow\texttt{random}\in\mathbb{R}^{N_{r}\times D}$ ;

/* Initialize $D$ GSVs */

${\bm{Q}},{\bm{R}}\leftarrow\text{QR}({\bm{U}})$ ;

/* Othonormalize GSVs */

${\bm{U}}\leftarrow{\bm{Q}}\in\mathbb{R}^{N_{r}\times D}$;

$N^{QR}\leftarrow(N^{test}-N^{w})/m$ ;

/* Number of QR decompositions */

Save the timeseries of ${\bm{R}}$ and ${\bm{Q}}$ for CLVs calculation

Initialize $\tilde{{\bm{R}}}\leftarrow\bm{0}\in\mathbb{R}^{D\times D\times
N^{QR}}$;

Initialize $\tilde{{\bm{Q}}}\leftarrow\bm{0}\in\mathbb{R}^{N_{r}\times D\times
N^{QR}}$;

Initialize $\Lambda\leftarrow\bm{0}\in\mathbb{R}^{D\times N^{QR}}$ ;

/* Save the FTLEs */

${\mathbf{W}}^{c}=\mathbf{W}^{T}_{\mathrm{in}}\mathbf{W}^{T}_{\mathrm{out}}+\mathbf{W}^{T}$
;

/* Constant matrices of Eq. (18) */

Evolve the hidden state and GSVs simultaneously.

Skip a transient initial $N^{w}$ steps for warm-up.

$n\leftarrow 0$ ;

/* Increments the number of QR decompositions */

for  _$i=0:N^{test}$_ do

${\bm{r}}({\rm{t}}_{i+1})\leftarrow\tanh\left([\hat{{\mathbf{y}}}({\rm{t}}_{i});b_{\mathrm{in}}]^{T}\mathbf{W}_{\mathrm{in}}+{\bm{r}}({\rm{t}}_{i})^{T}\mathbf{W}\right)$;

${\mathbf{y}}({\rm{t}}_{i+1})\leftarrow[{\bm{r}}({\rm{t}}_{i+1});1]^{T}\textbf{W}_{\mathrm{out}}$;

${\bm{J}}\leftarrow(1-{\bm{r}}({\rm{t}}_{i+1})^{2}){\mathbf{W}}^{c}$ ;

/* The updated Jacobian */

${\bm{U}}\leftarrow{\bm{J}}{\bm{U}}$ ;

/* The variational equation */

if  _$\mod(i,m)=0$_ then

${\bm{Q}},{\bm{R}}\leftarrow\text{QR}({\bm{U}})$ ;

/* QR every $m$ steps */

${\bm{U}}\leftarrow{\bm{Q}}$;

if  _$i >N^{w}$_ then

$\Lambda[:,n]\leftarrow\log{\left(\text{diag}[{\bm{R}}]\right)}/dt$ ;

/* Save the FTLEs */

$\tilde{{\bm{R}}}[:,:,n]\leftarrow{\bm{R}}$ ;

/* Save ${\bm{R}}$ */

$\tilde{{\bm{Q}}}[:,:,n]\leftarrow{\bm{Q}}$ ;

/* Save ${\bm{Q}}$ */

$n\leftarrow n+1$;

end if

end if

end for

$\lambda_{j}\leftarrow\sum_{i=0}^{N^{QR}}{\Lambda[j,i]}/T^{test}$ ;

/* The jth Lyapunov exponent */

Algorithm 1 Algorithm to calculate the Lyapunov exponents of the echo state
network

Data: Given $\tilde{{\bm{R}}}$ and $\tilde{{\bm{Q}}}$

Set to $\mathbf{0}$ the matrices ${\bm{C}}$, ${\mathbf{D}}$ and ${\bm{V}}$

${\bm{C}}\leftarrow\bm{0}\in\mathbb{R}^{D\times D\times N^{QR}}$ ;

/* Coordinates of CLVs in the GSV basis */

${\mathbf{D}}\leftarrow\bm{0}\in\mathbb{R}^{D\times D\times N^{QR}}$ ;

/* Growth factors of CLVs */

$\Lambda^{c}\leftarrow\bm{0}\in\mathbb{R}^{D\times N^{QR}}$ ;

/* The $D$ FTCLEs */

${\bm{V}}\leftarrow\bm{0}\in\mathbb{R}^{N_{r}\times D\times N^{QR}}$ ;

/* Each column is a CLV */

Set final time index to identity $\mathbb{I}$ for the matrices ${\bm{C}}$ and
${\mathbf{D}}$

${\bm{C}}[:,:,N^{QR}]\leftarrow\mathbb{I}$;

${\mathbf{D}}[:,:,N^{QR}]\leftarrow\mathbb{I}$;

${\bm{V}}[:,:,N^{QR}]\leftarrow\tilde{{\bm{Q}}}[:,:,N^{QR}]{\bm{C}}[:,:,N^{QR}]$;

Evolve backwards and solve Eq. (11)

for  _$i=N^{QR}-1:0$_ do

${\bm{G}}\leftarrow\texttt{solve\\_triangular}(\tilde{{\bm{R}}}[:,:,i]{\bm{C}}[:,:,i+1])$;

Normalize each column of ${\bm{G}}$

for  _$j=0:{\bm{G}}.shape[1]$_ do

${\mathbf{D}}[:,:,i]\leftarrow\texttt{norm}({\bm{G}}[:,j])$;

${\bm{C}}[:,:,i]\leftarrow{\bm{G}}[:,j]/{\mathbf{D}}[:,:,i]$;

end for

${\bm{V}}[:,:,i]\leftarrow\tilde{{\bm{Q}}}[:,:,i]{\bm{C}}[:,:,i]$;

/* Calculate CLVs */

$\Lambda^{c}[:,i]\leftarrow\log{\left(\text{diag}[{\mathbf{D}}[:,:,i]]\right)}/dt$
;

/* Calculate FTCLEs */

end for

Algorithm 2 Algorithm to calculate the covariant Lyapunov vectorsGinelli2007 ;
Ginelli2013

## References

  * (1) E.N. Lorenz, Journal of Atmospheric Sciences 20(2), 130 (1963).  DOI 10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2.  URL https://journals.ametsoc.org/view/journals/atsc/20/2/1520-0469_1963_020_0130_dnf_2_0_co_2.xml
  * (2) E. Ott, _Chaos in Dynamical Systems_ , 2nd edn. (Cambridge University Press, 2002).  DOI 10.1017/CBO9780511803260
  * (3) S.H. Strogatz, _Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering_ (CRC press, 2015).  DOI 10.1201/9780429492563.  URL https://doi.org/10.1201/9780429492563
  * (4) D. Ruelle, Annals of the New York Academy of Sciences 357(1), 1 (1980).  DOI https://doi.org/10.1111/j.1749-6632.1980.tb29669.x.  URL https://nyaspubs.onlinelibrary.wiley.com/doi/abs/10.1111/j.1749-6632.1980.tb29669.x
  * (5) J.P. Eckmann, D. Ruelle, Rev. Mod. Phys. 57, 617 (1985).  DOI 10.1103/RevModPhys.57.617.  URL https://link.aps.org/doi/10.1103/RevModPhys.57.617
  * (6) V.I. Oseledets, Trudy Moskovskogo Matematicheskogo Obshchestva 19, 179 (1968)
  * (7) G. Benettin, L. Galgani, A. Giorgilli, J.M. Strelcyn, Meccanica 15(1), 9 (1980).  DOI https://doi.org/10.1007/BF02128236
  * (8) I. Shimada, T. Nagashima, Progress of Theoretical Physics 61(6), 1605 (1979).  DOI 10.1143/PTP.61.1605.  URL https://doi.org/10.1143/PTP.61.1605
  * (9) F. Ginelli, H. Chaté, R. Livi, A. Politi, Journal of Physics A: Mathematical and Theoretical 46(25), 254005 (2013).  DOI 10.1088/1751-8113/46/25/254005.  URL https://doi.org/10.1088/1751-8113/46/25/254005
  * (10) S. Schubert, V. Lucarini, Quarterly Journal of the Royal Meteorological Society 141(693), 3040 (2015).  DOI https://doi.org/10.1002/qj.2588.  URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.2588
  * (11) S. Vannitsem, V. Lucarini, Journal of Physics A: Mathematical and Theoretical 49(22), 224001 (2016).  DOI 10.1088/1751-8113/49/22/224001.  URL https://doi.org/10.1088/1751-8113/49/22/224001
  * (12) N. Sharafi, M. Timme, S. Hallerberg, Phys. Rev. E 96, 032220 (2017).  DOI 10.1103/PhysRevE.96.032220.  URL https://link.aps.org/doi/10.1103/PhysRevE.96.032220
  * (13) E.L. Brugnago, J.A.C. Gallas, M.W. Beims, Chaos: An Interdisciplinary Journal of Nonlinear Science 30(8), 083106 (2020).  DOI 10.1063/5.0009765.  URL https://doi.org/10.1063/5.0009765
  * (14) A. Wolf, J.B. Swift, H.L. Swinney, J.A. Vastano, Physica D: Nonlinear Phenomena 16(3), 285 (1985).  DOI https://doi.org/10.1016/0167-2789(85)90011-9.  URL https://www.sciencedirect.com/science/article/pii/0167278985900119
  * (15) M.T. Rosenstein, J.J. Collins, C.J. De Luca, Physica D: Nonlinear Phenomena 65(1), 117 (1993).  DOI https://doi.org/10.1016/0167-2789(93)90009-P.  URL https://www.sciencedirect.com/science/article/pii/016727899390009P
  * (16) F. Takens, in _Dynamical Systems and Turbulence, Warwick 1980_ , ed. by D. Rand, L.S. Young (Springer Berlin Heidelberg, Berlin, Heidelberg, 1981), pp. 366–381.  DOI 10.1007/BFb0091924.  URL https://doi.org/10.1007/BFb0091924
  * (17) S. Rasp, M.S. Pritchard, P. Gentine, Proceedings of the National Academy of Sciences 115(39), 9684 (2018).  DOI 10.1073/pnas.1810286115.  URL https://www.pnas.org/doi/abs/10.1073/pnas.1810286115
  * (18) P.D. Dueben, P. Bauer, Geoscientific Model Development 11(10), 3999 (2018).  DOI 10.5194/gmd-11-3999-2018.  URL https://gmd.copernicus.org/articles/11/3999/2018/
  * (19) G. Margazoglou, T. Grafke, A. Laio, V. Lucarini, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477(2250), 20210019 (2021).  DOI 10.1098/rspa.2021.0019.  URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.2021.0019
  * (20) S. Verma, G. Novati, P. Koumoutsakos, Proceedings of the National Academy of Sciences 115(23), 5849 (2018).  DOI 10.1073/pnas.1800923115.  URL https://www.pnas.org/doi/abs/10.1073/pnas.1800923115
  * (21) D. Kochkov, J.A. Smith, A. Alieva, Q. Wang, M.P. Brenner, S. Hoyer, Proceedings of the National Academy of Sciences 118(21), e2101784118 (2021).  DOI 10.1073/pnas.2101784118.  URL https://www.pnas.org/doi/abs/10.1073/pnas.2101784118
  * (22) F. Huhn, L. Magri, Phys. Rev. Fluids 7, 014402 (2022).  DOI 10.1103/PhysRevFluids.7.014402.  URL https://link.aps.org/doi/10.1103/PhysRevFluids.7.014402
  * (23) I. Goodfellow, Y. Bengio, A. Courville, _Deep Learning_ (MIT Press, 2016).  http://www.deeplearningbook.org
  * (24) A.M. Schäfer, H.G. Zimmermann, in _Artificial Neural Networks – ICANN 2006_ , ed. by S.D. Kollias, A. Stafylopatis, W. Duch, E. Oja (Springer Berlin Heidelberg, Berlin, Heidelberg, 2006), pp. 632–640
  * (25) L. Grigoryeva, J.P. Ortega, Neural Networks 108, 495 (2018).  DOI https://doi.org/10.1016/j.neunet.2018.08.025.  URL https://www.sciencedirect.com/science/article/pii/S089360801830251X
  * (26) R.T.Q. Chen, Y. Rubanova, J. Bettencourt, D.K. Duvenaud, in _Advances in Neural Information Processing Systems_ , vol. 31, ed. by S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, R. Garnett (Curran Associates, Inc., 2018), vol. 31.  URL https://proceedings.neurips.cc/paper/2018/file/69386f6bb1dfed68692a24c8686939b9-Paper.pdf
  * (27) Z. Lu, B.R. Hunt, E. Ott, Chaos: An Interdisciplinary Journal of Nonlinear Science 28(6), 061104 (2018).  DOI 10.1063/1.5039508.  URL https://doi.org/10.1063/1.5039508
  * (28) J. Pathak, Z. Lu, B.R. Hunt, M. Girvan, E. Ott, Chaos: An Interdisciplinary Journal of Nonlinear Science 27(12), 121102 (2017).  DOI 10.1063/1.5010300.  URL https://doi.org/10.1063/1.5010300
  * (29) J. Pathak, B. Hunt, M. Girvan, Z. Lu, E. Ott, Phys. Rev. Lett. 120, 024102 (2018).  DOI 10.1103/PhysRevLett.120.024102.  URL https://link.aps.org/doi/10.1103/PhysRevLett.120.024102
  * (30) P.R. Vlachas, W. Byeon, Z.Y. Wan, T.P. Sapsis, P. Koumoutsakos, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 474(2213), 20170844 (2018).  DOI 10.1098/rspa.2017.0844.  URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.2017.0844
  * (31) P. Vlachas, J. Pathak, B. Hunt, T. Sapsis, M. Girvan, E. Ott, P. Koumoutsakos, Neural Networks 126, 191 (2020).  DOI https://doi.org/10.1016/j.neunet.2020.02.016.  URL https://www.sciencedirect.com/science/article/pii/S0893608020300708
  * (32) F. Borra, A. Vulpiani, M. Cencini, Phys. Rev. E 102, 052203 (2020).  DOI 10.1103/PhysRevE.102.052203.  URL https://link.aps.org/doi/10.1103/PhysRevE.102.052203
  * (33) N. Doan, W. Polifke, L. Magri, Journal of Computational Science 47, 101237 (2020).  DOI https://doi.org/10.1016/j.jocs.2020.101237.  URL https://www.sciencedirect.com/science/article/pii/S1877750320305408
  * (34) N.A.K. Doan, W. Polifke, L. Magri, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477(2253), 20210135 (2021).  DOI 10.1098/rspa.2021.0135.  URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.2021.0135
  * (35) A. Racca, L. Magri, Neural Networks 142, 252 (2021).  DOI https://doi.org/10.1016/j.neunet.2021.05.004.  URL https://www.sciencedirect.com/science/article/pii/S0893608021001969
  * (36) Y. Kuramoto, Progress of Theoretical Physics Supplement 64, 346 (1978).  DOI 10.1143/PTPS.64.346.  URL https://doi.org/10.1143/PTPS.64.346
  * (37) M. Lukoševičius, _A Practical Guide to Applying Echo State Networks_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012), pp. 659–686.  DOI 10.1007/978-3-642-35289-8_36.  URL https://doi.org/10.1007/978-3-642-35289-8_36
  * (38) O. Rössler, Physics Letters A 57(5), 397 (1976).  DOI https://doi.org/10.1016/0375-9601(76)90101-8.  URL https://www.sciencedirect.com/science/article/pii/0375960176901018
  * (39) J.G. Charney, J.G. DeVore, Journal of Atmospheric Sciences 36(7), 1205 (1979).  DOI 10.1175/1520-0469(1979)036<1205:MFEITA>2.0.CO;2.  URL https://journals.ametsoc.org/view/journals/atsc/36/7/1520-0469_1979_036_1205_mfeita_2_0_co_2.xml
  * (40) E. Lorenz, in _Seminar on Predictability, 4-8 September 1995_ , vol. 1. ECMWF (ECMWF, Shinfield Park, Reading, 1995), vol. 1, pp. 1–18.  URL https://www.ecmwf.int/node/10829
  * (41) F. Ginelli, P. Poggi, A. Turchi, H. Chaté, R. Livi, A. Politi, Phys. Rev. Lett. 99, 130601 (2007).  DOI 10.1103/PhysRevLett.99.130601.  URL https://link.aps.org/doi/10.1103/PhysRevLett.99.130601
  * (42) G. Boffetta, M. Cencini, M. Falcioni, A. Vulpiani, Physics Reports 356(6), 367 (2002).  DOI https://doi.org/10.1016/S0370-1573(01)00025-4.  URL https://www.sciencedirect.com/science/article/pii/S0370157301000254
  * (43) F. Huhn, L. Magri, Journal of Fluid Mechanics 882, A24 (2020).  DOI 10.1017/jfm.2019.828
  * (44) D. Ruelle, Communications in Mathematical Physics 187(1), 227 (1997).  DOI 10.1007/s002200050134.  URL https://doi.org/10.1007/s002200050134
  * (45) V. Lucarini, R. Blender, C. Herbert, F. Ragone, S. Pascale, J. Wouters, Reviews of Geophysics 52(4), 809 (2014).  DOI https://doi.org/10.1002/2013RG000446.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2013RG000446
  * (46) P.J. Blonigan, P. Fernandez, S.M. Murman, Q. Wang, G. Rigas, L. Magri, (Center for Turbulence Research Proceedings of the Summer Program, 2017).  DOI 10.17863/CAM.35422.  URL https://doi.org/10.17863/CAM.35422
  * (47) C.L. Wormell, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 478(2261), 20210808 (2022).  DOI 10.1098/rspa.2021.0808.  URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.2021.0808
  * (48) S.P. Kuznetsov, _Possible Occurrence of Hyperbolic Attractors_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012), pp. 35–56.  DOI 10.1007/978-3-642-23666-2_2.  URL https://doi.org/10.1007/978-3-642-23666-2_2
  * (49) G. Gallavotti, E.G.D. Cohen, Phys. Rev. Lett. 74, 2694 (1995).  DOI 10.1103/PhysRevLett.74.2694.  URL https://link.aps.org/doi/10.1103/PhysRevLett.74.2694
  * (50) G. Gallavotti, E.G.D. Cohen, Journal of Statistical Physics 80(5), 931 (1995)
  * (51) J.L. Lebowitz, H. Spohn, Journal of Statistical Physics 95(1), 333 (1999).  DOI 10.1023/A:1004589714161.  URL https://doi.org/10.1023/A:1004589714161
  * (52) S. Lepri, R. Livi, A. Politi, Physica D: Nonlinear Phenomena 119(1), 140 (1998).  DOI https://doi.org/10.1016/S0167-2789(98)00076-1.  URL https://www.sciencedirect.com/science/article/pii/S0167278998000761.  Localization in Nonlinear Lattices
  * (53) H.l. Yang, K.A. Takeuchi, F. Ginelli, H. Chaté, G. Radons, Phys. Rev. Lett. 102, 074102 (2009).  DOI 10.1103/PhysRevLett.102.074102.  URL https://link.aps.org/doi/10.1103/PhysRevLett.102.074102
  * (54) A. Viennet, N. Vercauteren, M. Engel, D. Faranda.  Guidelines for data-driven approaches to study transitions in multiscale systems: the case of lyapunov vectors (2022).  DOI 10.48550/ARXIV.2203.10322.  URL https://arxiv.org/abs/2203.10322
  * (55) C. Martin, N. Sharafi, S. Hallerberg, Chaos: An Interdisciplinary Journal of Nonlinear Science 32(3), 033105 (2022).  DOI 10.1063/5.0078112.  URL https://doi.org/10.1063/5.0078112
  * (56) P. Virtanen, R. Gommers, T.E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S.J. van der Walt, M. Brett, J. Wilson, K.J. Millman, N. Mayorov, A.R.J. Nelson, E. Jones, R. Kern, E. Larson, C.J. Carey, İ. Polat, Y. Feng, E.W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E.A. Quintero, C.R. Harris, A.M. Archibald, A.H. Ribeiro, F. Pedregosa, P. van Mulbregt, SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).  DOI 10.1038/s41592-019-0686-2
  * (57) P.V. Kuptsov, S.P. Kuznetsov, Regular and Chaotic Dynamics 23(7), 908 (2018).  DOI 10.1134/S1560354718070079.  URL https://doi.org/10.1134/S1560354718070079
  * (58) S. Hochreiter, J. Schmidhuber, Neural Computation 9(8), 1735 (1997).  DOI 10.1162/neco.1997.9.8.1735.  URL https://doi.org/10.1162/neco.1997.9.8.1735
  * (59) K. Cho, B. van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, Y. Bengio, in _Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ (Association for Computational Linguistics, Doha, Qatar, 2014), pp. 1724–1734.  DOI 10.3115/v1/D14-1179.  URL https://aclanthology.org/D14-1179
  * (60) P. Werbos, Proceedings of the IEEE 78(10), 1550 (1990).  DOI 10.1109/5.58337
  * (61) H. Jaeger, H. Haas, Science 304(5667), 78 (2004).  DOI 10.1126/science.1091277.  URL https://www.science.org/doi/abs/10.1126/science.1091277
  * (62) F. Huhn, L. Magri, in _Computational Science – ICCS 2020_ , ed. by V.V. Krzhizhanovskaya, G. Závodszky, M.H. Lees, J.J. Dongarra, P.M.A. Sloot, S. Brissos, J. Teixeira (Springer International Publishing, Cham, 2020), pp. 124–132
  * (63) A. Racca, L. Magri.  Data-driven prediction and control of extreme events in a chaotic flow (2022).  DOI 10.48550/ARXIV.2204.11682.  URL https://arxiv.org/abs/2204.11682
  * (64) A.N. Tikhonov, A. Goncharsky, V. Stepanov, A.G. Yagola, _Numerical methods for the solution of ill-posed problems_ , vol. 328 (Springer Science & Business Media, 1995).  DOI doi.org/10.1007/978-94-015-8480-7.  URL https://doi.org/10.1007/978-94-015-8480-7
  * (65) M. Hoffman, E. Brochu, N. de Freitas, in _Proceedings of the Twenty-Seventh Conference on Uncertainty in Artificial Intelligence_ (AUAI Press, Arlington, Virginia, USA, 2011), UAI’11, p. 327–336
  * (66) M. Lukoševičius, A. Uselis, in _Artificial Neural Networks and Machine Learning – ICANN 2019: Workshop and Special Sessions_ , ed. by I.V. Tetko, V. Kůrková, P. Karpov, F. Theis (Springer International Publishing, Cham, 2019), pp. 121–133
  * (67) H. De Swart, Physica D: Nonlinear Phenomena 36(3), 222 (1989).  DOI https://doi.org/10.1016/0167-2789(89)90082-1.  URL https://www.sciencedirect.com/science/article/pii/0167278989900821
  * (68) D.T. Crommelin, J.D. Opsteegh, F. Verhulst, Journal of the Atmospheric Sciences 61(12), 1406 (2004).  DOI 10.1175/1520-0469(2004)061<1406:AMFARB>2.0.CO;2.  URL https://journals.ametsoc.org/view/journals/atsc/61/12/1520-0469_2004_061_1406_amfarb_2.0.co_2.xml
  * (69) T. Woollings, D. Barriopedro, J. Methven, S.W. Son, O. Martius, B. Harvey, J. Sillmann, A.R. Lupo, S. Seneviratne, Current climate change reports 4(3), 287 (2018).  DOI https://doi.org/10.1007/s40641-018-0108-z
  * (70) S. Schubert, V. Lucarini, Quarterly Journal of the Royal Meteorological Society 142(698), 2143 (2016).  DOI https://doi.org/10.1002/qj.2808.  URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.2808
  * (71) D. Ruelle, Communications in Mathematical Physics 87(2), 287 (1982).  DOI https://doi.org/10.1007/BF01218566
  * (72) A. Karimi, M.R. Paul, Chaos: An Interdisciplinary Journal of Nonlinear Science 20(4), 043105 (2010).  DOI 10.1063/1.3496397.  URL https://doi.org/10.1063/1.3496397
  * (73) P.V. Kuptsov, S.P. Kuznetsov, Phys. Rev. E 80, 016205 (2009).  DOI 10.1103/PhysRevE.80.016205.  URL https://link.aps.org/doi/10.1103/PhysRevE.80.016205
  * (74) A.V. Knyazev, M.E. Argentati, SIAM Journal on Scientific Computing 23(6), 2008 (2002).  DOI 10.1137/S1064827500377332.  URL https://doi.org/10.1137/S1064827500377332
  * (75) J.L. Kaplan, J.A. Yorke, in _Functional differential equations and approximation of fixed points_ (Springer, 1979), pp. 204–227