# Reconstruction, forecasting, and stability of chaotic dynamics from partial
data

Elise Özalp  Georgios Margazoglou  Department of Aeronautics, Imperial College
London, London SW7 2BX, United Kingdom  Luca Magri<EMAIL_ADDRESS>Department of Aeronautics, Imperial College London, London SW7 2BX, United
Kingdom  The Alan Turing Institute, London NW1 2DB, United Kingdom

###### Abstract

The forecasting and computation of the stability of chaotic systems from
partial observations are tasks for which traditional equation-based methods
may not be suitable. In this computational paper, we propose data-driven
methods to (i) infer the dynamics of unobserved (hidden) chaotic variables
(full-state reconstruction); (ii) time forecast the evolution of the full
state; and (iii) infer the stability properties of the full state. The tasks
are performed with long short-term memory (LSTM) networks, which are trained
with observations (data) limited to only part of the state: (i) the low-to-
high resolution LSTM (LH-LSTM), which takes partial observations as training
input, and requires access to the full system state when computing the loss;
and (ii) the physics-informed LSTM (PI-LSTM), which is designed to combine
partial observations with the integral formulation of the dynamical system’s
evolution equations. First, we derive the Jacobian of the LSTMs. Second, we
analyse a chaotic partial differential equation, the Kuramoto-Sivashinsky
(KS), and the Lorenz-96 system. We show that the proposed networks can
forecast the hidden variables, both time-accurately and statistically. The
Lyapunov exponents and covariant Lyapunov vectors, which characterize the
stability of the chaotic attractors, are correctly inferred from partial
observations. Third, the PI-LSTM outperforms the LH-LSTM by successfully
reconstructing the hidden chaotic dynamics when the input dimension is smaller
or similar to the Kaplan-Yorke dimension of the attractor. This work opens new
opportunities for reconstructing the full state, inferring hidden variables,
and computing the stability of chaotic systems from partial data.

††preprint: AIP/123-QED

> It has been shown that the chaotic dynamics of complex systems can be
> inferred by machine learning from data, but questions arise about the
> physical relevance of the surrogate models, especially when only partial
> observations on the system are available. We design data-driven approaches
> that can accurately infer hidden chaotic dynamics form partial data, and
> correctly learn the chaotic attractor and its stability properties.

##  I Introduction

Chaotic dynamics appear in applications in various fields, from meteorology
Lorenz (1963), through engineering and chemistry Strogatz (1994), to
propulsion Huhn and Magri (2020); Magri and Doan (2020). Chaos emerges due to
the system’s exponential sensitivity to initial conditions, as is the case
with turbulent fluid dynamics Takens (1981). However, the accurate prediction
of chaotic spatiotemporal behaviour is challenging because numerical
computations can be computationally expensive, and sensor measurements are
often limited, providing only a partial observation of the dynamic behaviour.
These partial observations can lead to a misrepresentation of the unobserved
(hidden) variables, or other long-term physical properties, which further
challenges the modelling and prediction of the systems.

The predictability and stability of a chaotic system are characterized by its
tangent space, which can be computed using the linearized dynamics provided by
the Jacobian. This computation allows for the derivation of quantities such as
the Lyapunov exponents (LEs), which measure the exponential rate of separation
of trajectories. A geometric characterization is provided by the covariant
Lyapunov vectors (CLVs), which constitute a covariant basis of the tangent
space, and point to directions of asymptotic expansion and contraction of the
dynamical system Ginelli  _et al._ (2013). Preserving these stability
properties is crucial when building surrogate models from limited observations
to a more comprehensive dataset Takeuchi  _et al._ (2011).

Neural networks are an expressive nonlinear representation of continuous
functions, which can extract patterns from data and, once trained, provide
computationally cheap predictions. Suitable for time series and dynamical
evolutions are recurrent neural networks (RNNs), which have shown promising
performance in the inference of dynamical systems with multi-scale, chaotic,
or turbulent behaviour Vlachas  _et al._ (2018). There are several classes of
RNNs, the main classes being reservoir computing, such as echo state networks
(ESNs) Jaeger (2001), and networks using backpropagation through time, such as
long short-term memory networks (LSTMs) Hochreiter and Schmidhuber (1997a).
Previous works compared these two networks for the short- and long-term
prediction horizons of chaotic systems Vlachas  _et al._ (2020); Racca and
Magri (2021a). ESNs have been shown to accurately learn (i.e. infer) the
stability properties for prototypical chaotic systems Margazoglou and Magri
(2023), but when the full state is available for training.

This work proposes two approaches to infer the hidden chaotic dynamics from
partial observations (data) with LSTMs. The first approach considers the
scenario in which the full-state data is only available for a limited period,
and partial data is available for the remaining time. This scenario arises
when generating high-resolution data is computationally expensive. For this,
the goal is to analyse how well LSTMs can infer the full-state dynamics from
partial data, especially in capturing small-scale statistics in partial-
differential equations. The second approach analyses the case in which only
partial state data is available, and the aim is to reconstruct (i.e., infer)
the hidden dynamics. To achieve this, knowledge of the underlying physical
system has to be incorporated by adding constraints based on the governing
equations. We propose a physics-informed loss function based on an integral
formulation of the dynamical system, which is versatile and suitable for the
architecture of the RNN.

The paper aims to achieve two primary objectives. First, it aims to accurately
infer the physical properties of the chaotic systems by comparing the long-
term statistics of the degrees of freedom, the LEs, angles of the CLVs, and
energy of the system. Secondly, both approaches are compared on partial inputs
from two prototypical chaotic systems: the Lorenz-96 system Lorenz (1996) and
the Kuramoto–Sivashinsky equation Kuramoto (1978); Hyman and Nicolaenko
(1986).

The paper is structured as follows. Section II introduces the LSTM
architecture to address the state reconstruction problem of a chaotic
dynamical system. We propose a data-driven and physics-informed loss
formulation in Section III. In Section IV, we provide an overview of the
computation of LEs and CLVs, emphasizing their importance for physically
accurate modelling. We derive the Jacobian of the LSTM. Section V discusses
the results for the Kuramoto–Sivashinsky system, followed by the findings for
the Lorenz-96 system. Finally, in Section VI, we summarize our work and
discuss future directions.

##  II Long short-term memory

A dynamical system’s solution is a time series, which provides an ordered
sequence of data over time. Recurrent neural networks (RNNs) are a data-driven
approach for modelling time series that utilize hidden states to encode the
input’s information history. Among RNNs, long short-term memory (LSTM)
networks are a common tool due to their internal gating mechanisms, which
mitigate the vanishing gradient problem Hochreiter and Schmidhuber (1997b).

LSTMs are commonly employed for time series forecasting as they embed the
time-delayed inputs in their higher-dimensional hidden states. From a
dynamical systems perspective, a delay embedding can provide a structure that
is topologically equivalent to the attractor Takens (1981). Given their
architecture, LSTMs have the potential to be a powerful tool for representing
the dynamics of a chaotic system. This enables them to be utilised for
reconstructing and forecasting dynamical systems, even in cases in which we
lack part of the system’s full state information, as it will be demonstrated
in the following.

###  II.1 State reconstruction for the inference of hidden dynamics

We consider a nonlinear autonomous dynamical system

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

where $\bm{x}(t)\in\mathbb{R}^{D}$ is the state vector of the physical system
and $f:\mathbb{R}^{D}\to\mathbb{R}^{D}$ is a smooth nonlinear vector function.
Experimental observations often come with a limited amount of information on a
system’s full state. For example, if in a wind tunnel experiment
$\mathcal{O}(10)$ sensors are placed in different positions, they may provide
adequate but scarce spatial information of the system’s evolution Adrian
(1991). Mathematically, if $\bm{x}(t)=[\bm{y}(t);\bm{\xi}(t)]$ is the full
state of a chaotic dynamical system, then $\bm{y}(t)\in\mathbb{R}^{D_{y}}$ are
the scarce observations and $\bm{\xi}(t)\in\mathbb{R}^{D_{\xi}}$ are the
unobserved (hidden) variables with $D=D_{y}+D_{\xi}$. Specifically, let us
assume that $\bm{y}(t_{i})$ is measured at times $t_{i}=i\Delta t$ with
$i=0,\dots N_{t}$ and constant time step $\Delta t$. Based on these
observations, we wish to reconstruct and predict the full state
$\bm{x}(t_{i})=[\bm{y}(t_{i}),\bm{\xi}(t_{i})]$.

LSTMs have been successfully applied to time series forecasting of dynamical
systems when full observations are available Sangiorgio, Dercole, and Guariso
(2021); Vlachas  _et al._ (2018, 2020). In this work, we use the LSTM as a
model of partially observed chaotic timeseries to both reconstruct the full
state and perform accurate autonomous temporal evolution (forecasting).

###  II.2 Architecture

The long short-term memory networks are characterized by a cell state
$\bm{c}_{i}\in\mathbb{R}^{N_{h}}$ and a hidden state
$\bm{h}_{i}\in\mathbb{R}^{N_{h}}$ that are updated at each step. In the case
of partial observations, the states are updated by using the observed
variables $\bm{y}(t_{i})$ as

$\displaystyle\bm{i}_{i+1}$
$\displaystyle=\sigma\left(\bm{W}^{i}[\bm{y}(t_{i});\bm{h}_{i}]+\bm{b}^{i}\right),$
(2) $\displaystyle\bm{f}_{i+1}$
$\displaystyle=\sigma\left(\bm{W}^{f}[\bm{y}(t_{i});\bm{h}_{i}]+\bm{b}^{f}\right),$
(3) $\displaystyle\bm{o}_{i+1}$
$\displaystyle=\sigma\left(\bm{W}^{o}[\bm{y}(t_{i});\bm{h}_{i}]+\bm{b}^{o}\right),$
(4) $\displaystyle\bm{\tilde{c}}_{i+1}$
$\displaystyle=\tanh{\left(\bm{W}^{\bm{g}}[\bm{y}(t_{i});\bm{h}_{i}]+\bm{b}^{\bm{g}}\right)},$
$\displaystyle\bm{c}_{i+1}$
$\displaystyle=\bm{f}_{i+1}*\bm{c}_{i}+\bm{i}_{i+1}*\bm{\tilde{c}}_{i+1},$ (5)
$\displaystyle\bm{h}_{i+1}$
$\displaystyle=\tanh{\left(\bm{c}_{i+1}\right)}*\bm{o}_{i+1},$

where $*$ denotes the elementwise multiplication. In the LSTM,
$\bm{i}_{i+1},\bm{f}_{i+1},\bm{o}_{i+1}\in\mathbb{R}^{N_{h}}$ are the input,
forget and output gates (Fig. 1). The matrices $\bm{W}^{i},$ $\bm{W}^{f},$
$\bm{W}^{o},$ $\bm{W}^{\bm{g}}\in\mathbb{R}^{N_{h}\times(D+N_{h})}$ are the
corresponding weight matrices, and $\bm{b}^{i},$ $\bm{b}^{f},$ $\bm{b}^{o},$
$\bm{b}^{\bm{g}}\in\mathbb{R}^{N_{h}}$ are the biases.

$\sigma$$\sigma$Tanh$\sigma$$\times$+$\times$$\times$Tanh$c_{i}$previous cell
state$h_{i}$previous hidden state$y(t_{i})$input$c_{i+1}$new cell
state$h_{i+1}$new hidden
state$h_{i+1}$output$\begin{bmatrix}\tilde{y}(t_{i+1})\\\
\tilde{\xi}(t_{i+1})\end{bmatrix}$prediction

Figure 1: Schematic representation of the LSTM cell structure

The full prediction on the next time step,
$\bm{\hat{x}}(t_{i+1})=[\bm{\hat{y}}(t_{i+1}),\bm{\hat{\xi}}(t_{i+1})]$ is
obtained by applying a dense layer to the hidden state $\bm{h}_{i+1}$, i.e.

$\displaystyle\begin{bmatrix}\bm{\hat{y}}(t_{i+1})\\\
\bm{\hat{\xi}}(t_{i+1})\end{bmatrix}=\bm{W}^{dense}\bm{h}_{i+1}+\bm{b}^{dense},$

where $\bm{W}^{dense}\in\mathbb{R}^{D\times N_{h}}$ and
$\bm{b}^{dense}\in\mathbb{R}^{D}$.

$\mathbf{LSTM}$$\bm{y}(t_{0})$Dense$\bm{h}_{1}$$\mathbf{LSTM}$$\bm{y}(t_{1})$Dense$\bm{h}_{2}$$\begin{bmatrix}\bm{\tilde{y}}(t_{1})\\\
\bm{\tilde{\xi}}(t_{1})\end{bmatrix}$ $\begin{bmatrix}\bm{\tilde{y}}(t_{2})\\\
\bm{\tilde{\xi}}(t_{2})\end{bmatrix}$ …$\mathbf{LSTM}$
$\bm{y}(t_{n-1})$…$\bm{c}_{1},\bm{h}_{1}$$\bm{c}_{2},\bm{h}_{2}$$\bm{c}_{n-1},\bm{h}_{n-1}$Dense
$\bm{h}_{n}$$\begin{bmatrix}\bm{\tilde{y}}(t_{n})\\\
\bm{\tilde{\xi}}(t_{n})\end{bmatrix}$ $\bm{c}_{0},\bm{h}_{0}$

Figure 2: LSTM in open-loop configuration: Each cell receives input from the
training data.

LSTMs are universal approximators for a continuous target functionHornik,
Stinchcombe, and White (1989); however, practically, the network’s performance
depends on the parameters, such as weights and biases. To determine these
parameters, the available data is divided into three subsets: training,
validation, and testing. The training and validation data is split into
batches of fixed-length time windows. The training data is used for
backpropagation through time Werbos (1990), and during the forward pass, the
network output is computed and compared to a training label using a loss
function. In the backward pass, the network parameters are optimized by
computing the gradient of the loss function with respect to the parameters
with a gradient update via the Adam optimizer Kingma and Ba (2014). The
validation data is employed to optimize hyperparameters and determine
parameters before training (e.g. the dimension of the hidden and cell state
$N^{h}$), while test data is used to evaluate the final performance of the
model.

The network works in an open-loop configuration during training and validation
(Fig. 2) with the output at each time step depending on current and previous
inputs within the time window, and the LSTM’s states are reset at the start of
each input sequence. This is also known as "teacher-forced learning" Gers,
Schraudolph, and Schmidhuber (2002), which facilitates a straightforward
formulation of loss through backpropagation, as elaborated in Section III.
After training, the network is evaluated on test data with fixed weights and
biases, operating in a closed-loop configuration (Fig. 3). In this mode, the
network can make long-term predictions even if no data is available, and the
states of the LSTM are retained across time steps. The network predicts both
the observed and unobserved (hidden) variables. Subsequently, the observed
variables are provided as inputs for the next time step, which allows the
autonomous evolution of the LSTM.

##  III Tasks and loss functions

The LSTM’s internal dynamics, and therefore its prediction, are determined by
the weights and biases of the model, which are updated during the training.
The loss function plays a critical role in this process by defining the
optimization problem that guides the updating of these parameters. The choice
of the loss function depends on the specific application and can greatly
impact the accuracy and generalization performance of the model. In the
following, we propose two loss functions to guide the LSTM training: (i) a
data-driven loss, and (ii) a physics-informed loss.

$\mathbf{LSTM}$ $\bm{y}(t_{n})$$\bm{c}_{n},\bm{h}_{n}$Dense
$\bm{h}_{n+1}$$\begin{bmatrix}\bm{\tilde{y}}(t_{n+1})\\\
\bm{\tilde{\xi}}(t_{n+1})\end{bmatrix}$ $\mathbf{LSTM}$
$\tilde{\bm{y}}(t_{n+1})$$\bm{c}_{n+1},\bm{h}_{n+1}$Dense
$\bm{h}_{n+2}$$\begin{bmatrix}\bm{\tilde{y}}(t_{n+2})\\\
\bm{\tilde{\xi}}(t_{n+2})\end{bmatrix}$ $\mathbf{LSTM}$
$\tilde{\bm{y}}(t_{n+2})$$\bm{c}_{n+2},\bm{h}_{n+2}$Dense
$\bm{h}_{n+3}$$\begin{bmatrix}\bm{\tilde{y}}(t_{n+3})\\\
\bm{\tilde{\xi}}(t_{n+3})\end{bmatrix}$ $\bm{c}_{n+3},\bm{h}_{n+3}$

Figure 3: LSTM in closed-loop configuration: The network prediction is used as
input for the next cell.

###  III.1 Low-to-high resolution loss (LH-LSTM)

The first task is to infer the high-resolution dynamics from low-resolution
inputs, thereby reducing the cost of acquiring full measurement data. In fluid
dynamics, in which problems display turbulent behaviour, this task becomes
especially difficult as the prediction is sensitive due to the chaotic nature
of the data. The network predicts the full state
$\bm{\bm{\hat{x}}}(t_{i})=[\bm{\bm{\hat{y}}}(t_{i});\bm{\hat{\xi}}(t_{i})]$
given a partial input $\bm{y}(t_{i})$. In this section, we assume that the
full label $\bm{x}(t_{i})=[\bm{y}(t_{i});\bm{\xi}(t_{i})]$ is only available
for a short amount of time, which can be employed for training and validation,
and partial data is available for the remaining time. Hence, we can define a
data-driven loss using a mean-squared error (MSE)

$\displaystyle\mathcal{L}_{dd}(\bm{x},\bm{\hat{x}})=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\|\bm{x}(t_{i})-\bm{\hat{x}}(t_{i})\|^{2},$
(6)

where $||\cdot||$ is the $L^{2}$-norm. Although simple to implement, data-
driven loss functions have limitations. Specifically, they are sensitive to
overfitting, especially when the training data is noisy. When data is limited,
these loss functions may struggle to generalize to regions of the attractor
that are sparsely sampled during training. To regularize the problem, the
incorporation of the physical knowledge or constraints into the loss function
is an effective means Doan, Polifke, and Magri (2021); Özalp, Margazoglou, and
Magri (2023), which favour physical solutions (in contrast to traditional
Tikhonov regularization).

###  III.2 Physics-informed loss (PI-LSTM)

Physics-informed losses have shown success in feed-forward neural networks, in
which automatic differentiation can be exploited to accurately compute the
derivative with respect to time and space to machine precision Lagaris, Likas,
and Fotiadis (1998); Raissi, Perdikaris, and Karniadakis (2019a). Unlike feed-
forward neural networks, RNNs have a dynamic temporal structure that makes it
difficult to accurately differentiate the prediction with respect to time
because of the recurrence and the type of inputs. As a result, the governing
equations have to be incorporated in a manner that accounts for the temporal
structure of the RNN architecture. This approximation is typically achieved by
discretizing the time derivative in the loss function using finite
differences, which allows the physics-informed loss to be computed but it
relies on the numerical differentiation scheme Özalp, Margazoglou, and Magri
(2023). Hence, the autonomous dynamical system defined by Eq. (1) is
reformulated with its formal solution through the integral

$\displaystyle\bm{x}(t_{i+1})$
$\displaystyle=\int_{t_{0}}^{t_{i+1}}f(\bm{x}(t))dt\
=\bm{x}(t_{i})+\int_{t_{i}}^{t_{i+1}}f(\bm{x}(t))dt,\ i\geq 0,$ (7)

which enables the use of numerical quadrature methods to approximate the
integral $\int_{t_{i}}^{t_{i+1}}f(\bm{x}(t))dt$, instead of approximating the
derivative $\frac{d\bm{x}(t)}{dt}$ to construct an accurate surrogate model.
This approach is naturally compatible with explicit numerical schemes of
different orders of accuracy, including the explicit Runge-Kutta family or
pseudo-spectral methods, thereby providing a flexible framework for modelling
physical constraints into RNNs. Based on Eq. (7), we define the residual of
the dynamical system

$\displaystyle\mathcal{R}(\bm{x}(t_{i+1}))=\bm{x}(t_{i+1})-\left(\bm{x}(t_{i})+\int_{t_{i}}^{t_{i+1}}f(\bm{x}(t))dt\right).$
(8)

The true solution, $\bm{x}(t)$, of the dynamical system, is such that
$\mathcal{R}(\bm{x}(t_{i}))=0$, for all $t_{i}>t_{0}$. In this manner, the
predictions $\bm{\hat{x}}$ of the network favour physical solutions by
minimizing the physics-informed loss

$\displaystyle\mathcal{L}_{pi}(\bm{\hat{x}})=\frac{1}{N_{t}}\sum_{i=0}^{N_{t}-1}\|\mathcal{R}(\hat{\bm{x}}(t_{i+1}))\|^{2}.$
(9)

The proposed loss function addresses a distinct task from the data-driven loss
(6) by not requiring a full label, $\bm{x}(t_{i})$, making it particularly
advantageous in situations in which full labels are not available. Instead, it
enables the network’s prediction to be constrained based on the governing
equations, which allows the reconstruction and the forecasting of unobserved
(hidden) variables. By providing a regularization during the training of the
LSTM, we also ensure that the predictions satisfy the governing equation
defined in (1) within a numerical tolerance. Further, this regularization
helps to constrain the parameter space of the weights, which can mitigate
overfitting and improve generalization. We also include a data-driven loss for
the available data by computing the MSE between the network prediction on the
partial input $\bm{\hat{y}}(t_{i+1})$ and the partial label $\bm{y}(t_{i+1})$

$\displaystyle\mathcal{L}_{dd}(\bm{y},\bm{\hat{y}})=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\|\bm{y}(t_{i})-\bm{\hat{y}}(t_{i})\|^{2}.$
(10)

The loss of the physics-informed LSTM is computed by combining the data-driven
loss and weighing the physics-informed loss as

$\mathcal{L}(\bm{y},\bm{\hat{x}})=\mathcal{L}_{dd}(\bm{y},\bm{\hat{y}})+\alpha_{pi}\mathcal{L}_{pi}(\bm{\hat{x}}),\quad\alpha_{pi}\in\mathbb{R}^{+},$
(11)

where $\alpha_{pi}$ is a penalty hyperparameter, which acts as a
regularization factor. The selection of an appropriate hyperparameter for the
regularization term is a nontrivial task. For this, a grid search is employed.

##  IV Inferring the stability of the dynamical system from data

Thanks to its recurrent nature, the LSTM can use its own prediction as input
in closed-loop mode, which allows it to predict the system’s full state in
time beyond the input. This autonomous evolution of LSTMs defines a dynamical
system. Mathematically, we assess the LSTM stability properties by studying
their tangent space Margazoglou and Magri (2023). By introducing small
perturbations to the system’s trajectory and linearizing the system’s
equations, we can compute its stability in different directions. We determine
the properties of its linear tangent space, such as Lyapunov exponents (LEs)
and covariant Lyapunov vectors (CLVs). To gain insight into the network’s
dynamical behaviour, we apply these concepts to the LSTM by mathematically
deriving the LSTM Jacobian.

###  IV.1 Lyapunov exponents

The solution of the dynamical system in Eq. (1) is $\bm{x}(t)$, which defines
a trajectory in the phase space of dimension $D$. To analyse the properties of
the dynamical system and its attractor, we impose infinitesimal perturbations
to the trajectory of the system as

$\bm{x}+\bm{u},\quad\bm{x}\sim\mathcal{O}(1),\quad\bm{u}\sim\mathcal{O}(\epsilon),\quad\epsilon\to
0.$ (12)

By substituting (12) into (1) and linearizing the system, the perturbation is
governed by the tangent equation

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

where $J_{ij}=\frac{\partial f_{i}(\bm{x})}{\partial x_{j}}$ are the
components of the Jacobian $\bm{J}(\bm{x}(t))\in\mathbb{R}^{D\times D}$, which
is time-dependent on chaotic attractors. This shows that the perturbation
$\bm{u}(t)$ evolves along the tangent space with respect to $\bm{x}(t)$, and
must therefore be integrated alongside Eq. (1). Considering $K$ random
perturbations, we numerically integrate $K\leq D$ tangent vectors
$\bm{u}_{i}\in\mathbb{R}^{D}$ as columns of $\bm{U}\in\mathbb{R}^{D\times K}$

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

In chaotic systems, the tangent vectors align exponentially fast with the
leading Lyapunov vector, which can lead to an ill-conditioned matrix $\bm{U}$
Ginelli  _et al._ (2013). For practical purposes, $\bm{U}$ is periodically
orthonormalized through the Gram-Schmidt procedure Benettin  _et al._ (1980).
Hence, we decompose $\bm{U}$ into an orthogonal matrix $\bm{Q}$ and an upper-
triangular matrix $\bm{R}$ with QR-decomposition, i.e.
$\bm{U}(t)=\bm{Q}(t)\bm{R}(t)$. The columns of $\bm{Q}$ form an orthonormal
basis for the tangent space and are known as Gram-Schmidt vectors (GSVs). The
diagonal entries $R_{ii}(t)$ correspond to the local growth rates of the
corresponding GSV $\bm{q}_{i}(t)$. The Lyapunov exponents (LEs) are computed
by taking the time average of the logarithms of $R_{ii}(t)$

$\lambda_{i}=\lim_{T\to\infty}\frac{1}{T-t_{0}}\int_{t_{0}}^{T}\ln\left(R_{ii}(t)\right).$
(15)

The Lyapunov spectrum, $\lambda_{1}\geq\dots\geq\lambda_{D}$, provides
fundamental insight into the chaotic properties and geometry of an attractor
Ginelli  _et al._ (2013); Huhn and Magri (2020). If the leading Lyapunov
exponent $\lambda_{1}<0$, the perturbation decay and the attractor is a fixed
point. If $\lambda_{1}=0$ and the remaining exponents are negative, the
attractor is a periodic orbit. If $\lambda_{1}>0$, the perturbation grows
exponentially and, typically, the attractor is chaotic. In this case, the
Lyapunov time $\tau_{\lambda}=\frac{1}{\lambda_{1}}$ defines a characteristic
timescale for two nearby trajectories to separate, which gives an estimate of
the system’s predictability horizon Boffetta  _et al._ (2002); Magri and Doan
(2020). Further, the LEs provide an estimate for the attractor dimension
through the Kaplan-Yorke dimension Frederickson  _et al._ (1983); Kantz and
Schreiber (2003)

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

with $\sum_{i=1}^{k}\lambda_{i}>0$ and $\sum_{i=1}^{k+1}\lambda_{i}<0$.

###  IV.2 Covariant Lyapunov vectors

Because LEs are scalars, they cannot describe the geometry of the tangent
space. For this, we need to find a suitable basis that spans the tangent
space. The GSVs $\bm{Q}(t)$ provide an orthogonal basis of the tangent space,
which is not time-reversible because of the frequent orthogonalizations. On
the other hand, there exists a coordinate-independent, local decomposition of
the phase space into subspaces spanned by the covariant Lyapunov vectors
(CLVs) Oseledec (1968); Ruelle (1979). Unlike the GSVs, which are orthonormal
by construction, the subspaces spanned by the CLVs are generally non-
orthogonal to each other, which provide information on the local geometric
structure of the chaotic attractor. If the CLVs
$\bm{V}=\begin{bmatrix}\bm{v}_{1},\bm{v}_{2},\dots,\bm{v}_{k}\end{bmatrix}$
are uniquely defined (i.e. non degenerate), then each CLV
$\bm{v}_{i}\in\mathbb{R}^{D}$ describes the individual expansion and
contraction associated with the LE $\lambda_{i}$. The CLVs can be recovered
from the local growth rates $\bm{R}(t)$ and the GSVs $\bm{Q}(t)$ by evolving
backwards in time after the forward-in-time simulation is completed. We
provide a brief overview on the computation of the CLVs in the Appendix VII,
for further details on the computation, we refer the reader to Ginelli  _et
al._ (2007, 2013); Margazoglou and Magri (2023).

Two key pieces of information are contained in the CLVs. First, CLVs provide
information on the the hyperbolicity of the chaotic attractor Eckmann and
Ruelle (1985). An attractor is hyperbolic if there is a splitting of tangent
space at every point of the trajectory $\bm{x}(t)$ into three subspaces:
unstable $E^{U}_{\bm{x}}$, neutral $E^{N}_{\bm{x}}$, and stable
$E^{S}_{\bm{x}}$. These subspaces contain CLVs associated with positive, zero,
and negative LEs, respectively. This splitting has profound implications for
the geometric structure of the attractor and in hyperbolic systems when there
are no tangencies between these subspaces Huhn and Magri (2020). Therefore,
the distribution of the angles between these subspaces is bounded away from
zero, which is a property that can be analysed with the CLVs. Second, CLVs
provide information on the most important directions in phase space, which is
useful for reduced-order modelling. CLVs can split the tangent space of
spatially extended dissipative systems into a decomposition of a fixed, finite
number of “physical” modes and a remaining set of "spurious" modes Takeuchi
_et al._ (2011). The physical modes contain the relevant dynamics of the phase
space and define a finite-dimensional manifold in the tangent space, while the
spurious modes correspond to the negative LEs and increase in number with
increasing resolution. These spurious modes are hyperbolically decoupled from
the physical modes, because perturbations along them decay quickly and do not
propagate to the physical modes. A way to characterize the transition between
physical and spurious modes is by analyzing the statistics of the CLV angles
and locating a clear absence of tangencies in pairs of adjacent CLVsTakeuchi
_et al._ (2011). The number of physical modes provides a lower bound for the
degrees of freedom in the simulation.

###  IV.3 Jacobian of the LSTM

We perform stability analysis on the dynamical system defined by the LSTM. Key
to the analysis is the Jacobian of the system, therefore, the Jacobian of the
LSTM is required. To compute the Jacobian, we express in a compact form the
LSTM equations (2) during closed-loop as

$\displaystyle[\bm{c}_{i+1},\bm{h}_{i+1}]^{T}$
$\displaystyle=LSTM(\bm{\hat{x}}(t_{i}),\bm{c}_{i},\bm{h}_{i}),$ (17)
$\displaystyle\bm{\hat{x}}(t_{i+1})$
$\displaystyle=\bm{W}^{dense}\bm{h}_{i+1}+\bm{b}^{dense}.$ (18)

The Jacobian is the gradient of the internal states at a single timestep

$\displaystyle\bm{J}_{LSTM}(\bm{c}_{i},\bm{h}_{i})=\begin{bmatrix}\frac{\partial\bm{c}_{i}}{\partial\bm{c}_{i-1}}&\frac{\partial\bm{c}_{i}}{\partial\bm{h}_{i-1}}\\\\[8.61108pt]
\frac{\partial\bm{h}_{i}}{\partial\bm{c}_{i-1}}&\frac{\partial\bm{h}_{i}}{\partial\bm{h}_{i-1}}\end{bmatrix},$
(19)

which analytically is

$\displaystyle\begin{split}\frac{\partial\bm{c}_{i}}{\partial\bm{c}_{i-1}}&=\mathbf{I}*\bm{f}_{i},\\\
\frac{\partial\bm{c}_{i}}{\partial\bm{h}_{i-1}}&=\bm{c}_{i}*\bm{f}_{i}*(\mathbf{I}-\bm{f}_{i})\bm{W}^{f}+\bm{i}_{i}*(\mathbf{I}-\bm{i}_{i})\bm{W}^{i}*\tilde{\bm{c}}_{i}+\bm{i}_{i}*\left(\mathbf{I}-\bm{\tilde{C}}_{i}^{2}\right)\bm{W}^{\bm{g}},\\\
\frac{\partial\bm{h}_{i}}{\partial\bm{c}_{i-1}}&=\left(\mathbf{I}-\tanh^{2}{\left(\bm{c}_{i}\right)}\right)*\bm{o}_{i}*f_{i},\\\
\frac{\partial\bm{h}_{i}}{\partial\bm{h}_{i-1}}&=\bm{o}_{i}*(\mathbf{I}-\bm{o}_{i})*\tanh{(\bm{c}_{i})}+\left(\mathbf{I}-\tanh^{2}{(\bm{c}_{i})}\right)*\bm{o}_{i}*\\\
&\left(\bm{c}_{i}*\bm{f}_{i}*(\mathbf{I}-\bm{f}_{i})\bm{W}^{f}+\bm{i}_{i}*(\mathbf{I}-\bm{i}_{i})\bm{W}^{i}*\tilde{\bm{c}}_{i}+\bm{i}_{i}*\left(\mathbf{I}-\bm{\tilde{C}}_{i}^{2}\right)\bm{W}^{\bm{g}}\right).\end{split}$
(20)

where $\mathbf{I}\in\mathbb{R}^{N_{h}\times N_{h}}$ is the identity matrix and
$\bm{J}_{LSTM}\in\mathbb{R}^{(N_{h}+N_{h})\times(N_{h}+N_{h})}$. In chaotic
systems, the Jacobian matrix varies at each time step. To analyse the
stability from partial observations, the idea is to analyse the LSTM Jacobian
for the computation of the LEs and CLVs. With this, we can compute the
stability properties in the network in the closed-loop mode. From Eq. (20) we
expect a maximum of $2N_{h}>D$ Lyapunov exponents, of which the first $D$ are
physically relevant, whereas the remaining are spurious (i.e., they depend on
the network’s architecture). We provide a pseudocode for the calculation in
the Appendix VII.

##  V Results

The objective is to assess the capabilities of the two LSTM architectures
under investigation in the forecasting of both observed and unobserved
(hidden) chaotic timeseries. We perform a long autonomous temporal evolution
of the trained PI-LSTM and LH-LSTM networks, and compare the statistics with
the target chaotic dynamical systems. In Sec. V.1, we fix the number of
unobserved (hidden) variables in the example of the Kuramoto-Sivashinsky
equation. After reconstructing the full state, we compute the statistics and
spectrum of the kinetic energy, as well as key stability properties such as
the Lyapunov exponents, and the angles between selected CLVs. In Sec. V.2, we
employ the Lorenz 96 equation and test the capabilities of the LSTM networks
when a different number of hidden variables is considered.

###  V.1 Kuramoto–Sivashinsky

The Kuramoto–Sivashinsky equation, also known as the KS equation, was first
derived in Kuramoto (1978) to describe diffusion-induced space-time chaos in
reaction systems. In Sivashinsky (1977), the equation was independently
derived to model small diffusive-thermal instabilities of laminar flame plane
fronts. The KS equation is a fourth-order nonlinear partial differential
equation that can be written in one spatial dimension as

$\bm{\phi}_{t}(t,\mathrm{x})+\bm{\phi}_{\mathrm{x}\mathrm{x}}(t,\mathrm{x})+\bm{\phi}_{\mathrm{x}\mathrm{x}\mathrm{x}\mathrm{x}}(t,\mathrm{x})+\bm{\phi}(t,\mathrm{x})\bm{\phi}_{\mathrm{x}}(t,\mathrm{x})=0,$
(21)

with $\mathrm{x}\in[0,L)$ being the spatial direction. By assuming periodic
boundary conditions $\bm{\phi}(t,0)=\bm{\phi}(t,L)$ and setting $L=2\pi\cdot
10$, the solution of the KS equation displays chaotic behaviour Nicolaenko,
Scheurer, and Temam (1985); Hyman and Nicolaenko (1986) with
$\lambda_{1}=0.08$. We spatially discretize Eq. (21) with $128$ degrees of
freedom, corresponding to $\Delta\mathrm{x}=L/128$, and select a timestep
$\Delta t=0.25$. The equation is then solved with a fourth-order spectral
scheme for stiff PDEs Kassam and Trefethen (2005); Pathak  _et al._ (2017) up
to $T=2.5\cdot 10^{4}$, resulting in $10^{5}$ samples. To eliminate
transients, we disregard the initial $1000$ samples and partition the
remaining data into training, validation, and testing sets with sizes $5\cdot
10^{4}$, $2\cdot 10^{4}$, and $5\cdot 10^{4}$, respectively. The tangent space
calculation follows Pathak  _et al._ (2017), and the LEs and the CLVs are
computed, resulting in a Kaplan-Yorke dimension equal to $D_{KY}=15.02$. We
analyse the distribution of angles between the CLVs, and detect $j=29$
physical modes for the selected parameters in the present study; hence, the
modes $j\geq 30$ are spurious.

Figure 4: Comparison of a test data set’s (a) input, (b) label/reference
solution, (c) LH-LSTM’s closed-loop prediction, and (d) the absolute error
between the target and LH-LSTM prediction. Figure 5: Comparison of a test data
set’s (a) input, (b) label/reference solution, (c) PI-LSTM’s closed-loop
prediction, and (d) the absolute error between the target and PI-LSTM
prediction.

In the following, we present the results with $D_{y}=32$ observed inputs and
the network is trained to predict the full state with $D=128$. The batch size
is set to $128$, and we select the window size in $\\{10,20,50\\}$ and the
cell size in $\\{200,500\\}$ with a grid search. We do not observe significant
differences on the network performance for different choices of window and
cell sizes. For the physics-informed LSTM, the weighing $\alpha_{pi}$ is
selected from $\\{100,10,1\\}$.

To assess the forecast performance of the proposed LSTMs, we employ the
prediction horizon during closed-loop configuration

$\frac{\|x(t_{N_{PH}})-\hat{x}(t_{N_{PH}})\|}{\sqrt{\frac{1}{N_{PH}}\sum_{i={0}}^{N_{PH}}\|x(t_{i})\|^{2}}}<k,$
(22)

with $k=0.5$ as in Vlachas  _et al._ (2020). Starting from the same initial
conditions, a large prediction horizon means the LSTM can follow autonomously
the true chaotic evolution for long. We present the closed-loop prediction of
the two LSTMs over $6$ LTs, with the LH-LSTM in Fig. 4 and the PI-LSTM in Fig.
5. Both networks accurately predict the KS solution over a short time horizon,
with a prediction horizon of $2.6$ LT for the data-driven model and $2.86$ LT
for the physics-informed model.

We investigate the long-term physical behaviour of the network prediction by
examining the kinetic energy $E(t)=\frac{1}{2L}\int d{\rm x}|\phi(t,{\rm
x})|^{2}$ of the reference solution and the networks’ predictions in Fig. 6.
Figure 6(a) shows the kinetic energy of the first $10$ LT of the test set,
with an overlap of the energy during the prediction horizon. After the
prediction horizon, the kinetic energies gradually separate, due to the chaos,
but the LSTMs’ predictions have correct statistical behaviour in the kinetic
energy, suggesting that networks produce long-term physical solutions (Fig.
6(b)).

Figure 6: Kinetic energy of the Kuramoto–Sivashinsky system: Comparison of the
target (black line), PI-LSTM (red line) and LH-LSTM (blue line) for (a) short
time span on the test data and (b) the statistics over $500$ LT.

To investigate this further, we present in Fig. 7 the time-averaged energy
spectrum, defined as $E(k)=\frac{1}{T}\int dt|\phi(t,k)|^{2}$, which is the
distribution of energy at different wavenumbers $k$. High wavenumbers
correspond to small-scales, which are difficult to model accurately during the
numerical computation, as they involve a large range of spatial and temporal
scales. Despite passing only partial spatial observations as an input, both
LSTM models accurately capture the energy of the majority of wavenumbers,
showing that the models accurately extrapolate to smaller scales.

Figure 7: Energy spectrum of the Kuramoto–Sivashinsky system: Comparison of
the target (black line), PI-LSTM (red line) and LH-LSTM (blue line). The
distribution shows the energy transfer across different scales, with high
wavenumbers corresponding to small-scale behaviour.

Although the energy spectrum is a global quantity to analyse the distribution
of energy across different scales, stability analysis provides a detailed tool
to understand the chaotic and tangent dynamics of the system. In Fig. 8, we
compare the reference LEs (in black squares) with the LEs extracted from the
LH-LSTM (red circles) and the PI-LSTM (blue crosses). (As inVlachas  _et al._
(2020), the network does not capture two zero LEs ($\lambda_{9}$ and
$\lambda_{10}$) and the plot is augmented, accordingly.) By reconstructing the
full state the networks effectively reconstruct the tangent space, the
properties of which are encapsulated in the LEs. Both networks infer the first
$29$ LEs with high accuracy, with an error of $3.7\%$ (LH-LSTM) and $0.2\%$
(PI-LSTM) in $\lambda_{1}$. This means that the machine has inferred the
correct physical modes based on the approach in Takeuchi  _et al._ (2011). The
PI-LSTM case can also infer the hidden variables from without labelled data,
which provides further evidence of the physical accuracy and correct dynamics.

Figure 8: The first $35$ Lyapunov exponents of the Kuramoto–Sivashinsky system
for the reference data (black squares), PI-LSTM (red dots) and LH-LSTM (blue
crosses). Both networks infer correctly the first $30$ Lyapunov exponents,
which correspond to the physical behaviour (i.e., physical modes Takeuchi  _et
al._ (2011)) of the system.

We examine the principal angles between the dominant CLVs of the unstable
subspace $E_{x}^{U}$, the neutral subspace $E_{x}^{N}$, and the stable
subspace $E_{x}^{S}$ in Fig. 9, where

$\theta_{a,b}=\frac{180^{\circ}}{\pi}\cos^{-1}(|\bm{v}_{a}*\bm{v}_{a}|),$ (23)

$\theta_{a,b}\in[0^{\circ},90^{\circ}]$. In Fig. 9, the agreement of the angle
statistics between reference and the LSTMs is within a negligible numerical
error, which reflects the robust and accurate learning of the ergodic
properties from the partial input data. Further discussion on the CLVs is
provided in the Appendix VII.

Figure 9: The angle distribution of the Kuramoto–Sivashinsky system for the
leading covariant Lyapunov vectors of three different subspaces: (a) unstable-
neutral, (b) unstable-stable, and (c) neutral-stable. The black line
corresponds to the reference data, the red and blue lines indicate the results
obtained from the PI-LSTM and LH-LSTM models, respectively.

###  V.2 Lorenz-96

The Lorenz-96 model is a reduced-order model to describe the large-scale
behaviour of the mid-latitude atmosphere Lorenz (1996). It consists of a set
of coupled ordinary differential equations, which represent the variation of
an atmospheric quantity of interest such as temperature or vorticity, on a
periodic lattice representing a latitude circle on the earth

$\frac{d}{dt}x_{i}(t)=\left(x_{i+1}(t)-x_{i-2}(t)\right)x_{i-1}(t)-x_{i}(t)+F,\;\;\;\;i=1,\ldots,D,$
(24)

with
$\mathbf{x}=\begin{bmatrix}x_{1},\dots,x_{D}\end{bmatrix}\in\mathbb{R}^{D}$.
The periodic boundary conditions are $x_{1}(t)=x_{D+1}(t)$. Assuming a
constant external forcing $F=8$ and $D=20$, the system exhibits chaotic
behaviour with six positive LEsKarimi and Paul (2010). Both the numerical
solution and the reference LEs are computed with the fourth-order Runge-Kutta
method with a time step of $\Delta t=0.01$. The largest Lyapunov exponent is
$\lambda_{1}\approx 1.55$. The Kaplan-Yorke dimension is equal to
$D_{KY}=13.4$. The training set consists of $N_{t}=16000$ points, which is
equivalent to $250\tau_{\lambda}$. The hyperparameters can be found in the
Appendix VII. We deploy the PI-LSTM to reconstruct the full state in three
test cases: reconstruction of (i) $D_{\xi}=2$, (ii) $D_{\xi}=6$, (iii)
$D_{\xi}=10$ variables. We choose test cases (i) and (ii) to highlight the
capabilities of the PI-LSTM and select (iii) to demonstrate how the network
performs with significantly less measured variables for training than the
Kaplan-Yorke dimension. When the full state is available for training, both
LH-LSTM and PI-LSTM perform equally well in predicting the long-term
statistics and the LEs, but the PI-LSTM can performs well when the training
contains only partial obervations. The reconstruction is based on an
autonomous $1000\tau_{\lambda}$ long trajectory in closed-loop mode.

Figure 10: Lorenz-96: Statistics of reconstructed variables. Comparison of the
target (black line), PI-LSTM (red line) and LH-LSTM (blue line) probability
density functions (PDF) of (i) $D_{\xi}=2$, (ii) $D_{\xi}=6$, (iii)
$D_{\xi}=10$ variables over a $1000\tau_{\lambda}$ trajectory in closed-loop
configuration.

In Fig. 10, we show the statistics of the hidden variables, i.e. those
variables that are not used in the training. Despite training with the full
label, the LH-LSTM fails to predict the solution. In particular, the
corresponding delta-like or uniform distributions (in blue) indicate a fixed-
point or periodic behaviour. However, the PI-LSTM (in red) provides a markedly
more accurate statistical reconstruction of the target compared to the LH-
LSTM, with an average Wasserstein distance of (i) $0.05$, (ii) $0.15$, (iii)
$0.31$.

In Fig. 11, we compare the reference LEs (in black squares) with LEs extracted
from the data-driven LSTM (in blue circles) and PI-LSTM (in red circles) in
the three test cases. By reconstructing the variables the networks effectively
reconstruct the tangent space, the properties of which are encapsulated in the
LEs. In all cases, the LEs of the data-driven LSTM deviate significantly from
the reference LEs, differing more from the target when fewer observations are
available. In test case (i), the PI-LSTM correctly infers the LEs, with six
positive LEs. In case (ii), the PI-LSTM infers five positive LEs, decreasing
to three positive LEs for test case (iii). These findings indicate that the
Kaplan-Yorke dimension of the target chaotic system (here $D_{KY}=13.4$) plays
an important role in the reconstruction capabilities of PI-LSTM: as a
practical rule of thumb, at least $M\gtrsim D_{KY}$ independent timeseries are
required for a sufficient reconstruction of the full chaotic attractor, such
that accurate LEs statistics are extracted. This criterion is met in cases (i)
and (ii), as well as in the example of the KS equation in Sec. V.1. Still, the
PI-LSTM can sufficiently reconstruct the hidden variables statistics, as seen
in Fig. 10(iii), when the number of available timeseries is less than
$D_{KY}$.

Figure 11: Lorenz-96: Comparison of the target (black squares), PI-LSTM (red
dots) and LH-LSTM (blue dots) LEs for (i) $D_{\xi}=2$, (ii) $D_{\xi}=6$, (iii)
$D_{\xi}=10$ unobserved (hidden) variables.

####  V.2.1 Discussion on the traditional methods

As a benchmark, we reviewed a selection of methods for estimation of the
Lyapunov spectrum based on time-delay embeddingsWolf  _et al._ (1985);
Rosenstein, Collins, and De Luca (1993); Eckmann  _et al._ (1986), Extracting
LEs from time series from high dimensional chaotic attractors is a significant
numerical challenge, and it is estimatedEckmann and Ruelle (1992) that the
number of required data points, $N_{d}$, grows exponentially with the
attractor dimension $D_{KY}$ as $N_{d}\sim{\rm const}^{D_{KY}}$, where $\rm
const\sim\mathcal{O}(10)$ encapsulates the necessary parameters of those
methods. As demonstrated in Parlitz (2016) for Lorenz-96 at $D=6$ with
$D_{KY}=4.2$, additional care should be given to the chosen parameters,
requiring further experimentation. Based on these estimates, $N_{d}\geq
10^{14}$ samples are required for an adequate estimation of at least the
leading Lyapunov exponent in the case of Lorenz-96 at $D=20$ used here.
However, our approach achieves an accurate estimate of a large portion of the
Lyapunov spectrum with significantly fewer data points
($\mathcal{O}(10^{5})$).

##  VI Conclusion

We propose data-driven methods to infer the dynamics of unobserved (hidden)
chaotic variables – a task that is referred to as full-state reconstruction;
time forecast the evolution of the full state after the inference of the
hidden variables; and compute the stability properties of the reconstructed
full state. The tasks are performed from observations (data) limited to only
part of the state with long short-term memory networks (LSTMs), which are
versatile gated recurrent neural networks for sequential data (such as time
series). First, we analyse and propose architectures. The first architecture
is the low-to-high resolution LSTM (LH-LSTM), which takes partial observations
as training input, but it requires access to the full system state (full
labelled dataset). In order not to rely on the full labelled data set, we
design a physics-informed LSTM (PI-LSTM), which combines data with the
integral evolution equations. This allows for the use of existing numerical
integration methods, such as pseudospectral time-stepping schemes. The
proposed reformulation is beneficial when automatic differentiation with
respect to time is not readily available, and is therefore not limited to
recurrent neural networks. Second, we mathematically derive the Jacobian of
the LSTM. The Jacobian is key to computing the stability properties of the
chaotic attractor such as Lyapunov exponents and covariant Lyapunov vectors.
Third, we test both the LH-LSTM and the PI-LSTM on the chaotic Kuramoto-
Sivashinsky system. Our results demonstrate that both approaches correctly
perform short-term prediction and energy-spectrum inference on unseen
scenarios. The machines correctly learn the covariant Lyapunov vectors
(measured through the angles), and infer the Lyapunov spectrum of the physical
manifold. Fourth, we analyse the Lorenz-96 system. Using a purely data-driven
method with the LH-LSTM leads to markedly inaccurate long-term statistics and
unphysical Lyapunov exponents. One the other hand, the proposed PI-LSTM is
able to infer the long-term statistics and chaotic properties. Fifth, the PI-
LSTM outperforms the LH-LSTM by successfully reconstructing the hidden chaotic
dynamics when the input dimension is smaller than the Kaplan-Yorke dimension
of the attractor. This is because the missing information on the attractor is
indirectly embedded in the equations. This work opens new opportunities for
inferring hidden variables and computing the stability of dynamical systems by
combining prior knowledge of the equations and data. Current work is focused
on transferring the technology of this paper to experimental data and larger
dynamical systems.

###### Acknowledgements.

Elise Özalp thanks Alberto Racca for fruitful discussions on chaos and
recurrent neural networks. This research has received financial support from
the ERC Starting Grant No. PhyCo 949388.

## References

  * Lorenz (1963) E. N. Lorenz, “Deterministic Nonperiodic Flow.” Journal of Atmospheric Sciences 20, 130–148 (1963).
  * Strogatz (1994) S. H. Strogatz,  _Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering_ (Westview Press, 1994).
  * Huhn and Magri (2020) F. Huhn and L. Magri, “Stability, sensitivity and optimisation of chaotic acoustic oscillations,” Journal of Fluid Mechanics 882, A24 (2020).
  * Magri and Doan (2020) L. Magri and N. A. K. Doan, “Physics-informed data-driven prediction of turbulent reacting flows with lyapunov analysis and sequential data assimilation,” Data Analysis for Direct Numerical Simulations of Turbulent Combustion: From Equation-Based Analysis to Machine Learning , 177–196 (2020).
  * Takens (1981) F. Takens, “Detecting strange attractors in turbulence,” Lecture Notes in Mathematics 898, 366–381 (1981).
  * Ginelli  _et al._ (2013) F. Ginelli, H. Chaté, R. Livi, and A. Politi, “Covariant lyapunov vectors,” Journal of Physics A: Mathematical and Theoretical 46 (2013), 10.1088/1751-8113/46/25/254005.
  * Takeuchi  _et al._ (2011) K. A. Takeuchi, H.-l. Yang, F. Ginelli, G. Radons, and H. Chaté, “Hyperbolic decoupling of tangent space and effective dimension of dissipative systems,” Physical Review E 84, 046214 (2011).
  * Vlachas  _et al._ (2018) P. R. Vlachas, W. Byeon, Z. Y. Wan, T. P. Sapsis, and P. Koumoutsakos, “Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 474, 20170844 (2018).
  * Jaeger (2001) H. Jaeger, “The “echo state” approach to analysing and training recurrent neural networks,” GMD Report 148, 34–49 (2001).
  * Hochreiter and Schmidhuber (1997a) S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural computation 9, 1735–1780 (1997a).
  * Vlachas  _et al._ (2020) P. Vlachas, J. Pathak, B. Hunt, T. Sapsis, M. Girvan, E. Ott, and P. Koumoutsakos, “Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics,” Neural Networks 126, 191–217 (2020).
  * Racca and Magri (2021a) A. Racca and L. Magri, “Robust optimization and validation of echo state networks for learning chaotic dynamics,” Neural Networks 142, 252–268 (2021a).
  * Margazoglou and Magri (2023) G. Margazoglou and L. Magri, “Stability analysis of chaotic systems from data,” Nonlinear Dynamics 111, 8799–8819 (2023).
  * Lorenz (1996) E. N. Lorenz, “Predictability: a problem partly solved,” in  _Seminar on Predictability, 4-8 September 1995_, Vol. 1 (ECMWF, Shinfield Park, Reading, 1996) pp. 1–18.
  * Kuramoto (1978) Y. Kuramoto, “Diffusion-Induced Chaos in Reaction Systems,” Progress of Theoretical Physics Supplement 64, 346–367 (1978), https://academic.oup.com/ptps/article-pdf/doi/10.1143/PTPS.64.346/5293041/64-346.pdf .
  * Hyman and Nicolaenko (1986) J. M. Hyman and B. Nicolaenko, “The kuramoto-sivashinsky equation: A bridge between pde’s and dynamical systems,” Physica D: Nonlinear Phenomena 18, 113–126 (1986).
  * Hochreiter and Schmidhuber (1997b) S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural Computation 9, 1735–1780 (1997b).
  * Adrian (1991) R. J. Adrian, “Particle-imaging techniques for experimental fluid mechanics,” Annual review of fluid mechanics 23, 261–304 (1991).
  * Sangiorgio, Dercole, and Guariso (2021) M. Sangiorgio, F. Dercole, and G. Guariso, “Forecasting of noisy chaotic systems with deep neural networks,” Chaos, Solitons & Fractals 153, 111570 (2021).
  * Hornik, Stinchcombe, and White (1989) K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” Neural networks 2, 359–366 (1989).
  * Werbos (1990) P. J. Werbos, “Backpropagation through time: what it does and how to do it,” Proceedings of the IEEE 78, 1550–1560 (1990).
  * Kingma and Ba (2014) D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,”  (2014).
  * Gers, Schraudolph, and Schmidhuber (2002) F. A. Gers, N. N. Schraudolph, and J. Schmidhuber, “Learning precise timing with lstm recurrent networks,” Journal of machine learning research 3, 115–143 (2002).
  * Doan, Polifke, and Magri (2021) N. Doan, W. Polifke, and L. Magri, “Short-and long-term predictions of chaotic flows and extreme events: a physics-constrained reservoir computing approach,” Proceedings of the Royal Society A 477, 20210135 (2021).
  * Özalp, Margazoglou, and Magri (2023) E. Özalp, G. Margazoglou, and L. Magri, “Physics-informed long short-term memory for forecasting and reconstruction of chaos,” arXiv preprint arXiv:2302.10779 (2023).
  * Lagaris, Likas, and Fotiadis (1998) I. Lagaris, A. Likas, and D. Fotiadis, “Artificial neural networks for solving ordinary and partial differential equations,” IEEE Transactions on Neural Networks 9, 987–1000 (1998).
  * Raissi, Perdikaris, and Karniadakis (2019a) M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” Journal of Computational Physics 378, 686–707 (2019a).
  * Benettin  _et al._ (1980) G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, “Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. part 2: Numerical application,” Meccanica 15 (1980).
  * Boffetta  _et al._ (2002) G. Boffetta, M. Cencini, M. Falcioni, and A. Vulpiani, “Predictability: a way to characterize complexity,” Physics reports 356, 367–474 (2002).
  * Frederickson  _et al._ (1983) P. Frederickson, J. L. Kaplan, E. D. Yorke, and J. A. Yorke, “The liapunov dimension of strange attractors,” Journal of differential equations 49, 185–207 (1983).
  * Kantz and Schreiber (2003) H. Kantz and T. Schreiber,  _Nonlinear Time Series Analysis_, 2nd ed. (Cambridge University Press, 2003).
  * Oseledec (1968) V. Oseledec, “A multiplicative ergodic theorem. lyapunov characteristic numbers for dynamical systems,” Trans. Moscow Math. Soc. 19, 197–231 (1968).
  * Ruelle (1979) D. Ruelle, “Ergodic theory of differentiable dynamical systems,” Publications Mathématiques de l’IHÉS 50, 27–58 (1979).
  * Ginelli  _et al._ (2007) F. Ginelli, P. Poggi, A. Turchi, H. Chaté, R. Livi, and A. Politi, “Characterizing dynamics with covariant lyapunov vectors,” Phys. Rev. Lett. 99, 130601 (2007).
  * Eckmann and Ruelle (1985) J. P. Eckmann and D. Ruelle, “Ergodic theory of chaos and strange attractors,” Rev. Mod. Phys. 57, 617–656 (1985).
  * Sivashinsky (1977) G. Sivashinsky, “Nonlinear analysis of hydrodynamic instability in laminar flames—i. derivation of basic equations,” Acta Astronautica 4, 1177–1206 (1977).
  * Nicolaenko, Scheurer, and Temam (1985) B. Nicolaenko, B. Scheurer, and R. Temam, “Some global dynamical properties of the kuramoto-sivashinsky equations: Nonlinear stability and attractors,” Physica D: Nonlinear Phenomena 16, 155–183 (1985).
  * Kassam and Trefethen (2005) A.-K. Kassam and L. N. Trefethen, “Fourth-order time-stepping for stiff pdes,” SIAM Journal on Scientific Computing 26, 1214–1233 (2005).
  * Pathak  _et al._ (2017) J. Pathak, Z. Lu, B. R. Hunt, M. Girvan, and E. Ott, “Using machine learning to replicate chaotic attractors and calculate lyapunov exponents from data,” Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 121102 (2017), https://doi.org/10.1063/1.5010300 .
  * Karimi and Paul (2010) A. Karimi and M. R. Paul, “Extensive chaos in the lorenz-96 model,” Chaos: An interdisciplinary journal of nonlinear science 20, 043105 (2010).
  * Wolf  _et al._ (1985) A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano, “Determining lyapunov exponents from a time series,” Physica D: Nonlinear Phenomena 16, 285–317 (1985).
  * Rosenstein, Collins, and De Luca (1993) M. T. Rosenstein, J. J. Collins, and C. J. De Luca, “A practical method for calculating largest lyapunov exponents from small data sets,” Physica D: Nonlinear Phenomena 65, 117–134 (1993).
  * Eckmann  _et al._ (1986) J.-P. Eckmann, S. O. Kamphorst, D. Ruelle, and S. Ciliberto, “Liapunov exponents from time series,” Physical Review A 34, 4971 (1986).
  * Eckmann and Ruelle (1992) J.-P. Eckmann and D. Ruelle, “Fundamental limitations for estimating dimensions and lyapunov exponents in dynamical systems,” Physica D: Nonlinear Phenomena 56, 185–187 (1992).
  * Parlitz (2016) U. Parlitz, “Estimating lyapunov exponents from time series,” in  _Chaos Detection and Predictability_, edited by C. H. Skokos, G. A. Gottwald, and J. Laskar (Springer Berlin Heidelberg, Berlin, Heidelberg, 2016) pp. 1–34.
  * Sandri (1996a) M. Sandri, “Numerical calculation of lyapunov exponents,” Math. J. 6 (1996a).
  * Nastac  _et al._ (2017) G. Nastac, J. W. Labahn, L. Magri, and M. Ihme, “Lyapunov exponent as a metric for assessing the dynamic content and predictability of large-eddy simulations,” Phys. Rev. Fluids 2, 094606 (2017).
  * Raissi, Perdikaris, and Karniadakis (2019b) M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” Journal of Computational Physics 378, 686–707 (2019b).
  * Scarselli and Tsoi (1998) F. Scarselli and A. C. Tsoi, “Universal approximation using feedforward neural networks: A survey of some existing methods, and some new results,” Neural Networks 11, 15–37 (1998).
  * Libbrecht and Noble (2015) M. W. Libbrecht and W. S. Noble, “Machine learning applications in genetics and genomics,” Nature Reviews Genetics 16, 321–332 (2015).
  * Komura and Ishikawa (2018) D. Komura and S. Ishikawa, “Machine learning methods for histopathological image analysis,” Computational and Structural Biotechnology Journal 16, 34–42 (2018).
  * Lopez (2008) A. Lopez, “Statistical machine translation,” ACM Comput. Surv. 40, 8 (2008).
  * Pathak  _et al._ (2018) J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, “Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach,” Phys. Rev. Lett. 120, 024102 (2018).
  * Yu  _et al._ (2019) Y. Yu, X. Si, C. Hu, and J. Zhang, “A Review of Recurrent Neural Networks: LSTM Cells and Network Architectures,” Neural Computation 31, 1235–1270 (2019).
  * Sandri (1996b) M. Sandri, “Numerical calculation of lyapunov exponents,” Math. J. 6 (1996b).
  * Baydin  _et al._ (2018) A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, “Automatic differentiation in machine learning: a survey,” Journal of Marchine Learning Research 18, 1–43 (2018).
  * Rosenstein, Collins, and De Luca (1993) M. T. Rosenstein, J. J. Collins, and C. J. De Luca, “A practical method for calculating largest lyapunov exponents from small data sets,” Physica D: Nonlinear Phenomena 65, 117–134 (1993).
  * Schäfer and Zimmermann (2007) A. M. Schäfer and H.-G. Zimmermann, “Recurrent neural networks are universal approximators,” International journal of neural systems 17, 253–263 (2007).
  * Sherstinsky (2020) A. Sherstinsky, “Fundamentals of recurrent neural network (rnn) and long short-term memory (lstm) network,” Physica D: Nonlinear Phenomena 404, 132306 (2020).
  * Cai  _et al._ (2021) S. Cai, Z. Wang, S. Wang, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks for heat transfer problems,” Journal of Heat Transfer 143 (2021).
  * Cai  _et al._ (2022) S. Cai, Z. Mao, Z. Wang, M. Yin, and G. E. Karniadakis, “Physics-informed neural networks (pinns) for fluid mechanics: A review,” Acta Mechanica Sinica , 1–12 (2022).
  * Mao, Jagtap, and Karniadakis (2020) Z. Mao, A. D. Jagtap, and G. E. Karniadakis, “Physics-informed neural networks for high-speed flows,” Computer Methods in Applied Mechanics and Engineering 360, 112789 (2020).
  * Wang, Sankaran, and Perdikaris (2022) S. Wang, S. Sankaran, and P. Perdikaris, “Respecting causality is all you need for training physics-informed neural networks,” arXiv preprint arXiv:2203.07404 (2022).
  * Zhang and Xiao (2000) J.-S. Zhang and X.-C. Xiao, “Predicting chaotic time series using recurrent neural network,” Chinese Physics Letters 17, 88 (2000).
  * Madondo and Gibbons (2018) M. Madondo and T. Gibbons, “Learning and modeling chaos using lstm recurrent neural networks,” in  _MICS 2018 proceedings_ (2018).
  * Racca and Magri (2021b) A. Racca and L. Magri, “Automatic-differentiated physics-informed echo state network (api-esn),” in  _International Conference on Computational Science_ (Springer, 2021) pp. 323–329.
  * Che  _et al._ (2018) Z. Che, S. Purushotham, K. Cho, D. Sontag, and Y. Liu, “Recurrent neural networks for multivariate time series with missing values,” Scientific Reports 8, 6085 (2018).
  * Tian  _et al._ (2018) Y. Tian, K. Zhang, J. Li, X. Lin, and B. Yang, “Lstm-based traffic flow prediction with missing data,” Neurocomputing 318, 297–305 (2018).
  * Kelshaw, Rigas, and Magri (2022) D. Kelshaw, G. Rigas, and L. Magri, “Physics-informed cnns for super-resolution of sparse observations on dynamical systems,” (2022), arXiv:2210.17319 [physics.flu-dyn] .
  * Magri, Schmid, and Moeck (2023) L. Magri, P. J. Schmid, and J. P. Moeck, “Linear flow analysis inspired by mathematical methods from quantum mechanics,” Annual Review of Fluid Mechanics 55, 541–574 (2023).
  * Haber and Ruthotto (2017) E. Haber and L. Ruthotto, “Stable architectures for deep neural networks,” Inverse problems 34, 014004 (2017).

##  VII Appendix

Algorithm to compute the Lyapunov exponents of an LSTM

Start state $x_{t}$,

Number of time steps to compute Lyapunov times $N^{lyap}$,

Number of transient initial steps to skip for warm-up $N^{w}$,

Number of steps until QR decomposition is performed again $N^{norm}$

$U\leftarrow\textit{random}\in\mathbb{R}^{(N_{h}+N_{h})\times d}$
$\triangleright$ Initialize $d$ Gram Schmidt Vectors

$Q,R\leftarrow QR(U)$ $\triangleright$ Orthonormalize GSVs

$U\leftarrow Q$

$N^{QR}\leftarrow(N^{lyap}-N^{w})/N^{Norm}$

$\lambda\leftarrow 0\in\mathbb{R}^{d\times N^{QR}}$ Evolve the LSTM with its
cell and hidden state and GSV simultaneously for $N^{w}$ steps.

for $i=N^{w}:N^{lyap}$ do

$x(t_{i+1}),c_{i+1},h_{i+1}=LSTM(x(t_{i}),c_{i},h_{i})$ $\triangleright$ Next
LSTM step

$J\leftarrow Jac_{LSTM}(c_{i+1},h_{i+1})$ $\triangleright$ Update Jacobian

$U\leftarrow JU$ $\triangleright$ The variational equation

if $mod(i,N^{norm})=0$ then $\triangleright$ Orthonormalize every $N^{norm}$
steps

$Q,R\leftarrow QR(U)$

$U\leftarrow Q$

$\lambda[:,i/N^{norm}]\leftarrow log(diag(R))/dt$

end if

end for

Output: $\lambda_{j}\leftarrow\sum_{i=0}^{N_{QR}}\lambda[j,i]/T_{lyap}$
$\triangleright$ $j$-th Lyapunov exponent

### Computation of covariant Lyapunov vectors

Let ${\bm{M}}(t,\Delta t)=\exp\left(\int_{t}^{t+\Delta
t}{\bm{J}}(x(\tau))d\tau\right)$ (the exponential should be considered as a
path-ordered exponential Magri, Schmid, and Moeck (2023)) be the system’s
tangent evolution operator. Each bounded non-zero CLV ${\mathbf{v}}_{i}$ is
evolved by the tangent dynamics ${\bm{J}}(x(t))$ Huhn and Magri (2020)

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

where the extra term $-\lambda_{i}{\mathbf{v}}_{i}$ ensures that the norm is
bounded. The vectors ${\mathbf{v}}_{i}$ are covariant, meaning the ith CLV at
time $t_{1}$, ${\mathbf{v}}_{i}(t_{1})$, maps to the ith CLV at time $t_{2}$,
${\mathbf{v}}_{i}(t_{2})$ and vice versa. This translates to
${\bm{M}}(t,\Delta t){\mathbf{v}}_{i}(t)={\mathbf{v}}_{i}(t+\Delta t)$ and
time-invariance follows directly as ${\bm{M}}(t,-\Delta
t){\mathbf{v}}_{i}(t+\Delta t)={\mathbf{v}}_{i}(t)$. The CLVs can be recovered
from the local growth rates ${\bm{R}}(t)$ and the GSVs, ${\bm{Q}}(t)$, by
evolving backwards in time after the forward-in-time simulation is completed.
We provide a brief overview on the computation of the CLVs, for further
details on the computation, we refer the reader to Margazoglou and Magri
(2023); Ginelli  _et al._ (2013). Recall that the GSV ${\bm{Q}}(t)$ and growth
rates ${\bm{R}}(t)$ are computed by solving (1) and (14) simultaneously with
periodic orthonormalization. After $\Delta t$, the GSV evolution is given by

$\displaystyle{\bm{M}}(t,\Delta t){\bm{Q}}(t)={\bm{Q}}(t+\Delta
t){\bm{R}}(t,\Delta t)$ (26)

and the corresponding CLVs ${\bm{V}}(t)$ can be written as

$\displaystyle{\bm{V}}(t)={\bm{Q}}(t){\bm{C}}(t),$ (27)

where ${\bm{C}}(t)$ is an upper triangular matrix containing the CLV expansion
coefficients. From the CLV evolution equation

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

we replace ${\bm{V}}(t)$ using (27) and then rearrange. It then follows that
the CLV expansion coefficients can be computed by inverting the upper
triangular matrix ${\bm{R}}(t,\Delta t)$

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

after the forward-in-time computation of the GSVs is completed. This equation
is evolved backwards in time starting from the end of the forward-in-time
simulation. We employ the solve_triangular routine of scipy to invert
${\bm{R}}(t,\Delta t)$ and solve with respect to ${\bm{C}}(t)$. The diagonal
matrix ${\mathbf{D}}(t,\Delta t)$ contains the CLV local growth factors,
similar to ${\bm{R}}(t,\Delta 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. (27), to ensure that they are
converged.

### Hyperparameters

Table 1: Hyperparameters for the Kuramoto–Sivashinsky equation Hidden and cell state $N^{h}$ | $200,500$
---|---
Batch size | $128$
Physics-informed weighing $\alpha_{pi}$ | $0,10,100$
Tikhonov regularization |  $10^{-8}$, $10^{-9}$
Window size | $25$
Epochs (early stopping) | $2000$
Table 2: Hyperparameters for the Lorenz-96 model Hidden and cell state $N^{h}$ | $100,500$
---|---
Batch size | $128$
Physics-informed weighing $\alpha_{pi}$ | $0,1,10,100$
Tikhonov regularization |  $10^{-8}$, $10^{-9}$
Window size | $20,50$
Epochs (early stopping) | $2000$

### CLVs of the Kuramoto–Sivashinsky equation

In Sec. V.1, we have demonstrated that both the LH-LSTM and the PI-LSTM
accurately infer the angle statistics measured between the leading CLV for the
unstable $E^{U}_{\bm{x}}$, neutral $E^{N}_{\bm{x}}$, and stable
$E^{S}_{\bm{x}}$ subspaces of the Kuramoto-Sivashinksy equation. Specifically,
as depicted in Fig. 12(a)-(c) and (e), the angle distributions exhibit a
pronounced peak at $\pi/2$, then rapidly decline near $0$ and $\pi$. A change
of this behaviour can be found in Fig. 12(d), in which the angle distribution
spans the whole $[0,\pi]$ interval. This shows the transition between physical
and spurious modes, as discussed in Takeuchi  _et al._ (2011).

Figure 12: The angle distribution of the Kuramoto–Sivashinsky system for the
CLV pairs (a) ${\mathbf{v}}_{1},{\mathbf{v}}_{20}$, (b)
${\mathbf{v}}_{13},{\mathbf{v}}_{21}$, (c)
${\mathbf{v}}_{28},{\mathbf{v}}_{29}$, (d)
${\mathbf{v}}_{29},{\mathbf{v}}_{30}$ and (e)
${\mathbf{v}}_{30},{\mathbf{v}}_{32}$.