# Data-Driven Modeling and Forecasting of Chaotic Dynamics on Inertial
Manifolds
Constructed as Spectral Submanifolds

Aihui Liu  Institute for Mechanical Systems, ETH Zürich, Leonhardstrasse 21,
8092 Zürich, Switzerland.  Joar Axås  Institute for Mechanical Systems, ETH
Zürich, Leonhardstrasse 21, 8092 Zürich, Switzerland.  George Haller
<EMAIL_ADDRESS>Institute for Mechanical Systems, ETH Zürich,
Leonhardstrasse 21, 8092 Zürich, Switzerland.

###### Abstract

We present a data-driven and interpretable approach for reducing the
dimensionality of chaotic systems using spectral submanifolds (SSMs).
Emanating from fixed points or periodic orbits, these SSMs are low-dimensional
inertial manifolds containing the chaotic attractor of the underlying high-
dimensional system. The reduced dynamics on the SSMs turn out to predict
chaotic dynamics accurately over a few Lyapunov times and also reproduce long-
term statistical features, such as the largest Lyapunov exponents and
probability distributions, of the chaotic attractor. We illustrate this
methodology on numerical data sets including a delay-embedded Lorenz
attractor, a nine-dimensional Lorenz model, and a Duffing oscillator chain. We
also demonstrate the predictive power of our approach by constructing an SSM-
reduced model from unforced trajectories of a buckling beam, and then
predicting its periodically forced chaotic response without using data from
the forced beam.

††preprint:

> In high-dimensional systems, inertial manifolds containing the low-
> dimensional chaotic attractor play a significant role in understanding and
> predicting chaotic behaviors. Such low-dimensional, attracting inertial
> manifolds contain the essential dynamics of the system and therefore offer a
> tool for model reduction. We construct such inertial manifolds as spectral
> submanifolds (SSMs) emanating from an unstable steady state and develop SSM-
> reduced models on them. SSM-reduced models of chaotic systems give accurate
> short-term predictions and simultaneously reproduce long-term statistical
> properties. We further demonstrate the predictability of such SSM-reduced
> models on several numerical data sets including the delay-embedded Lorenz
> attractors, forced Duffing oscillator chain and a periodically forced
> buckling beam.

##  I Introduction

Chaos is ubiquitous in a variety of systems, varying from the double pendulum
to turbulent flows Strogatz (2000); Baker and Gollub (1996); Devaney (2021).
Most of these chaotic systems are high-dimensional but their chaotic dynamics
take place on low-dimensional attractors. The dynamics within these strange
attractors provides an ideal reduced model for the full system, as long as
appropriate localized coordinates can be constructed along the attractor. This
typically requires the construction of a smooth, low-dimensional manifold into
which the attractor is properly embedded.

Based on Whitney’s Whitney (1936) and Takens’s Takens (1981) embedding
theorems, several algorithms have been developed for attractor reconstruction
and for learning appropriate coordinates near those attractors from delay-
embedded data Gilpin (2020); Kraemer  _et al._ (2021); Lu, Hunt, and Ott
(2018). In general, a higher embedding dimension preserves more information at
the cost of increasing system dimensionality. To infer the optimal embedding
dimension and delay embedding time lag, using average mutual information (AMI)
Fraser and Swinney (1986) and the false nearest neighbor (FNN) method Rhodes
and Morari (1997) has become a standard practice. For high-dimensional time
series, the dimensionality of the measurements is typically higher than that
of the intrinsic system dynamics, requiring dimensional reduction for system
identification and modeling.

For this purpose, it is customary to first perform model reduction, then fit a
low-dimensional chaotic model in the reduced space. Frequently used linear
model reduction techniques, such as the principal component analysis (PCA)
Scherl  _et al._ (2020) or the dynamic mode decomposition (DMD) and its
variants Schmid (2022), cannot capture intrinsically non-linearizable
phenomena such as chaos. To counter this limitation, several nonlinear model
reduction techniques have been developed, including kernel PCA and local
linear embedding van der Maaten, Postma, and Herik (2007); Touzé, Vizzaccaro,
and Thomas (2021). Other popular methods for data-driven model reduction
include autoencoders, deep learning and neural networks. These methods aim to
extract hidden information with multiple levels of representation from
nonlinear and complex dynamical systems. Various algorithms for black-box
models have been developed, compared, and applied in various fields, including
fluid dynamics and neural science Sangiorgio, Dercole, and Guariso (2021);
Shahi, Fenton, and Cherry (2022).

To improve the physical interpretability of models produced by these
approaches, previous studies have utilized manifold learning Lee and Carlberg
(2020) or the concept of inertial manifolds Foias, Sell, and Temam (1988).
Inertial manifolds are smooth, finite-dimensional invariant manifolds that
contain the global attractor and attract all solutions exponentially. In
principle, therefore, inertial manifolds are the perfect tool for reducing
complex systems. To utilize this tool, a recent method has been developed
involving training neural networks to learn an atlas of low-dimensional
inertial manifolds in high-dimensional spatio-temporal chaotic systems Linot
and Graham (2020); Floryan and Graham (2022). These approaches show
significant promise but tend to generate intricate models.

Subsequent to model reduction, modeling chaos in the reduced coordinate is a
further challenge. Such models can be classified into local models and global
models Abarbanel  _et al._ (1993). Local models focus on capturing the
dynamics within a small region of the state space for short-term prediction.
An example of local chaos modeling is the $k$-th nearest neighbor (kNN) method
Farmer and Sidorowich (1987), which only relies on neighboring trajectory
information to make forecasts. In contrast, global models aim to capture the
overall behavior across the whole phase space. They are typically more complex
and contain more parameters. Popular global modeling techniques include neural
ODEs Linot and Graham (2022) and recurrent neural networks, such as LSTM
Vlachas  _et al._ (2018); Eivazi  _et al._ (2020). While these neural network
models have demonstrated their ability to produce forecasts with respect to
short Lyapunov times, they often require a large time domain for training
data, extensive parameter tuning, and trial-and-error model selection.
Moreover, these models are generally not trained to capture essential
statistical properties of the system, such as the largest Lyapunov exponent
Shahi, Fenton, and Cherry (2022). For chaotic neural network models, the
optimal predictor model and the model most accurately reproducing the
statistical properties of chaos are generally different Sangiorgio, Dercole,
and Guariso (2021).

A data-driven modeling alternative that offers simplicity and interpretability
is the Sparse Identification of Nonlinear Dynamical Systems (SINDy) approach,
which applies sparse regression to fit the right-hand side of a system of ODEs
to the time-history of select observables that are assumed to form a closed,
reduced dynamical system Brunton, Proctor, and Kutz (2016). While this
methodology has been successful in several settings Fasel  _et al._ (2022);
Kaptanoglu  _et al._ (2023), its results are dependent on the choice of the
selected function library, require an a priori knowledge of a reduced
observable set, and suffer from sensitivity to noise. For increased
interpretability, recent developments combine DMD or SINDy with physical
principles Baddoo  _et al._ (2023); de Silva  _et al._ (2019).

In this work, we introduce a data-driven model reduction method for chaotic
dynamics using spectral submanifolds. A spectral submanifold (SSM) is the
smoothest nonlinear continuation of a non-resonant spectral subspace of the
linear part of a system at a stationary state Haller and Ponsioen (2016). This
stationary state can be a fixed point, periodic orbit, or quasiperiodic orbit.
SSM-based model reduction has so far been used in the data-driven modeling of
non-linearizable and deterministic dynamics in mechanical systems Cenedese
_et al._ (2022a). SSM-reduction has also been successfully combined with
model-predictive control of soft robots Alora  _et al._ (2022); Mahlknecht
_et al._ (2022) and augmented with a multivariate delay-embedding technique
Axås and Haller (2023). In recent work, the mathematical theory behind SSM-
reduction has been extended to fractional and mixed-mode SSMs Haller  _et al._
(2023), with the latter type of SSM also containing unstable modes in addition
to the slowest stable modes.

This paper aims to further extend data-driven SSM-reduction to chaotic
systems. Given an observable space that has a dimension higher than the
Lyapunov dimension of the strange attractor, we seek to reconstruct a mixed-
mode SSM that contains the chaotic attractor. Such an SSM then constitutes an
explicitly constructed, very low-dimensional inertial manifold that is
anchored to a stationary state. We illustrate on examples that mixed-mode,
SSM-reduced models of chaotic systems give accurate short-term predictions
while they simultaneously reproduce long-term statistical properties, such as
the Lyapunov exponent and probability distribution density.

##  II Chaotic attractors on spectral submanifolds

###  II.1 Spectral submanifold theory

We consider a dynamical system of the form

$\dot{\bm{x}}=\bm{F}(\bm{x})=\mathcal{A}\bm{x}+\bm{f}(\bm{x}),\;\;\bm{x}\in\mathbb{R}^{n},\;\;\mathcal{A}\in\mathbb{R}^{n\times
n},\;\;\bm{f}=\mathcal{O}(|\bm{x}|^{2})\in\mathcal{C}^{\infty}.\;$ (1)

We assume that $\mathcal{A}\in\mathbb{R}^{n\times n}$ is diagonalizable and
$\bm{x}=\bm{0}$ is a hyperbolic fixed point, i.e.,

$0\notin\text{Re}\left[\text{spect}(\mathcal{A})\right],$ (2)

where
$\text{spect}(\mathcal{A})=\\{\lambda_{1},\lambda_{2},...,\lambda_{n}\\}$
denotes the spectrum of the linear part of the dynamics at the origin. We
further assume the following non-resonance conditions within
$\text{spect}(\mathcal{A})$:

$\lambda_{j}\neq\sum^{n}_{k=1}m_{k}\lambda_{k},\;\;m_{k}\in\mathbb{N},\;\;\sum^{n}_{k=1}m_{k}\geq
2,\;\;j=1,...,n.$ (3)

We then select a $d$-dimensional spectral subspace $E\subset\mathbb{R}^{n}$
that is the direct sum of a family of real eigenspaces of $\mathcal{A}$.
Without loss of generality, we assume that the spectrum of $\mathcal{A}$
within $E$, $\text{spect}(\mathcal{A}|_{E})$, contains $p$ purely real
eigenvalues and $q$ pairs of complex conjugate eigenvalues, where $p+2q=d$.
After a linear change of coordinates, the linear part of the system (1) within
the spectral subspace $E$ can be assumed diagonal. We denote the coordinates
along the real eigenvectors by $\bm{u}\in\mathbb{R}^{p}$ and their
corresponding Jordan block by

$\bm{A}=\text{diag}[\lambda_{1},...,\lambda_{p}]\in\mathbb{R}^{p\times p}.$
(4)

Similarly, we denote the coordinates along the complex conjugate $q$ pairs of
eigenvectors by $\bm{z}\in\mathbb{C}^{q}$, and their real Jordan block by

$\bm{B}_{i}=\begin{pmatrix}\alpha_{i}&\omega_{i}\\\
-\omega_{i}&\alpha_{i}\end{pmatrix}\in\mathbb{R}^{2\times 2},\;\;i=1,...,q.$
(5)

We finally assume that the rest of the spectrum of $\mathcal{A}$ outside
$\text{spect}(\mathcal{A}|_{E})$ comprises $r$ purely real eigenvalues and $s$
pairs of complex conjugate eigenvalues that correspond to the coordinates

$\bm{v}\in\mathbb{R}^{r},\;\;\bm{w}\in\mathbb{C}^{s},$ (6)

respectively, within our final set of coordinates
$(\bm{u},\bm{z},\bm{v},\bm{w})$.

Spectral submanifolds (SSMs) are invariant manifolds $\mathcal{W}(E)$ that
serve as nonlinear continuations of a $(p+2q)$-dimensional spectral subspace

$E=\\{(\bm{u},\bm{z},\bm{v},\bm{w}):\bm{v}=\bm{0},\bm{w}=\bm{0}\\}$ (7)

of the linearized dynamics at a fixed point. When the linear part of system
(1) is asymptotically stable and the non-resonance conditions (3) hold, then
there is a unique, primary $\mathcal{W}(E)$ that is as smooth as the full
system Haller and Ponsioen (2016); Cenedese  _et al._ (2022b). Recent work by
Haller  _et al._ (2023) extended this result to general hyperbolic fixed
points with general (like-mode or mixed-mode) spectral subspaces. (Like-mode
spectral subspaces contain only eigenspaces of the same stability type, in
contrast to their mixed-mode counterparts.) In this paper, we focus only on
the unique, $\mathcal{C}^{\infty}$-smooth primary SSMs, $\mathcal{W}(E)$, that
are tangent to mixed-mode spectral subspaces $E$ of hyperbolic fixed points.

By the existence theory of such primary SSMs Haller  _et al._ (2023), for any
positive integer $\mathcal{K}\geq 2$, there exists a unique set of
coefficients stored in a matrix $\bm{G}$ such that in the
$(\bm{u},\bm{z},\bm{v},\bm{w})$ coordinates introduced above, the SSM,
$\mathcal{W}(E)$, tangent to a spectral subspace $E$ can locally be written
near $\bm{x}=\bm{0}$ as a $\mathcal{C}^{\infty}$ graph over the
$(\bm{u},\bm{z})$ variables:

$\displaystyle(\bm{v},\bm{w})^{\mathrm{T}}$
$\displaystyle=\mathcal{G}(\bm{u},\bm{z})=\bm{G}\cdot(\bm{u},\bm{z},\overline{\bm{z}})^{1:\mathcal{K}}+\mathcal{O}\left(|(\bm{u},\bm{z})|^{\mathcal{K}}\right),$
(8) $\displaystyle\bm{G}$
$\displaystyle=[\bm{G}_{1},\bm{G}_{2},\dots,\bm{G}_{\mathcal{K}}]\in\mathbb{C}^{(n-d)\times\sum^{\mathcal{K}}_{k=1}d_{k}},\;\;\bm{G}_{k}\in\mathbb{C}^{(n-d)\times
d_{k}}.\;$

Here the superscript $(\cdot)^{l:r}$ denotes a vector of all monomials of the
variables listed in $(\cdot)$ from order $l$ to order $r$ and $d_{k}$ is the
total number of $d$-variate monomials composed of those variables at order
$k$. For instance, if $p=2,\;q=0$ in Eqs. (4) and (5), and the polynomial
order we select to approximate the SSM is $\mathcal{K}=3$, then we have
$(\bm{u},\bm{z},\overline{\bm{z}})=[u_{1},u_{2}]^{\mathrm{T}}$, $d_{2}=3$ and
$d_{3}=4$, generating the set of monomials

$(\bm{u},\bm{z},\overline{\bm{z}})^{2:3}=[u_{1}^{2},u_{1}u_{2},u_{2}^{2},u_{1}^{3},u_{1}^{2}u_{2},u_{1}u_{2}^{2},u_{2}^{3}]^{\mathrm{T}},$
(9)

and the coefficient matrix $\bm{G}_{2:3}\in\mathbb{R}^{(n-2)\times(3+4)}$ is
real.

The reduced dynamics on the SSM, $\mathcal{W}(E)$, are given by

$\begin{pmatrix}\dot{\bm{u}}\\\
\dot{\bm{z}}\end{pmatrix}=\begin{pmatrix}\bm{A}&0&0&0\\\ 0&\bm{B}_{1}&0&0\\\
0&0&\ddots&0\\\ 0&0&0&\bm{B}_{q}\end{pmatrix}\begin{pmatrix}\bm{u}\\\
\bm{z}\end{pmatrix}+f^{\bm{u},\bm{z}}(\bm{u},\bm{z},\mathcal{G}(\bm{u},\bm{z})),$
(10)

where $f^{\bm{u},\bm{z}}$ denotes the $(\bm{u},\bm{z})$ coordinate component
of $f=(f^{\bm{u},\bm{z}},f^{\bm{v}},f^{\bm{w}})$.

The open source Matlab software SSMTool Jain  _et al._ (2023) can compute a
Taylor series approximation of the SSM, $\mathcal{W}(E)$, up to any required
order when system (1) is explicitly known. When system (1) is only known from
data, then the open source SSMLearn Cenedese  _et al._ (2022b) algorithm can
approximate $\mathcal{W}(E)$ purely from observed trajectories.

###  II.2 SSM reconstruction via delay embedding

In practice, physical systems in which all state variables can be
simultaneously monitored are rare. This necessitates the use of invariant
manifold reconstruction methods relying only on a limited number of
observables. Specifically, using the Takens delay-embedding theorem Takens
(1981); Sauer, Yorke, and Casdagli (1991), it is possible to embed an
invariant manifold of dimension $d$ in the space of $2d+1$ (or higher) time-
shifted samples of a single, generic scalar observable. Recently, Ref. Axås
and Haller (2023) discussed the use of delay embedding to recover specifically
SSMs and their tangent spaces in such observable spaces.

To recall this methodology, we let $\mathcal{W}$ be a $d$-dimensional spectral
submanifold of (1), and denote the reduced dynamics on $\mathcal{W}$ by
$\dot{\bm{\eta}}=\bm{R}(\bm{\eta})$, $\bm{\eta}\in\mathbb{R}^{d}$, which
generates the reduced flow map
$\bm{R}^{t}:\bm{\eta_{0}}\mapsto\bm{\eta}(t,\bm{\eta_{0}})$. Then there exists
a smooth SSM parametrization $\bm{M}:\mathbb{R}^{d}\mapsto\mathbb{R}^{n}$,
such that the reduced flow map $\bm{R}^{t}$ is conjugate to the restriction of
the full flow map $\bm{F}^{t}:\mathbb{R}^{n}\mapsto\mathbb{R}^{n}$ of system
(1) to $\mathcal{W}\subset\mathbb{R}^{n}$:

$\bm{F}^{t}\circ\bm{M}=\bm{M}\circ\bm{R}^{t}.$ (11)

We then consider time-resolved, scalar-valued observations
$s(t)=\zeta(\bm{x}(t))$ of system (1) from a generic observable
$\zeta:\mathbb{R}^{n}\mapsto\mathbb{R}$. We define a delay coordinate map with
$m$ delays: $\bm{\Psi}:\mathbb{R}^{n}\mapsto\mathbb{R}^{m}$ by stacking $m$
subsequent entries of the time series $s(t)$, sampled uniformly in time with
sampling time $T_{s}$, in a vector:

$\bm{y}(t)=\bm{\Psi}(\bm{x}(t))=[s(t),s(t+T_{s}),\cdots,s(t+(m-1)T_{s})]^{\mathrm{T}}\in\mathbb{R}^{m}.$
(12)

The map $\bm{\Psi}$ is thus a mapping from the phase space of system (1) to
the delay embedded space $\mathbb{R}^{m}$. Let
$\bm{F}_{\Psi}:\mathbb{R}^{m}\mapsto\mathbb{R}^{m}$ be the flow map induced by
the flow of system (1) on this delay embedding space. By construction, if the
observation map $\zeta$ maps $\bm{0}$ to $0$, the $\bm{x}=\bm{0}$ equilibrium
point of system (1) is mapped into a fixed point $\bm{y}=\bm{0}$, i.e.,
$\bm{F}_{\Psi}(\bm{0})=\bm{0}$. According to the Takens embedding theorem
Takens (1981); Sauer, Yorke, and Casdagli (1991), for a generic observable
function $\zeta$ and under certain nondegeneracy conditions on the sampling
time $T_{s}$, if $m>2d$, then the delay coordinate map $\bm{\Psi}$ restricted
to $\mathcal{W}$ provides an embedding of $\mathcal{W}$ into $\mathbb{R}^{m}$.
As a consequence, $\tilde{\mathcal{W}}=\bm{\Psi}(\mathcal{W})$ is
diffeomorphic to the manifold $\mathcal{W}$.

The delay-embedded dynamics $\bm{F}_{\Psi}$ on $\tilde{\mathcal{W}}$ and the
reduced dynamics restricted to the invariant manifold $\mathcal{W}$ are
conjugate, i.e.,

$\bm{F}_{\Psi}\circ\bm{\Psi}=\bm{\Psi}\circ\bm{F}^{t}|_{\mathcal{W}},$ (13)

which provides the theoretical basis for identifying an SSM attached to the
$\bm{y}=\bm{0}$ origin of the $m$-dimensional delay-embedding space.

For dynamical systems with (generally non-smooth) global attractors, inertial
manifolds have been defined as finite-dimensional, forward-invariant manifolds
that are at least Lipschitz and contain the global attractor Foias, Sell, and
Temam (1988); Temam (1990). Therefore, attracting SSMs with an attractor in
their reduced dynamics function as inertial manifolds, at least locally near a
stationary state. If the attractor within the SSM is globally unique, then an
SSM is, in fact, an inertial manifold. Unlike general inertial manifolds,
however, SSMs can be precisely located and computed up to arbitrary order of
accuracy, which makes them the perfect tools for a mathematically precise
reduction of high-dimensional systems to very low-dimensional models.

In order to contain the attractor, the dimension of the SSM should be larger
than the Hausdorff dimension of the attractor Robinson (2010). In a data-
driven setting with no prior information on the attractor dimension, the
necessary embedding dimension is often estimated via the false nearest
neighbors (FNN) method Strogatz (2000); Kantz and Schreiber (2003). This
method is based on the idea that if the embedding dimension is too low,
neighboring points in the embedding space may appear closer than they actually
are, leading to false neighbors. Specifically, the FNN algorithm examines data
points in the delay-embedding space, and checks if they are indeed close to
each other in all directions. The percentage of the false nearest neighbors
goes down when the embedding dimension is increased. When the percentage
becomes lower than a predefined value, the data disentangle in the embedding
space, thus giving a correct approximation for the attractor dimension. The
SSMs we construct will have a dimension of $d$ that is at least as large as
the estimate obtained from the FNN algorithm.

##  III Learning spectral submanifolds from data

According to the SSM theory discussed in Section II.1, the primary SSM,
$\mathcal{W}(E)$, tangent to the spectral subspace $E$ at the origin, can
locally be written as a polynomial expansion over $E$. The open source Matlab
and Python implementations of SSMLearn ssm  compute the SSM parametrization
and its reduced dynamics purely from observed data Cenedese  _et al._ (2022b,
a). In this section, we discuss the computation of the SSM-reduced dynamics
and predictions from such reduced models of the full system dynamics.

###  III.1 Data-driven SSM-based dimensionality reduction

For simplicity, we only consider here examples wherein $q=0$, i.e., the chosen
spectral subspace $E$ only contains $p$ real eigenvalues. Our procedure,
however, extends directly to spectral subspaces with complex eigenvalues, such
as those arising in Haller  _et al._ (2023).

Let the coordinate $\bm{y}\in\mathbb{R}^{\rho}$ denote data obtained from
system (1) either via direct observations or delay-embedding (see Section
II.2). Also, let $\bm{\eta}\in\mathbb{R}^{d}$ denote coordinates along the
corresponding spectral subspace $E$ of $D\bm{\Psi}(E)$, respectively. As a
graph over this spectral subspace, the embedded SSM can be approximated by a
polynomial expansion. According to the $\mathcal{C}^{\infty}$ version of the
SSM theory outlined in Section II.1, for any positive integer order
$\mathcal{K}\geq 2$, the coefficients of this polynomial expansion are unique.
We, therefore, seek a parametrization
$\bm{M}:\mathbb{R}^{d}\mapsto\mathbb{R}^{\rho}$ of the embedded SSM,
$\tilde{\mathcal{W}}(E)\subset\mathbb{R}^{\rho}$, as a $\mathcal{K}$-th order
polynomial in the form

$\displaystyle\bm{M}(\bm{\eta})$
$\displaystyle=\bm{V}\bm{\eta}^{1:\mathcal{K}}=\bm{V}_{1}\bm{\eta}+\bm{V}_{2:\mathcal{K}}\bm{\eta}^{2:\mathcal{K}},$
(14) $\displaystyle\bm{V}$
$\displaystyle=[\bm{V}_{1},\bm{V}_{2},\dots,\bm{V}_{\mathcal{K}}]\in\mathbb{R}^{\rho\times\sum^{\mathcal{K}}_{k=1}d_{k}},\;\;\bm{V}_{k}\in\mathbb{R}^{\rho\times
d_{k}},$

with the reduced coordinates $\bm{\eta}=\bm{V}_{1}^{\mathrm{T}}\bm{y}$, where
the matrix $\bm{V}_{1}\in\mathbb{R}^{\rho\times d}$ has orthonormal columns
that span the tangent space of the manifold
$\tilde{\mathcal{W}}(E)\subset\mathbb{R}^{\rho}$ at $\bm{y}=\bm{0}$.

To find the constant in the above parametrization from data, we use SSMLearn
to solve the constrained minimization problem

$\displaystyle\bm{V}^{*}=[\bm{V}_{1}^{*},\bm{V}_{2:\mathcal{K}}^{*}]=$
$\displaystyle\operatorname*{arg\,min}_{\bm{V}_{1},\bm{V}_{2:\mathcal{K}}}\quad\sum_{j}||\bm{y}_{j}-\bm{V}_{1}\bm{V}_{1}^{\mathrm{T}}\bm{y}_{j}-\bm{V}_{2:\mathcal{K}}\bm{\eta}^{2:\mathcal{K}}||,$
(15)
$\displaystyle\bm{V}_{1}^{\mathrm{T}}\bm{V}_{1}=\bm{I},\quad\bm{V}_{1}^{\mathrm{T}}\bm{V}_{2:\mathcal{K}}=\bm{0},$

where the last constraint represents a basic nonlinear extension of the
principal component analysis.

To quantify the accuracy of this SSM approximation, we define the invariance
error as follows. From a specific observation $\bm{y}_{i}\neq\bm{0}$, we
compute the corresponding reduced coordinate
$\bm{\eta}_{i}=\bm{V}_{1}^{\mathrm{T}}\bm{y}_{i}$ as the projection of
$\bm{y}_{i}$ onto the tangent space $\bm{V}_{1}$. By lifting $\bm{\eta}_{i}$
to the full embedding space via the manifold parametrization (14) and
comparing it with the original data point $\bm{y}_{i}$, we obtain a
measurement of the manifold fitting accuracy for a given choice of the
approximation order $\mathcal{K}$. For a total of $P$ observations, we denote
the lifted data points by
$\bm{M}(\bm{\eta})=\bm{M}(\bm{V}_{1}^{\mathrm{T}}\bm{y}_{i})$ and define the
normalized invariance error as

$\text{Invariance
Error}=\frac{1}{||\underline{\bm{y}}||}\frac{1}{P}\sum_{i=1}^{P}||\bm{y}_{i}-\bm{M}(\bm{V}_{1}^{\mathrm{T}}\bm{y}_{i})||,\qquad||\underline{\bm{y}}||=\max_{i}||\bm{y}_{i}||.$
(16)

A larger manifold parametrization order $\mathcal{K}$ will generally decrease
the invariance error, but excessively large $\mathcal{K}$ values will lead to
overfitting. In the examples from this paper, $\mathcal{K}$ is often chosen at
the local minimum of the invariance error. The value of the invariance error
depends on specific problems, but keeping it below $0.5\%$ suffices in our
experience.

###  III.2 Prediction from SSM-reduced order models

We obtain a model from the reduced dynamics on the SSM by fitting the right-
hand side of an ODE to the numerically determined, SSM-reduced vector field in
the reduced coordinates. In this paper, we employ two alternative
approximations for the reduced dynamics models: a polynomial fit and a nearest
neighbor interpolation.

In our first prediction method, we approximate the reduced dynamics in the
form

$\dot{\bm{\eta}}=\bm{R}_{1}\bm{\eta}+\bm{R}_{2:\mathcal{K}}\bm{\eta}^{2:\mathcal{K}},$
(17)

where the coefficient matrix
$\bm{R}=[\bm{R}_{1},\bm{R}_{2:\mathcal{K}}]\in\mathbb{R}^{d\times\sum^{\mathcal{K}}_{k=1}d_{k}}$
is obtained via regression as

$\bm{R}^{*}=\operatorname*{arg\,min}_{\bm{R}}\;\sum_{j}||\dot{\bm{\eta}}_{j}-\bm{R}\bm{\eta}_{j}^{1:\mathcal{K}}||.$
(18)

We then diagonalize the linear part of the reduced dynamics (17) by a linear
transformation $\bm{\xi}=\bm{W}^{-1}\bm{\eta}$, where
$\bm{\xi}\in\mathbb{R}^{d}$ denotes modal coordinates and
$\bm{W}\in\mathbb{R}^{d\times d}$ is the matrix of the eigenvectors of
$\bm{R}_{1}$. In those new coordinates, the ODE (17) becomes

$\dot{\bm{\xi}}=\bm{\Lambda}\bm{\xi}+\bm{N}_{2:\mathcal{K}}\bm{\xi}^{2:\mathcal{K}},$
(19)

where $\bm{\Lambda}\in\mathbb{R}^{d\times d}$ is a diagonal matrix containing
the eigenvalues of $\bm{R}^{*}$, and
$\bm{N}_{2:\mathcal{K}}\in\mathbb{R}^{d\times\sum^{\mathcal{K}}_{k=2}d_{k}}$
is the coefficient matrix of nonlinear monomials of $\bm{\xi}$ from order $2$
to $\mathcal{K}$. The coefficient matrix $\bm{N}_{2:\mathcal{K}}$ can be
specified via a recursive sequence of normal form transformations that
preserve the dynamics, as explained in Cenedese  _et al._ (2022b). Here we do
not perform this step because subsequent normal form transformations are
generally defined on increasingly smaller neighborhoods of $\bm{y}=\bm{0}$.
This would limit our ability to capture larger-sized attractors on the SSM.

A global polynomial representation of the reduced dynamics might require very
high-order polynomials, and yet fail to reproduce the chaotic dynamics
Abarbanel  _et al._ (1993). As an alternative, for complicated strange
attractors, we turn to the simplest and earliest method of local forecasting
Ruppert (1996); Su (2010): we fit maps locally to data points, and then use
the resulting collection of local polynomials to form a global model.

This idea brings us to our second prediction method, the local $k$-th nearest
neighbor (kNN) approach. The kNN prediction method Farmer and Sidorowich
(1987); Yankov, Decoste, and Keogh (2006) identifies the $k$ nearest neighbors
in the training data set and then, by using various interpolation or
regression methods, makes local predictions based on those neighbors Martínez
_et al._ (2019). This method is suitable for predicting highly nonlinear or
even chaotic behavior Finkenstädt (1995). In this paper, for a prediction
based at a point $\bm{\eta}$, we first search for its $k$-th nearest neighbor
$\bm{\eta}_{j}$, $j=1,...,k$ in the training data set that minimizes the
Euclidean norm $\|\bm{\eta}-\bm{\eta}_{j}\|_{2}$. We then use a zeroth-order
linear fit to give the prediction $\bm{\eta}^{*}$ at the next time step, i.e.,
define the weight $w_{j}$ as

$w_{j}=\frac{\|\bm{\eta}-\bm{\eta}_{j}\|_{2}}{\sum^{k}_{j=1}\|\bm{\eta}-\bm{\eta}_{j}\|_{2}},$
(20)

and take the weighted sum

$\bm{\eta}^{*}=\sum^{k}_{j=1}w_{j}\|\bm{\eta}-\bm{\eta}_{j}\|_{2}.$ (21)

Although the kNN method, as other local chaos prediction methods, has the
limitation of discontinuities, large numbers of adjustable parameters, and the
ability to predict only within the training data range, as long as the
training data is sufficiently dense on the attractor, the method still has
good predictive power.

In this paper, we only apply the global polynomial fit and the kNN method, but
it is also possible to combine SSM-based model reduction with other
forecasting techniques. For instance, some global prediction methods rely on
local neighborhood information to provide the weighting functions Crutchfield
and McNamara (1987); Giona, Lentini, and Cimagalli (1991), use other function
bases such as radial basis functions Casdagli (1989), or employ machine
learning and deep learning methods Karunasinghe and Liong (2006); Sangiorgio
(2022). An advantage of SSM-based reduction over these other methods is that
it can capture statistical properties of the chaotic dynamics from relatively
limited training data, as we will see.

###  III.3 Reconstruction criteria

For non-chaotic dynamics, the prediction errors from SSM-reduced models are
quantified by the normalized mean-trajectory-error (NMTE) Cenedese  _et al._
(2022b), which is the normalized difference between the test trajectory and
the predicted trajectory. Chaotic dynamics, however, have sensitive dependence
on initial conditions, which implies that the predictions for individual
trajectories will invariably fail for longer times. Indeed, even a well-
constructed reduced model can only be expected to capture statistical measures
of the chaotic dynamics in long-term predictions. Accordingly, in this paper,
we will evaluate the ability of SSM-reduced models to reproduce both the
Lyapunov exponent and the probability density distribution of the attractors.

First, we consider the Lyapunov characteristic exponents Wolf  _et al._
(1985); Lathrop (2015), the most often used indicators of chaotic flows. Let
$\bm{x}(t;\bm{x}_{0})$ denote a trajectory of system (1) starting from the
initial condition $\bm{x}_{0}$ at $t=0$. An initial perturbation
$\bm{\varepsilon}(0)=\bm{\varepsilon}_{0}$ to $\bm{x_{0}}$ will evolve along
$\bm{x}(t;\bm{x}_{0})$ as
$\bm{\varepsilon}_{t}=\bm{\Phi}(t;\bm{x}_{0})\bm{\varepsilon}_{0}$, where
$\bm{\Phi}(t;\bm{x}_{0})$ is the fundamental matrix solution of the linearized
system (1) along the trajectory $\bm{x}(t;\bm{x}_{0})$. The maximal
exponential rate of expansion along $\bm{x}(t;\bm{x}_{0})$ is then the maximal
Lyapunov exponent

$\lambda_{\text{max}}(\bm{x}_{0})=\lim_{t\rightarrow\infty}\frac{1}{2t}\log\lambda_{\text{max}}\left[\bm{\Phi}^{\mathrm{T}}(t;\bm{x}_{0})\bm{\Phi}(t;\bm{x}_{0})\right]$
(22)

if the limit exists.

As the maximum rate of separation of nearby trajectories,
$\lambda_{\text{max}}(\bm{x}_{0})$ defines a local measure of predictability
for a dynamical system. The Lyapunov exponents are independent of the initial
condition $\bm{x}_{0}$ and depend only on the trajectory
$\bm{x}(t;\bm{x}_{0})$. If $\bm{x}(t;\bm{x}_{0})$ is contained in a chaotic
(or strange) attractor, then $\lambda_{\text{max}}(\bm{x}_{0})$ is constant
for almost all trajectories along the attractor Tél, Gruiz, and Kulacsy
(2006).

One way to approximate $\lambda_{\text{max}}(\bm{x}_{0})$ is to compute the
exponent of the rate at which trajectories from nearby initial conditions
diverge from the trajectory $\bm{x}(t;\bm{x}_{0})$, i.e., calculate

$\lambda(t;\bm{x}_{0},\bm{\epsilon}_{0})=\log\frac{|\bm{x}(t;\bm{x}_{0}+\bm{\epsilon}_{0})-\bm{x}(t;\bm{x}_{0})|}{|\bm{\epsilon}_{0}|}$
(23)

for a randomly chosen initial perturbation $\bm{\epsilon}_{0}$ to $\bm{x}_{0}$
and for large enough times $t$. For a chaotic system,
$\lambda(t;\bm{x}_{0},\bm{\epsilon}_{0})$ is approximately a linear function
of time Rosenstein, Collins, and Luca (1993). When a good linear fit to
$\lambda(t;\bm{x}_{0},\bm{\epsilon}_{0})$ is not possible, the system is
generally not chaotic Dingwell (2006). While
$\lambda(t;\bm{x}_{0},\bm{\epsilon}_{0})/t$ depends on the initial condition
$\bm{x}_{0}$ and the perturbation $\bm{\epsilon}_{0}$, this ratio will
converge to $\lambda_{\text{max}}(\bm{x}_{0})$ for generic choices of
$\bm{x}_{0}$ and $\bm{\epsilon}_{0}$ Tél, Gruiz, and Kulacsy (2006). Based on
this observation, we randomly pick two nearby initial conditions on chaotic
attractors, and use a linear fit to their temporal separation rate to estimate
the Lyapunov exponent $\lambda_{\text{max}}$ of the attractor.

A related measure of predictability is the Lyapunov time, defined as the
reciprocal of $\lambda_{\text{max}}$ of the attractor Gaspard (1998), i.e.,

$\text{Lyapunov}\;\text{Time}=\frac{1}{\lambda_{\text{max}}}.$ (24)

This quantity gives the characteristic time scale of chaotic dynamics. We call
a prediction a short-time prediction if it is carried out over times that are
low-order multiples of the Lyapunov time.

Along with the Lyapunov exponent, the probability density Ruelle (1989);
Berliner (1992) is also an invariant property of attractors and hence is often
used for comparisons between the original and modeled chaotic attractors
Dörfle and Graham (1983); Graham and Scholz (1980). If a chaotic attractor
contains the orbit $\bm{x}(t)$, its probability measure $\sigma$ is defined as
the asymptotic average of Dirac deltas along $\bm{x}(t)$ Eckmann and Ruelle
(1985), i.e.,

$\sigma=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\delta_{x(t)}dt,$ (25)

where the integral of $\delta_{x(t)}$ over the entire trajectory is equal to
one.

On a chaotic attractor, the invariant measure $\sigma$ is independent of the
initial conditions. To approximate the probability density of the chaotic
attractors, we will take a sampling time that is long enough for $\sigma$ to
converge. Along each axis of the system, we can either count the number of
data points between intervals and plot the histogram, or estimate the
probability distribution $\sigma$ based on a normal kernel function Abarbanel
_et al._ (1993). In our examples, we will compare the probability density of
the original data points with those of the predictions made by SSM-reduced
model using the same kernel function. The closeness of the two plots indicates
how statistically accurate our long-term SSM-based predictions are.

In summary, we will test SSM-based prediction quality for short times via
trajectory error calculations, and for long times via Lyapunov exponent and
probability density distribution calculations.

##  IV Applications to chaotic systems

In this section, we perform SSM-based model reduction on four numerical data
sets containing chaotic attractors.

In our first two examples, we observe only a scalar time series and delay
embed an SSM containing the attractor. In our subsequent examples, we will
assume that all phase space variables can be fully observed. In our autonomous
ODE examples, the SSM-reduction outlined in Section II.1 will be directly
applicable. In our non-autonomous (time-periodic) ODE example, we will carry
out the discrete version of the same procedure for the autonomous chaotic map.
Specifically, we take the Poincaré section and perform model reduction on the
Poincaré map, attaching the SSM to an unstable periodic orbit.

###  IV.1 Delay-embedded Lorenz attractor

As our first example, we consider the Lorenz model Lorenz (1963)

$\displaystyle\dot{x}$ $\displaystyle=\sigma(y-x),$ (26)
$\displaystyle\dot{y}$ $\displaystyle=x(\rho-z)-y,$ $\displaystyle\dot{z}$
$\displaystyle=xy-\beta z,\;$

with the classic parameter values $\sigma=10$, $\rho=28$ and $\beta=8/3$. This
system has an unstable fixed point at $(x,y,z)=\bm{0}$ and hence falls in the
setting of this paper. We will only use a single scalar signal in our analysis
and hence our methodology is fully unaware of the dimension and form of the
ODE (26). As an SSM containing a chaotic attractor must be at least three-
dimensional, Takens’s theorem guarantees success (with probability one) in
embedding the SSM in an at least seven-dimensional space. We then seek to
reconstruct the Lorenz attractor within the 3D SSM.

To compare the SSM-reduced model results with the broadly used SINDy algorithm
mentioned in the Introduction, we will adopt the same setting as in the Matlab
code given by Brunton, Proctor, and Kutz (2016). Accordingly, we will observe
only the $x$ coordinate of system (26) to generate a scalar observable time
series. We generate two trajectories (one for training and the other for
testing) over the interval $(0,100)$ with a time step of $0.001$. The test
data starts with the classic initial condition $(-8,8,27)^{\mathrm{T}}$, while
the initial condition of the training data is slightly separated
$(-8,8,27+\epsilon)^{\mathrm{T}}$, where $\epsilon=10^{-9}$. After truncating
the data over the first short period of the time interval $(0,1)$, each of
these two individual trajectories lies on the low-dimensional attractor of
system (26).

We use FNN to estimate the dimension $d$ of the SSM that prevails in the data.
The different choices of model reduction parameters and corresponding manifold
fitting results are shown in Table 1. The high FNN percentage, $31.4\%$, in
the first row of Table 1 indicates that the SSM containing the strange
attractor has a dimension higher than two. This prompts us to select the SSM
dimension as $d=3$, Based on the Takens embedding theorem, a 3D SSM can
generically be delay-embedded in a $\rho=2\times d+1=7$ dimensional space.
Table 1 shows the accuracy of SSM reconstruction in this setting under
different polynomial approximation orders for the 3D SSM.

SSM dimension ($d$) | 2 | 3 | 3 | 3
---|---|---|---|---
FNN percentage | 31.4% | 0.0% | 0.0% | 0.0%
Delay embedding dimension ($\rho$) | 5 | 7 | 7 | 7
Manifold expansion order ($\mathcal{K}$) | - | 2 | 3 | 4
Invariance Error | - | 3.85% | 1.25% | 1.45%
Table 1: Different choices of SSM dimension, false nearest neighbor percentage
(FNN), order of manifold parametrization, and the corresponding invariance
error. Due to the high FNN percentage for $d=2$, we did not compute an SSM
approximation for that case.

In delayed coordinates, the origin is a fixed point corresponding to the
origin of system (26). We therefore seek a 3D SSM attached to the origin in
the 7-dimensional embedding space and determine the reduced dynamics on that
SSM. By computing the normalized invariance error defined in Eq. (16) and
identifying its minimum, we choose the order $\mathcal{K}=3$ for the
polynomial expansion for the SSM. We find that the 3D reduced dynamics on the
SSM does admit a chaotic attractor, as shown in Fig. 1.

Figure 1: Reconstruction and prediction results of Lorenz attractor in 3D. The
inertial manifold is reconstructed by observing only $x$ coordinate as a
scalar time serie, delay embedded into 7D, and performing model reduction to
3D. Dynamics on SSM is represented by a fifth-order polynomial.

In the modal coordinates (19) that diagonalize the linear part of the SSM-
reduced dynamics, we obtain the polynomial reduced-order model on the 3D SSM
in the form

$\displaystyle\dot{u}=$
$\displaystyle-11.0194u-0.0001u^{2}-0.0005uv-0.0025uw-0.0005v^{2}+0.0026w^{2}+0.0317u^{3}$
(27)
$\displaystyle+0.0622u^{2}v-0.0011u^{2}w+0.0421uv^{2}+0.0102uvw-0.0143uw^{2}+0.0058v^{3}$
$\displaystyle+0.0235v^{2}w-0.0349vw^{2}+0.0108w^{3},$ $\displaystyle\dot{v}=$
$\displaystyle-2.3859v+0.0036uw+0.0002v^{2}+0.0006vw-0.0025w^{2}-0.0273u^{3}-0.0550u^{2}v$
$\displaystyle-0.0001u^{2}w-0.0391uv^{2}-0.0096uvw+0.0117uw^{2}-0.0059v^{3}-0.0231v^{2}w$
$\displaystyle+0.0357vw^{2}-0.0131w^{3},$ $\displaystyle\dot{w}=$
$\displaystyle+9.3871w-0.0001u^{2}-0.0004uv+0.0032uw+0.0009vw-0.0016w^{2}-0.0144u^{3}$
$\displaystyle-0.0303u^{2}v-0.0001u^{2}w-0.0233uv^{2}-0.0054uvw+0.0058uw^{2}-0.0041v^{3}$
$\displaystyle-0.0147v^{2}w+0.0246vw^{2}-0.0112w^{3}.\;$

One might wonder what we have achieved here, given that system (27) is
substantially more complicated than the original system (26) and has the same
dimension as that system. One should not forget, however, that the input to
our procedure was a single scalar observable and we used no information about
the dimension or the form of the Lorenz system (26). The procedure we follow
here is the same for any system that is only known from a scalar observable.
This procedure will generally yield a smooth manifold, the SSM, that contains
the attractor of the system.

(a)

(b)

(c)

Figure 2: Time series forecasting results for delay-embedded Lorenz attractor,
and its reconstructed Lyapunov exponent. Subplots (2(a)) and (2(b)) show more
accurate predictions from SSM-reduction than from SINDy over approximately
2.97 Lyapunov times (NMTE less than 0.1). Subplot (2(c)) plots the separation
of nearby trajectories against time. The red prediction made by the SINDy
model has the maximal Lyapunov exponent 0.23 compared with the real MLE 0.90
of system (26). A blue linear fit to the SSM prediction confirms that
trajectories separate exponentially, and the system is chaotic with the
maximal Lyapunov exponent 0.89. At the same time, the red prediction made by
SINDy suggests that the system is not chaotic.

With the SSM-reduced model, we can now make predictions for test trajectories
not used in constructing this model, and compare it with predictions obtained
from the SINDy algorithm. From the same training data, SINDy Brunton, Proctor,
and Kutz (2016) constructs a Hankel matrix by stacking 10 delayed time series
as rows, and uses singular value decomposition (SVD) to select the first three
dominant modes in these time-delay coordinates. Each column of the data is
then trained with a different sparsification parameter
($\lambda_{\text{SINDy}}=0.01,0.2,2$) to fit a third-order polynomial right-
hand side of an ODE to the data.

Figure 2 shows the comparison of both short- and long-term predictions from
SINDy and SSM reduction. Both methods are able to forecast for short times but
our approach shows a more accurate result (Fig. 2(a) and Fig. 2(b)). We then
perform a Lyapunov analysis by plotting the separation of two trajectories
against time in a log plot in Fig. 2(c). The predictions given by SSMLearn
admit a clear linear fit which indicates our model is chaotic, while the
trajectories reconstructed by SINDy are not chaotic. Indeed, as noted in
Section 4.5 in the Appendix of Ref. Brunton, Proctor, and Kutz (2016),
although the skeleton of the Lorenz attractor is captured by SINDy, the
attractor is reconstructed as a non-chaotic, quasi-periodic orbit. In
contrast, the Lyapunov exponent from the SSM-based model (27) is $0.89$, which
is very close to the real value of $0.90$ computed directly from the Lorenz
system (26).

###  IV.2 A nine-dimensional Lorenz model

We consider a 3D viscous fluid layer that is uniformly heated from below, and
let $t$ and $\bm{r}$ be the time and space variables. The density field
$\rho(\bm{r},t)$ describes the mass distribution of the fluid, and the
physical state is determined by the pressure $\bm{p}(\bm{r},t)$, the
temperature $\bm{T}(\bm{r},t)$ and the velocity field $\bm{v}(\bm{r},t)$. The
dynamic behavior is governed by three nonlinear PDEs, namely the equation of
continuity, the Navier-Stokes equation, and the equation of thermal
conductivity:

$\displaystyle\frac{\partial\bm{v}}{\partial
t}+(\bm{v}\cdot\nabla)\bm{v}=\frac{\rho}{\rho_{0}}\bm{g}-\frac{1}{\rho_{0}}\nabla\bm{p}+\nu\nabla^{2}\bm{v},\quad\mathrm{div}\;\bm{v}=0,$
(28) $\displaystyle\frac{\partial\bm{T}}{\partial
t}+(\bm{v}\cdot\nabla)\bm{T}=\chi\nabla^{2}\bm{T},\;$

where $\bm{g}=(0,0,-g)$ denotes the vector of gravity, $\rho_{0}$ the standard
density, $\nu$ the kinematic viscosity and $\chi$ the thermal conductivity.

Lorenz applied a double Fourier expansion in deriving his 3D governing
equations (26) in Ref. Lorenz (1963). Ref. Reiterer  _et al._ (1998) applied a
similar approach but extended the double Fourier expansion to triple
expansion. After truncation and selection of the leading-order-terms in the
expansion, Ref. Reiterer  _et al._ (1998) obtains the nine-dimensional ODE

$\displaystyle\dot{C}_{1}$ $\displaystyle=-\sigma
b_{1}C_{1}-C_{2}C_{4}+b_{4}C_{4}^{2}+b_{3}C_{3}C_{5}-\sigma b_{2}C_{7},$ (29)
$\displaystyle\dot{C}_{2}$ $\displaystyle=-\sigma
C_{2}+C_{1}C_{4}-C_{2}C_{5}+C_{4}C_{5}-\frac{\sigma}{2}C_{9},$
$\displaystyle\dot{C}_{3}$ $\displaystyle=-\sigma
b_{1}C_{3}+C_{2}C_{4}-b_{4}C_{2}^{2}-b_{3}C_{1}C_{5}+\sigma b_{2}C_{8},$
$\displaystyle\dot{C}_{4}$ $\displaystyle=-\sigma
C_{4}+C_{2}C_{3}-C_{2}C_{5}+C_{4}C_{5}+\frac{\sigma}{2}C_{9},$
$\displaystyle\dot{C}_{5}$ $\displaystyle=-\sigma
b_{5}C_{5}+\frac{1}{2}C_{2}^{2}-\frac{1}{2}C_{4}^{2},$
$\displaystyle\dot{C}_{6}$ $\displaystyle=-b_{6}C_{6}+C_{2}C_{9}-C_{4}C_{9},$
$\displaystyle\dot{C}_{7}$
$\displaystyle=-b_{1}C_{7}-rC_{1}+2C_{5}C_{8}-C_{4}C_{9},$
$\displaystyle\dot{C}_{8}$
$\displaystyle=-b_{1}C_{8}+rC_{3}-2C_{5}C_{7}+C_{2}C_{9},$
$\displaystyle\dot{C}_{9}$
$\displaystyle=-C_{9}-rC_{2}+rC_{4}-2C_{2}C_{6}+2C_{4}C_{6}+C_{4}C_{7}-C_{2}C_{8},\;$

where $\sigma$ is the Prandtl number, $r$ is the Rayleigh number, and the
constant parameters $b_{i}$ represent a measure of the geometry of the square
cell.

Depending on $r$ and $\sigma$, the dynamics can be stationary, periodic,
chaotic or hyperchaotic, with the latter referring to the case wherein
multiple positive Lyapunov exponents exist. For certain parameters, four
symmetric low-dimensional chaotic attractors are observed in system (29). From
Reiterer  _et al._ (1998), we choose the system parameters as

$\displaystyle\sigma$ $\displaystyle=\frac{1}{2},\;r=14.2,\;a=\frac{1}{2},$
(30) $\displaystyle b_{1}$
$\displaystyle=\frac{4(1+a^{2})}{1+2a^{2}},\;b_{2}=\frac{1+2a^{2}}{2(1+a^{2})},\;b_{3}=\frac{2(1-a^{2})}{1+a^{2}},$
$\displaystyle b_{4}$
$\displaystyle=\frac{a^{2}}{1+a^{2}},\;b_{5}=\frac{8a^{2}}{1+2a^{2}},\;b_{6}=\frac{4}{1+2a^{2}}.\;$

$\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\lambda_{4}$ | $\lambda_{5}$ | $\lambda_{6}$ | $\lambda_{7}$ | $\lambda_{8}$ | $\lambda_{9}$
---|---|---|---|---|---|---|---|---
1.9263 | -0.2741 | -0.2741 | -0.5000 | -0.6667 | -2.6667 | -3.4263 | -4.7259 | -4.7259
Table 2: Eigenvalues of the 9D Lorenz attarctor.

The origin is a fixed point of system (29) with the eigenvalues listed in
Table (2). As in our previous example, no information about these eigenvalues
will be used in our data-driven SSM reduction; we only list these eigenvalues
for reference. Note that the 1:1 resonances between $\lambda_{2}$ and
$\lambda_{3}$ and between $\lambda_{8}$ and $\lambda_{9}$ do not violate the
non-resonance condition (3). There is also a 5th-order resonance between
$\lambda_{4}$, $\lambda_{5}$ and $\lambda_{6}$, which violate Eq. (3).
However, as earlier work on SSMs shows Haller and Ponsioen (2016); Cenedese
_et al._ (2022b), a polynomial expansion for an SSM, $\mathcal{W}(E)$, can
still be computed as long as the linearized spectrum within $E$ does not
contain $\lambda_{4}$, $\lambda_{5}$ and $\lambda_{6}$. We will choose $E$ in
this fashion.

As in the 3D Lorenz example, we observe only the $C_{1}$ coordinate of system
(29) as a scalar time series and delay-embed the signal. We generate two
random initial conditions (one for training and the other for testing) within
the unit sphere of the origin in the 9D phase space and observe the first
coordinate $C_{1}$. We then select and truncate the trajectories to make them
lie close to one of the four symmetric strange attractors. Assuming no
knowledge of the original system or the dimension of the chaotic attractor, we
again use FNN to estimate the dimension $d$ of the SSM, and pick $d=3$
according to the false neighbor percentage in Table 3. By the Takens embedding
theorem, we delay-embed this 3D SSM into a $\rho=2\times d+1=7$ dimensional
space. The training and testing trajectories in this seven-dimensional delay-
embedded space are shown in Fig. 3.

Figure 3: The chaotic attractor of the 9D Lorenz model (29) in the 7D delay-embedded space, projected onto the coordinates $y_{1}$, $y_{3}$ and $y_{7}$. SSM dimension ($d$) | 2 | 3 | 3 | 3 | 3 | 3 | 3
---|---|---|---|---|---|---|---
FNN percentage | 41.6% | 0.0% | 0.0% | 0.0% | 0.0% | 0.0% | 0.0%
Delay embedding dimension ($\rho$) | 5 | 7 | 7 | 7 | 7 | 7 | 7
Manifold expansion order ($\mathcal{K}$) | - | 2 | 3 | 4 | 5 | 6 | 7
Invariance Error | - | 2.24% | 0.59% | 0.22% | 0.09% | 0.05% | 0.04 %
Table 3: Different choices of SSM dimension, false nearest neighbor percentage
(FNN), order of manifold parametrization, and the corresponding invariance
error of the delay-embedded 9D Lorenz model. Due to the high FNN percentage
for $d=2$, we did not compute an SSM approximation for that case.

Using SSMLearn, we approximate a 3D SSM attached to the origin of the 7D
delay-embedded space. We minimize the invariance error over different choices
of SSM dimension and polynomial approximation orders. As shown in Table 3,
these considerations lead us to construct a 3D SSM via a sixth-order
polynomial expansion. This construct produces an invariance error of order
$0.05\%$, which is low enough to guarantee high accuracy for our SSM-reduced
model.

In the reduced coordinates, we find that a polynomial representation of the
SSM-reduced dynamics fails to capture the complex dynamics with sufficient
accuracy. Instead, we turn to the $k$-th nearest neighbor method discussed in
Section III.2. The number of neighbors in kNN algorithm is often chosen to be
larger than the data dimension Farmer and Sidorowich (1987). Here, we select
$k=4$ and perform the kNN method on a training set containing $1.5$ million
data points. Note here that the training data for SSM fitting ($10^{5}$
points) is in general much less than the number of points required for reduced
dynamics fitting. This is because the kNN prediction algorithm requires data
to fully cover the attractor in its phase space. Fig. 4(a) shows the
prediction results from this nearest neighbor method. Since chaotic dynamics
have sensitive dependence on initial conditions, the forecast is only accurate
for 5.94 Lyapunov times.

We attempted to perform the SINDy algorithm on the same training data to
identify a seven-dimensional SINDy model and compare it to SSM-reduced model
prediction results. We tested several combinations of different model orders
(from two to six) and sparsity parameters ($\lambda_{\text{SINDy}}=0.0001$ to
$0.1$). We also tried to follow the SINDy tutorial for Kaptanoglu  _et al._
(2022) to tune the hyperparameter $\lambda_{\text{SINDy}}$ by minimizing RMSE
of predicted trajectories for different choices of $\lambda_{\text{SINDy}}$.
However, the resulting SINDy models do not contain a global attractor in
general, and their solutions tend to blow up after some time of integration.
As an example, Fig. 4 shows results from a third-order SINDy model with
$\lambda_{\text{SINDy}}=0.001$. One could also try to set different sparsity
parameters for each coordinate, as for the delay-embedded Lorenz example in
Brunton, Proctor, and Kutz (2016), but we did not pursue that exploration
here.

(a)

(b)

Figure 4: Trajectory prediction in the time domain (plot (4(a))) and in the
delay-embedding space (plot (4(b))) from the reduced models for an initial
condition not contained in the training data. The length of the SSM prediction
(blue) interval is approximately 5.94 Lyapunov times, while the SINDy model
(red) fails for trajectory integration at around 1.82 Lyapunov times. We
obtain the prediction of the SSM-reduced model by making predictions on 3D
reduced coordinates using the kNN method, then projecting the trajectory back
in 7D delay-embedded space. The training set contains $1.5$ million points.

(a)

(b)

Figure 5: Trajectory prediction error (5(a)) and the Lyapunov exponent
analysis (5(b)) on test data of the 3D SSM-reduced model. We compute the time
average of the NMTE and the nearby trajectory separation rate over 200 test
trajectories, but only plot 20 of them in this figure. The average trajectory
separation rate of system (29) is 0.032, which is closely approximated by the
mean Lyapunov exponent 0.033 computed from subplot (5(b)) of the SSM-reduced
model.

To test the statistical properties of the SSM-reduced model, we further
compute the NMTE and Lyapunov exponent of 200 pieces of the test trajectory
whose initial conditions are not contained in the training data. After adding
a small perturbation, we plot the deviation of the original and the predicted
model trajectory against time on a log scale. For 200 trajectories of the
original system (29), this analysis yields the Lyapunov exponent 0.032 as the
mean of the observed separation exponents. Remarkably, this exponent is
closely approximated by the value 0.033 yielded by our SSM-reduced model in
Fig. 5(b).

###  IV.3 Duffing oscillator chain

We now consider the classic forced-damped Duffing oscillator (Kovacic and
Brennan (2011)) given by

$\ddot{x}+\delta\dot{x}+\alpha x+\beta x^{3}=\gamma\cos{\omega t},$ (31)

where the $\delta$, $\alpha$, $\beta$ and $\gamma$ are the coefficients for
damping, linear stiffness, nonlinear stiffness, and forcing amplitude,
respectively. We use this nonlinear oscillator as the third element in an
otherwise linear chain of oscillators, as shown in Fig. 6. We denote the
displacement of the $i$-th degree of freedom by the coordinate $x_{i}$. We
then apply periodic forcing along $x_{3}$ as illustrated in Fig. 6.

Figure 6: Forced harmonic oscillator chain with a single Duffing oscillator
added. We apply periodic forcing to this Duffing oscillator.

We choose the third mass as $m_{3}=1$ kg and the remaining mass as
$m_{1}=m_{2}=m_{4}=0.1$ kg each; the third linear damping coefficient attached
to the ground is $c_{3}=0.1$ N$\cdot$s/m; the other damping coefficients are
$c_{1}=c_{2}=c_{3}=c_{4}=0.75$ N$\cdot$s/m. The linear spring stiffness is
$k=1$ N/m and the third nonlinear spring has the force relation
$F_{3}(x_{3})=k_{3}x_{3}+\beta x_{3}^{3}$, where $k_{3}=3$ N/m and
$\beta=-0.25$ N/$\text{m}^{3}$. The periodic forcing that acts on the Duffing
oscillator has an amplitude of 2.25 N and a frequency of 2 rad/s. The
resulting equations of motion is given by a system of second-order ODEs of the
form

$\displaystyle 0.1\ddot{x}_{1}+0.75\dot{x}_{1}+2x_{1}-x_{2}$
$\displaystyle=0,$ (32) $\displaystyle
0.1\ddot{x}_{2}+0.75\dot{x}_{2}-x_{1}+2x_{2}-x_{3}$ $\displaystyle=0,$
$\displaystyle\ddot{x}_{3}+0.1\dot{x}_{3}-x_{2}-x_{3}-x_{4}+0.25x_{3}^{3}$
$\displaystyle=2.25\cos{2t},$ $\displaystyle
0.1\ddot{x}_{4}+0.75\dot{x}_{4}-x_{3}+3x_{4}$ $\displaystyle=0.\;$

When written in first-order form, these equations define an eight-dimensional,
non-autonomous dynamical system. For the autonomous part of the dynamics, the
origin is an unstable fixed point. This unstable equilibrium has eight
eigenvalues, as listed in Table 4. Again, these are listed only for reference
and will not be used in our SSM construction.

$\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3,4}$ | $\lambda_{5,6}$ | $\lambda_{7,8}$
---|---|---|---|---
1.2204 | -5.8047 | $-1.5713\pm 0.6269i$ | $-3.6865\pm 4.0004i$ | $-3.7500\pm 3.9922i$
Table 4: Eigenvalues of the autonomous part of the Duffing oscillator chain.

When we add the periodic forcing term to the autonomous part of the equation,
the unstable fixed point at the origin bifurcates into an unstable periodic
orbit. As long as the amplitude of the time-dependent term in system (32) is
moderate, the corresponding $T=\pi$-periodic Poincaré map will also have an
unstable fixed point. This hyperbolic fixed point is the anchor point for our
SSM computation. Our strategy will be to locate an SSM for this map that
contains an observed chaotic attractor of the system. As in our previous
examples, this SSM will act as an inertial manifold whose restricted dynamics
provide a reduced model for the dynamics on and near the attractor. The
corresponding geometry is illustrated in Fig. 7.

Figure 7: SSM construction for the Poincaré map of system (32) in the phase
space.

(a)

(b)

Figure 8: (8(a)) Chaotic trajectories (blue) and a saddle-type periodic orbit
of system (32) projected onto the $(x_{3},\dot{x}_{3})$ plane. (8(b)) The same
trajectories and periodic orbits, as they appear as a fixed point and a
trajectory on the attractor of the Poincaré map, projected again to the
$(x_{3},\dot{x}_{3})$ plane.

For the eight-dimensional chaotic system (32), we consider all eight degrees
of freedom as observables. We select seven initial conditions (six for
training and one for testing) randomly within the unit cycle around the
origin, and generate the corresponding trajectories for $10^{5}$ times the
forcing period. We track the evolution of the origin under increasing forcing
amplitude into a hyperbolic periodic orbit using the numerical continuation
package COCO Ahsan  _et al._ (2022). In Fig. 8(a), we plot the coordinates
$x_{3}$ and $\dot{x}_{3}$ for the original continuous dynamics in blue, and
the unstable periodic orbit computed by COCO in red. When we pass to the
Poincaré map, the continuous orbits become discrete dots in the Poincaré
section as in Fig. 8(b). Based on the theory of mixed-mode SSMs Haller  _et
al._ (2023), SSMLearn computes a polynomial expansion for an SSM attached to
the hyperbolic fixed point (red dot in Fig. 8(b)). Note that our construct
relies on the existence of the hyperbolic fixed point, which we compute using
COCO from system (32). However, even when the exact location of this anchor
point is not specifically known, SSMLearn can still provide an approximation
for the SSM automatically by regressing the manifold to the available data.
The location of the fixed point itself will then follow as a prediction from
the SSM-reduced dynamics.

As in the previous examples, we use the invariance error as a criterion for
determining the correct SSM dimension. Table 5 shows different invariance
errors obtained for different choices of SSM dimension and expansion order.
Based on these results, we choose a 3D SSM construction via a 7th-order
polynomial expansion.

SSM dimension ($d$) | 2D SSM | 3D SSM
---|---|---
SSM expansion orders $(\mathcal{K})$ | 3 | 5 | 7 | 3 | 5 | 7
Invariance Error | 8.26% | 4.24% | 1.89% | 0.061% | 0.005% | 0.0001%
Table 5: Different choices of SSM dimension, manifold expansion order and the
corresponding invariance error for the oscillator chain (32). Figure 9:
Predictions for the strange attractor of the forced Duffing oscillator chain
in the 8-dimensional full state space, using the kNN prediction method on a
three-dimensional SSM. Here $k=6$ and the training set contains $6\times
10^{5}$ points. Figure 10: Probability distribution for the chaotic attractor
of the full system (32) (red) and for the attractor of the 3D, SSM-reduced
model in various sets of coordinates (blue).

After obtaining a data-driven, reduced map on the 3D SSM, we again use the kNN
algorithm discussed in Section III.2 to reconstruct the observed chaotic
attractor of the system on this SSM. We set $k=6$ and use a training set with
the total size of $6\times 10^{5}$ points. The prediction is generated by the
3D reduced model, then lifted to the eight-dimensional state space via the SSM
parametrization. In Fig. 9, we visualize the test data (red) and its
prediction (blue) by showing them in various combinations of three of the
eight coordinates. From Figure 9, we see that the red dots and the blue dots
coincide in the 8D state space, and the shape of the chaotic attractor is
accurately reconstructed. To verify the accurate reproduction of long-term
statistical features of the SSM-reduced model, we further plot the probability
distribution of the original and the reconstructed dynamics in all of the
coordinates separately. From Fig. 10, the statistical properties of the
reduced-order model are very close to those of the original attractor.

We also attempted to apply the map version of SINDy Bramburger and Kutz (2020)
to the same training eight-dimensional data (six trajectories) for maps.
However, after testing different combinations of various values of polynomial
order (from three to six) and sparsity parameter
($\lambda_{\text{SINDy}}=0.001$ to $0.05$), we failed to identify a feasible
model using SINDy and the prediction results invariably blew up after some
time. This is likely due to the high complexity of the Poincaré map dynamics.

###  IV.4 Periodically forced finite-element beam

In our final and most complex example, we use data-driven model reduction to a
mixed-mode SSM to identify a low-order model for the chaotic attractor of a
periodically forced buckling beam. The nonlinear finite-element model we use
for a von Kármán beam is a second-order approximation of the kinematics of a
beam, taking into account the deformation of the cross-section for both
bending and stretching. We specifically adopt here the finite-element Matlab
code developed in Jain, Tiso, and Haller (2018) and the same buckling beam
setting as in Haller  _et al._ (2023). The boundary condition of the beam is
set as pinned-pinned, as shown in Fig. 11, and an axial compressive force is
applied to the right node.

Figure 11: Dynamic buckling of a pinned-pinned von Kármán beam, with the
subsequent addition of a time-periodic vertical load at its midpoint.

The finite-element model of the von Kármán beam consists of 13 nodes, each
with three degrees of freedom (DOFs): axial direction, transverse direction
and rotation. Under the above boundary conditions, the beam model has 12
elements and 36 degrees of freedom. With the displacement and velocity of each
DOF counted, the full model is a 72-dimensional dynamical system in its phase
space. We adopt the physical parameters listed in Table 6 from Ref. Jain,
Tiso, and Haller (2018).

Symbol | Meaning | value | [Unit]
---|---|---|---
$t$ | time | - | [s]
$L$ | length of the beam | 2 | [m]
$h$ | height of the beam | 1e-2 | [m]
$b$ | width of the beam | 5e-2 | [m]
$A$ | width of the beam | 5e-4 | [$\mathrm{m}^{2}$]
$E$ | Young’s modulus | 190e9 | [Pa]
$\kappa$ | viscous damping rate of material | 7e6 | [Pa$\cdot$s]
$\rho$ | density | 7,850 | [kg/$\mathrm{m}^{3}$]
Table 6: Physical parameters for the von Kármán beam

Euler’s critical load for buckling is given by
$P_{n}=\frac{n^{2}\pi^{2}EI}{l^{2}}$. For the first buckling mode $n=1$, the
critical buckling force is

$F_{\mathrm{buckling}}=\frac{\pi^{2}EI}{l^{2}}=\frac{\pi^{2}Ebh^{3}}{12l^{2}}.$
(33)

For the buckling force value $F=1.1F_{\mathrm{buckling}}$, the system has
three fixed points: one unstable fixed point corresponding to purely axial
displacement, and two stable fixed points corresponding to the upper and lower
buckled states. The first 10 eigenvalues of the unstable fixed points are
listed in Table 7. Therefore, each unstable steady state has one pair of real
eigenvalues $11.06$ and $-11.10$, the other pairs correspond to stable
spirals. The positive real eigenvalue $11.06$ has the physical meaning of the
beam leaving the unstable fixed point and approaching the upper or lower
buckled states, and the negative real eigenvalue $-11.10$ corresponds to a
trajectory settling down at the unstable state. The other complex pairs of
eigenvalues represent the dynamics of higher modes. According to mixed-mode
SSM theory Haller  _et al._ (2023), there exists a two-dimensional SSM tangent
to the spectral subspaces corresponding to the two real eigenvalues. As in our
earlier examples, these eigenvalues will not be used in our SSM construction.

$\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3,4}$ | $\lambda_{5,6}$ | $\lambda_{7,8}$ | $\lambda_{9,10}$
---|---|---|---|---|---
11.06 | -11.10 | $-0.36\pm 119.36i$ | $-1.83\pm 295.56i$ | $-5.80\pm 541.50i$ | $-14.19\pm 858.19i$
Table 7: First ten eigenvalues of the buckled von Kármán beam for reference.

We generate four training trajectories from the autonomous beam model (without
the vertical forcing shown in Fig. 11), either by adding a small perturbation
to the unstable stationary state, or by applying a large transverse force
($10^{4}$ N) at the middle node. After applying this force, the beam is
further down than the lower buckled state, and will first oscillate around the
two stable fixed points with a large amplitude, then converge to one of them.
The higher modes corresponding to the complex eigenvalue pairs die out, and
the trajectories converge to a 2D mixed-mode SSM, which we locate from the
training trajectories using SSMLearn. Converged portions of two training
trajectories and the SSM are shown in Fig. 12. Based on invariance error and
NMTE for non-chaotic test data (one trajectory) in Table. 8, we construct a
7th-order SSM expansion and a 9th-order polynomial representation of the
reduced dynamics in SSMLearn.

Manifold expansion order $(\mathcal{K})$ | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---
Invariance error | 0.10% | 0.10% | 0.14% | 0.14% | 0.097%
| Reduced model order
---
(SSM expansion order $\mathcal{K}$ = 7)
7 | 8 | 9 | 10 | 11
NMTE for non-chaotic test data | 0.24% | 0.24% | 0.12% | 0.12% | 0.16%
Table 8: Different choices of SSM expansion order and their corresponding
invariance error for non-chaotic test data of the buckling beam. After
selecting the manifold expansion order to be $\mathcal{K}=7$, we choose the
SSM-model polynomial order to be 9th according to their NMTE. Figure 12: Non-
chaotic training trajectories of the buckling beam model on a 2D, autonomous
SSM.

Next, we apply transverse periodic forcing on the middle node of the beam, as
shown in Fig. 11 which makes the beam oscillate chaotically. For moderate
forcing amplitudes, the autonomous SSM, $\mathcal{W}(E)$, is known to persist
in the form of a nearby, time-periodic SSM, $\mathcal{W}(E,t)$ Haller  _et
al._ (2023). The leading-order forced dynamics on $\mathcal{W}(E,t)$ can be
obtained by simply adding to the reduced dynamics of $\mathcal{W}(E)$ the
projection of the time-periodic forcing on $E$ Haller  _et al._ (2023). This
fact enables one to make predictions for the forced response based solely on a
model trained on unforced data.

(a)

(b)

Figure 13: (13(a)) Chaotic response of the beam under periodic vertical
forcing in the full phase space. This chaotic trajectory lies on the 2D SSM
that we computed earlier from the non-chaotic data. (13(b)) Chaotic
trajectories of one testing data in two-dimensional reduced coordinates.

This procedure is implemented in SSMLearn Cenedese  _et al._ (2022b); ssm  and
will be discussed here. The periodic forcing in Fig. 11 has the amplitude
$F_{0}=3.4$ N and the frequency $\Omega=15.4$ rad/s. We then generate the
chaotic response of the buckling beam to this periodic force, and this one
chaotic trajectory as the test data for our reduced-order modeling. Figure 13
shows this test trajectory lies entirely on the 2D SSM, and can be projected
onto the two-dimensional reduced coordinates.

Specifically, up to leading order, the forced, SSM-reduced model on 2D SSM,
$\mathcal{W}(E,t)$, is in the form (Haller  _et al._ (2023))

$\dot{\bm{\eta}}=\bm{R}(\bm{\eta})+\begin{pmatrix}A\cos{\Omega t}\sin{\phi}\\\
A\cos{\Omega t}\cos{\phi}\end{pmatrix}.$ (34)

After calibrating the forcing amplitude $A$ and forcing phase $\phi$ to a
single forced trajectory (see Cenedese  _et al._ (2022b) for details), we find
that the SSM-reduced model computed from the autonomous system is indeed able
to forecast the chaotic forced response, as shown in Fig. 14. The Lyapunov
exponent of the reconstructed system $1.732$ is very close to its counterpart
$1.742$ calculated for the full, 72-dimensional system. Thus, the statistical
properties of the chaotic attractor are successfully predicted by an SSM-
reduced model in this example, too.

(a)

(b)

Figure 14: (14(a)) Prediction for a chaotic beam trajectory from an SSM-
reduced model trained on unforced, non-chaotic data, then appended with the
leading-order forcing term. The prediction lasts for about 3.31 Lyapunov time.
(14(b)) Separation of nearby trajectories against time. The original
72-dimensional system (red) and the reduced 2-dimensional model (blue) exhibit
nearly equal Lyapunov exponents (1.742 and 1.732).

Even though a short-term prediction shown in Fig. 14(a) quickly loses its
accuracy, one should not forget: This prediction by the 2D, SSM-reduced model
is made entirely based on unforced, non-chaotic trajectories and yet it does
approximate both the geometry and the long-term statistics (as described by
the maximal Lyapunov exponent) with notable accuracy for the chaotic attractor
of the full, 72-dimensional system. This shows the high physical
interpretability of SSM-reduced models, which other data-driven model
reduction techniques generally do not provide. For instance, the SINDy
algorithm we used in some of our previous examples was not realistic to run on
this high-dimensional example and the algorithm would not be equipped to make
forced predictions from unforced training data.

##  V Discussion

Using the recent theory of mixed-mode spectral submanifolds Haller  _et al._
(2023), we have developed a methodology to derive SSM-reduced models for
dynamical systems with chaotic attractors. Specifically, we combined the
SSMLearn algorithm with the method of nearest neighbors to obtain an accurate
enough representation of the SSM-reduced dynamics on large enough domains that
capture the attractors of the system. In order to contain the chaotic
attractor and the dominant behavior of the underlying system, the SSM must
have a dimension higher than the attractor dimension. In a strictly data-
driven setting, the SSM dimension can be approximated using the false nearest
neighbor method, then further determined by computing the invariance error.

This methodology is generally applicable to systems in which an inertial
manifold containing the chaotic attractor is, in fact, an SSM emanating from
an unstable steady state. In the examples treated here, that steady state was
either a fixed point or a periodic orbit. Building the inertial manifold
specifically as a spectral submanifold enables us to construct mathematically
justified, low-dimensional models for the dynamics on and near the attractors.
The models are either explicit polynomial ODEs or ODEs with right-hand sides
interpolated using nearest neighbor points on the SSM. The core part of our
approach is the open-source SSMLearn algorithm ssm , a general Matlab (Python)
package for SSM-based model reduction. We additionally deploy the false
nearest neighbor (FNN) algorithm Abarbanel and Kennel (1993) to identify the
optimal SSM dimension from data.

Our SSM-based attractor models have yielded accurate local predictions over
short Lyapunov times and also reproduced closely the statistical and ergodic
properties of chaotic attractors over longer times. Examples supporting these
conclusions included the classic Lorenz model, its nine-dimensional version, a
forced oscillator chain, and a periodically forced buckling beam. We also
compare our approach with the data-driven nonlinear modeling method SINDy
Brunton, Proctor, and Kutz (2016). While SINDy has advantages, such as faster
computation time and the ability to use less training data, it fails to
forecast complex chaotic dynamics or reconstruct the statistical properties of
attractors in our examples.

Importantly, since SSM-reduced models are fully dynamics-based, they can
predict forced chaotic responses based solely on the knowledge of unforced
(non-chaotic) trajectory data as demonstrated by our forced buckling beam
example. We are unaware of any other data-driven model reduction approach that
has this capacity.

Future applications of our methodology will include the Kuramoto–Sivashinsky
equation Cvitanović, Davidchack, and Siminos (2010) and the 2D Kolmogorov flow
Platt, Sirovich, and Fitzmaurice (1991); Chandler and Kerswell (2013), both of
which have finite-dimensional attractors and coexisting unstable steady
states. Further extensions of mixed-mode SSM theory to more general anchoring
steady states (such as invariant tori and other bounded invariant sets) are
needed for the present methodology to become applicable to inertial manifolds
that emanate from invariant sets other than fixed points or periodic orbits.

## Data availability

The code used to generate the numerical results included in this paper is
available as part of the open-source Matlab package SSMLearn, available at
https://www.dropbox.com/scl/fi/fodbixmmwb07z51g5zgtq/SSMLearn_chaotic.zip?rlkey=uua7q3bjjdpcjbkgw6rc3zmpm&dl=0.

## Acknowledgments

We acknowledge helpful discussions with Mattia Cenedese, Mingwu Li and Shobhit
Jain.

## Author Declearations

The authors declare no conflict of interest.

## References

  * Strogatz (2000) S. H. Strogatz,  _Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering_ (Westview Press, 2000).
  * Baker and Gollub (1996) G. L. Baker and J. P. Gollub,  _Chaotic Dynamics: An Introduction_, 2nd ed. (Cambridge University Press, 1996).
  * Devaney (2021) R. L. Devaney,  _An Introduction to Chaotic Dynamical Systems_ (Chapman and Hall/CRC, 2021).
  * Whitney (1936) H. Whitney, “Differentiable manifolds,” Ann. Math. 37, 645–680 (1936).
  * Takens (1981) F. Takens, “Detecting strange attractors in turbulence,” in  _Lecture Notes in Mathematics_ (Springer Berlin Heidelberg, 1981) pp. 366–381.
  * Gilpin (2020) W. Gilpin, “Deep reconstruction of strange attractors from time series,” (2020), arXiv:2002.05909 [cs.LG] .
  * Kraemer  _et al._ (2021) K. H. Kraemer, G. Datseris, J. Kurths, I. Kiss, J. Ocampo-Espindola, and N. Marwin, “A unified and automated approach to attractor reconstruction,” New J. Phys. 23, 033017 (2021).
  * Lu, Hunt, and Ott (2018) Z. Lu, B. R. Hunt, and E. Ott, “Attractor reconstruction by machine learning,” Chaos 28, 061104 (2018).
  * Fraser and Swinney (1986) A. M. Fraser and H. L. Swinney, “Independent coordinates for strange attractors from mutual information,” Phys. Rev. A 33, 1134–1140 (1986).
  * Rhodes and Morari (1997) C. Rhodes and M. Morari, “The false nearest neighbors algorithm: An overview,” Comput. Chem. Eng. 21, S1149–S1154 (1997).
  * Scherl  _et al._ (2020) I. Scherl, B. Strom, J. K. Shang, O. Williams, B. L. Polagye, and S. L. Brunton, “Robust principal component analysis for modal decomposition of corrupt fluid flows,” Phys. Rev. Fluids 5, 054401 (2020).
  * Schmid (2022) P. J. Schmid, “Dynamic mode decomposition and its variants,” Annu. Rev. Fluid Mech. 54, 225–254 (2022).
  * van der Maaten, Postma, and Herik (2007) L. van der Maaten, E. Postma, and H. Herik, “Dimensionality reduction: A comparative review,” JMLR 10 (2007).
  * Touzé, Vizzaccaro, and Thomas (2021) C. Touzé, A. Vizzaccaro, and O. Thomas, “Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques,” Nonlinear Dyn. 105, 1141–1190 (2021).
  * Sangiorgio, Dercole, and Guariso (2021) M. Sangiorgio, F. Dercole, and G. Guariso,  _Deep Learning in Multi-step Prediction of Chaotic Dynamics_ (Springer International Publishing, 2021).
  * Shahi, Fenton, and Cherry (2022) S. Shahi, F. H. Fenton, and E. M. Cherry, “Prediction of chaotic time series using recurrent neural networks and reservoir computing techniques: A comparative study,” Mach. Learn. Appl. 8, 100300 (2022).
  * Lee and Carlberg (2020) K. Lee and K. T. Carlberg, “Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders,” J. Comput. Phys. 404, 108973 (2020).
  * Foias, Sell, and Temam (1988) C. Foias, G. R. Sell, and R. Temam, “Inertial manifolds for nonlinear evolutionary equations,” J. Differ. Equations 73, 309–353 (1988).
  * Linot and Graham (2020) A. J. Linot and M. D. Graham, “Deep learning to discover and predict dynamics on an inertial manifold,” Phys. Rev. E 101, 062209 (2020).
  * Floryan and Graham (2022) D. Floryan and M. Graham, “Data-driven discovery of intrinsic dynamics,” Nat. Mach. Intell. 4 (2022), 10.1038/s42256-022-00575-4.
  * Abarbanel  _et al._ (1993) H. D. I. Abarbanel, R. Brown, J. J. Sidorowich, and L. S. Tsimring, “The analysis of observed chaotic data in physical systems,” Rev. Mod. Phys. 65, 1331–1392 (1993).
  * Farmer and Sidorowich (1987) J. D. Farmer and J. J. Sidorowich, “Predicting chaotic time series,” Phys. Rev. Lett. 59, 845–848 (1987).
  * Linot and Graham (2022) A. J. Linot and M. D. Graham, “Data-driven reduced-order modeling of spatiotemporal chaos with neural ordinary differential equations,” Chaos 32 (2022), 10.1063/5.0069536, 073110.
  * 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,” Proc. R. Soc. A: Math. Phys. Eng. Sci. 474, 20170844 (2018).
  * Eivazi  _et al._ (2020) H. Eivazi, H. Veisi, M. H. Naderi, and V. Esfahanian, “Deep neural networks for nonlinear model order reduction of unsteady flows,” Phys. Fluids 32, 105104 (2020).
  * Brunton, Proctor, and Kutz (2016) S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” PNAS 113, 3932–3937 (2016).
  * Fasel  _et al._ (2022) U. Fasel, J. N. Kutz, B. W. Brunton, and S. L. Brunton, “Ensemble-sindy: Robust sparse model discovery in the low-data, high-noise limit, with active learning and control,” Proc. R. Soc. A: Math. Phys. Eng. Sci. 478, 20210904 (2022).
  * Kaptanoglu  _et al._ (2023) A. A. Kaptanoglu, L. Zhang, Z. G. Nicolaou, U. Fasel, and S. L. Brunton, “Benchmarking sparse system identification with low-dimensional chaos,” Nonlinear Dyn. (2023), 10.1007/s11071-023-08525-4.
  * Baddoo  _et al._ (2023) P. J. Baddoo, B. Herrmann, B. J. McKeon, J. Nathan Kutz, and S. L. Brunton, “Physics-informed dynamic mode decomposition,” Proc. R. Soc. A: Math. Phys. Eng. Sci. 479, 20220576 (2023).
  * de Silva  _et al._ (2019) B. de Silva, D. M. Higdon, S. L. Brunton, and J. N. Kutz, “Discovery of physics from data: Universal laws and discrepancies,” Front. Artif. Intell. 3 (2019).
  * Haller and Ponsioen (2016) G. Haller and S. Ponsioen, “Nonlinear normal modes and spectral submanifolds: existence, uniqueness and use in model reduction,” Nonlinear Dyn. 86, 1493–1534 (2016).
  * Cenedese  _et al._ (2022a) M. Cenedese, J. Axås, H. Yang, M. Eriten, and G. Haller, “Data-driven nonlinear model reduction to spectral submanifolds in mechanical systems,” Philos. Trans. Royal Soc. A PHILOS T R SOC A 380, 20210194 (2022a).
  * Alora  _et al._ (2022) J. I. Alora, M. Cenedese, E. Schmerling, G. Haller, and M. Pavone, “Data-driven spectral submanifold reduction for nonlinear optimal control of high-dimensional robots,” 2023 IEEE Int. Conf. Robot. Autom. , 2627–2633 (2022).
  * Mahlknecht  _et al._ (2022) F. Mahlknecht, J. I. Alora, S. Jain, E. Schmerling, R. Bonalli, G. Haller, and M. Pavone, “Using spectral submanifolds for nonlinear periodic control,” in  _2022 IEEE 61st Conference on Decision and Control (CDC)_ (2022) pp. 6548–6555.
  * Axås and Haller (2023) J. Axås and G. Haller, “Model reduction for nonlinearizable dynamics via delay-embedded spectral submanifolds,” (2023), arXiv:2301.04560 [math.DS] .
  * Haller  _et al._ (2023) G. Haller, B. Kaszás, A. Liu, and J. Axås, “Nonlinear model reduction to fractional and mixed-mode spectral submanifolds,” Chaos 33 (2023), 10.1063/5.0143936.
  * Cenedese  _et al._ (2022b) M. Cenedese, J. Axås, B. Bäuerlein, K. Avila, and G. Haller, “Data-driven modeling and prediction of non-linearizable dynamics via spectral submanifolds,” Nat. Commun. 13 (2022b).
  * Jain  _et al._ (2023) S. Jain, T. Thurnher, Mingwu Li, and G. Haller, “SSMTool: Computation of invariant manifolds in high-dimensional mechanics problems,”  (2023).
  * Sauer, Yorke, and Casdagli (1991) T. Sauer, J. A. Yorke, and M. Casdagli, “Embedology,” J. Stat. Phys. 65, 579–616 (1991).
  * Temam (1990) R. Temam, “Inertial manifolds,” Math. Intell. 12, 68–74 (1990).
  * Robinson (2010) J. C. Robinson,  _Dimensions, Embeddings, and Attractors_, Cambridge Tracts in Mathematics (Cambridge University Press, 2010).
  * Kantz and Schreiber (2003) H. Kantz and T. Schreiber,  _Nonlinear Time Series Analysis_, 2nd ed. (Cambridge University Press, 2003).
  * (43) “SSMLearn,” https://www.dropbox.com/scl/fi/fodbixmmwb07z51g5zgtq/SSMLearn_chaotic.zip?rlkey=uua7q3bjjdpcjbkgw6rc3zmpm&dl=0.
  * Ruppert (1996) D. Ruppert, “Local polynomial regression and its applications in environmental statistics,” Statist. Environ. 3 (1996).
  * Su (2010) L.-y. Su, “Prediction of multivariate chaotic time series with local polynomial fitting,” Comput. Math. Appl. 59, 737–744 (2010).
  * Yankov, Decoste, and Keogh (2006) D. Yankov, D. Decoste, and E. Keogh,  _Ensembles of Nearest Neighbor Forecasts_, Vol. 4212 (2006) pp. 545–556.
  * Martínez  _et al._ (2019) F. Martínez, M. Frías, M. Pérez, and A. Rivera Rivas, “A methodology for applying k-nearest neighbor to time series forecasting,” Artif. Intell. Rev. 52 (2019), 10.1007/s10462-017-9593-z.
  * Finkenstädt (1995) B. Finkenstädt, “A nearest neighbor approach to forecast nonlinear time series,” in  _Nonlinear Dynamics in Economics: A Theoretical and Statistical Approach to Agricultural Markets_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 1995) pp. 122–147.
  * Crutchfield and McNamara (1987) J. P. Crutchfield and B. S. McNamara, “Equations of motion from a data series,” Complex Syst. 1 (1987).
  * Giona, Lentini, and Cimagalli (1991) M. Giona, F. Lentini, and V. Cimagalli, “Functional reconstruction and local prediction of chaotic time series,” Phys. Rev. A 44, 3496–3502 (1991).
  * Casdagli (1989) M. Casdagli, “Nonlinear prediction of chaotic time series,” Phys. D: Nonlinear Phenom. 35, 335–356 (1989).
  * Karunasinghe and Liong (2006) D. S. Karunasinghe and S.-Y. Liong, “Chaotic time series prediction with a global model: Artificial neural network,” J. Hydrol. 323, 92–105 (2006).
  * Sangiorgio (2022) M. Sangiorgio, “Deep learning in multi-step forecasting of chaotic dynamics,” in  _Special Topics in Information Technology_, edited by L. Piroddi (Springer International Publishing, Cham, 2022) pp. 3–14.
  * Wolf  _et al._ (1985) A. Wolf, J. Swift, H. Swinney, and J. Vastano, “Determining Lyapunov exponents from a time series,” Phys. D: Nonlinear Phenom. 16, 285–317 (1985).
  * Lathrop (2015) D. Lathrop, “Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering,” Phys. Today 68, 54–55 (2015).
  * Tél, Gruiz, and Kulacsy (2006) T. Tél, M. Gruiz, and K. Kulacsy,  _Chaotic dynamics: An introduction based on classical mechanics_ (2006) pp. 1–408.
  * Rosenstein, Collins, and Luca (1993) M. T. Rosenstein, J. J. Collins, and C. J. D. Luca, “A practical method for calculating largest lyapunov exponents from small data sets,” Phys. D: Nonlinear Phenom. 65, 117–134 (1993).
  * Dingwell (2006) J. B. Dingwell, “Lyapunov exponents,” in  _Wiley Encyclopedia of Biomedical Engineering_ (John Wiley & Sons, Ltd, 2006).
  * Gaspard (1998) P. Gaspard,  _Chaos, Scattering and Statistical Mechanics_, Cambridge Nonlinear Science Series (Cambridge University Press, 1998).
  * Ruelle (1989) D. Ruelle,  _Chaotic Evolution and Strange Attractors_, Lezioni Lincee (Cambridge University Press, 1989).
  * Berliner (1992) L. M. Berliner, “Statistics, Probability and Chaos,” Stat. Sci. 7, 69 – 90 (1992).
  * Dörfle and Graham (1983) M. Dörfle and R. Graham, “Probability density of the lorenz model,” Phys. Rev. A 27, 1096–1105 (1983).
  * Graham and Scholz (1980) R. Graham and H. J. Scholz, “Analytic approximation of the lorenz attractor by invariant manifolds,” Phys. Rev. A 22, 1198–1204 (1980).
  * Eckmann and Ruelle (1985) J.-P. Eckmann and D. Ruelle, “Ergodic theory of chaos and strange attractors,” Rev. Mod. Phys. 57, 617–656 (1985).
  * Lorenz (1963) E. N. Lorenz, “Deterministic Nonperiodic Flow.” J. Atmos. Sci. 20, 130–148 (1963).
  * Reiterer  _et al._ (1998) P. Reiterer, C. Lainscsek, F. Schuerrer, C. Letellier, and J. Maquet, “A nine-dimensional Lorenz system to study high-dimensional chaos,” J. Phys. A 31, 7121–7139 (1998).
  * Kaptanoglu  _et al._ (2022) A. Kaptanoglu, B. de Silva, U. Fasel, K. Kaheman, A. Goldschmidt, J. Callaham, C. Delahunt, Z. Nicolaou, K. Champion, J.-C. Loiseau, J. Kutz, and S. Brunton, “PySINDy: A comprehensive python package for robust sparse system identification,” J. Open Source Softw. 7, 3994 (2022).
  * Kovacic and Brennan (2011) I. Kovacic and M. J. Brennan, “Background: On georg duffing and the duffing equation,” in  _The Duffing Equation_ (John Wiley & Sons, Ltd, 2011) Chap. 1, pp. 1–23.
  * Ahsan  _et al._ (2022) Z. Ahsan, H. Dankowicz, M. Li, and J. Sieber, “Methods of continuation and their implementation in the COCO software platform with application to delay differential equations,” Nonlinear Dyn. 107, 3181–3243 (2022).
  * Bramburger and Kutz (2020) J. J. Bramburger and J. N. Kutz, “Poincaré maps for multiscale physics discovery and nonlinear floquet theory,” Phys. D: Nonlinear Phenom. 408, 132479 (2020).
  * Jain, Tiso, and Haller (2018) S. Jain, P. Tiso, and G. Haller, “Exact nonlinear model reduction for a von Kármán beam: Slow-fast decomposition and spectral submanifolds,” J. Sound Vib. 423, 195–211 (2018).
  * Abarbanel and Kennel (1993) H. D. I. Abarbanel and M. B. Kennel, “Local false nearest neighbors and dynamical dimensions from observed chaotic data,” Phys. Rev. E 47, 3057–3068 (1993).
  * Cvitanović, Davidchack, and Siminos (2010) P. Cvitanović, R. L. Davidchack, and E. Siminos, “On the state space geometry of the kuramoto–sivashinsky flow in a periodic domain,” SIAM J. Appl. Dyn. Syst. 9, 1–33 (2010).
  * Platt, Sirovich, and Fitzmaurice (1991) N. Platt, L. Sirovich, and N. Fitzmaurice, “An investigation of chaotic kolmogorov flows,” Phys. Fluids A 3, 681–696 (1991).
  * Chandler and Kerswell (2013) G. J. Chandler and R. R. Kerswell, “Invariant recurrent solutions embedded in a turbulent two-dimensional kolmogorov flow,” J. Fluid Mech. 722, 554–595 (2013).