# Discovering Governing Equations from Partial Measurements with Deep Delay
Autoencoders

Joseph Bakarji1∗, Kathleen Champion2†, J. Nathan Kutz2, Steven L. Brunton1
1 Department of Mechanical Engineering, University of Washington, Seattle, WA
98195, United States
2 Department of Applied Mathematics, University of Washington, Seattle, WA
98195, United States

###### Abstract

A central challenge in data-driven model discovery is the presence of hidden,
or latent, variables that are not directly measured but are dynamically
important. Takens’ theorem provides conditions for when it is possible to
augment these partial measurements with time delayed information, resulting in
an attractor that is diffeomorphic to that of the original full-state system.
However, the coordinate transformation back to the original attractor is
typically unknown, and learning the dynamics in the embedding space has
remained an open challenge for decades. Here, we design a custom deep
autoencoder network to learn a coordinate transformation from the delay
embedded space into a new space where it is possible to represent the dynamics
in a sparse, closed form. We demonstrate this approach on the Lorenz, Rössler,
and Lotka-Volterra systems, learning dynamics from a single measurement
variable. As a challenging example, we learn a Lorenz analogue from a single
scalar variable extracted from a video of a chaotic waterwheel experiment. The
resulting modeling framework combines deep learning to uncover effective
coordinates and the sparse identification of nonlinear dynamics (SINDy) for
interpretable modeling. Thus, we show that it is possible to simultaneously
learn a closed-form model and the associated coordinate system for partially
observed dynamics.

_Keywords–_ Model discovery, machine learning, deep learning, dynamical
systems, time delays, Takens embedding, sparse nonlinear modeling,
autoencoder, attractor reconstruction

††∗ Corresponding author (jbakarji@uw.edu).††† Work performed independently of
employment at Amazon.

##  1 Introduction

Scientific progress has been driven by the discovery of simple and predictive
mathematical models from observations. Interpretable, parsimonious governing
equations have been especially valuable as they typically have allowed for
greater engineering insight, simple parametrizations, and improved
extrapolation capabilities. Increasingly, these models are learned from data
using machine learning algorithms, such as genetic programming [1, 2, 3, 4]
and the sparse identification of nonlinear dynamics (SINDy) [5, 6]. When only
partial observations of a dynamical system are available, so that some states
remain hidden, it is generally not possible to formulate a closed-form model
in these variables, using either analytic or data-driven techniques. Time-
delay embedding provides an approach to augment these limited measurements,
and under certain conditions, given by Takens’ embedding theorem [7], the
delay-augmented state yields an attractor that is diffeomorphic (i.e.
continuously differentiably deformable) to the underlying, though unmeasured,
full-state attractor. However, the coordinate transformation back to the
original attractor is unknown and may be arbitrarily complex to represent.
Learning interpretable and generalizable dynamical systems models in the
embedding space has remained an open challenge since Takens introduced the
embedding theorem in 1981 [7], with several promising recent developments [4,
8, 9, 10]. The central goal of this work is to learn parsimonious, closed-form
models from partially observed systems, leveraging the SINDy algorithm for
model discovery and a custom deep autoencoder to learn an appropriate
coordinate transformation from a time-delay embedding.

Time-delay embedding has a rich history in data-driven modeling [11, 12, 13,
14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30],
especially for making predictions when analytical models are not available.
Notable applications of the theory include distinguishing noise from chaotic
dynamics [11, 15], and recent work has shown that delay embedding can improve
machine learning models of time-series data [31, 32]. Numerous attempts have
been made to learn dynamical systems models directly in the delay embedded
space. For example, when SINDy was applied to a time-delay embedding for
chaotic systems, it was shown that the simplest model that describes the data
is linear, establishing a strong connection to Koopman operator theory [22].
This connection between delay embedding and Koopman has a longer history [33,
34, 22, 35, 36, 25, 26, 29, 37, 38], including the use of dynamic mode
decomposition (DMD) [39, 40] on time-delay coordinates [41, 42]. However, for
chaotic systems, these linear models in delay coordinates suffer from closure
issues and fall short of the goal of parsimonious nonlinear models. In related
work, attempts have been made to extend SINDy to partially measured systems,
using the variational annealing data assimilation technique [10] or applying
SINDy to higher derivatives of the measured variables [8], which is closely
related to time-delay embedding.

Figure 1: Summary of workflow. ${\bm{\phi}}$, ${\bm{\psi}}$ and $\mathbf{f}$
are parametric functions to be discovered.

Because the data in the delay-embedded space is diffeomorphic to the full-
state attractor, where a parsimonious nonlinear model often exists, it is
reasonable to try and learn a coordinate transformation starting from a delay-
embedded state. Figure 1 shows a schematic of the modeling task, involving
learning a coordinate transform into a space where an analytical model may be
obtained. Deep neural network autoencoders have shown great promise in
representing arbitrarily complex coordinate transformations, especially for
dynamical systems [43, 44, 45, 46, 32, 47, 48, 49, 50, 51, 52]. Importantly,
Gilpin [50] introduced a false nearest neighbors loss function in a deep
neural network to promote embedded attractors that are similar to the full-
state attractor. Similarly, Wehmeyer and Noé [32] showed that autoencoders may
be combined with time-delayed measurements to learn effective slow reaction
coordinates for molecular kinetics. However, neither of these studies learned
parsimonious models in the autoencoder coordinates. Champion et al. [47]
showed that it is possible to combine autoencoder networks with SINDy to
simultaneously learn a nonlinear coordinate transformation and discover a
parsimonious dynamics model, making this a promising technique for learning
models in delay-embedded coordinates.

Figure 2: General approach to data-driven discovery. ${\bm{\phi}}$,
${\bm{\psi}}$ and $\mathbf{f}$ are the encoder, decoder and analytical model
respectively.

In this work, we leverage the SINDy-autoencoder framework [47] to learn
parsimonious governing dynamical systems from delay-embedding of partial
observations, as summarized in Fig. 2. In general, model discovery from full-
state measurements is an ill-posed inverse problem, and model discovery from
limited measurements compounds this challenge. As with other custom
autoencoders that are used to learn dynamics [44, 47], there are several
critical steps, including choosing the architecture, designing the loss
function, initializing the parameters, and tuning the optimization algorithm
used to train the model. In addition to demonstrating this algorithm to learn
closed-form models from partial observations of several dynamical systems, we
describe how this inverse modeling problem may be broken down into sub-
problems. In general, starting with a problem where the answer is known is
important to guide the design of the architecture and loss functions. Thus, we
show how the architecture and individual terms in the loss function affect the
learning process with different levels of assumed knowledge. We also
investigate what preprocessing should be applied to the data, such as whether
or not to apply the singular value decomposition to the delay coordinates
before the autoencoder network. Incredibly, we show that for the Lorenz
system, from a partial observation of only the first coordinate, it is
possible to learn a coordinate system in which the dynamics may be represented
by a differential equation with even fewer terms than the original Lorenz
system; the original system has 7 linear and quadratic terms, while our model
only has 6 terms. We also show that it is possible to learn a Lorenz analogue
from a video sequence of a chaotic waterwheel experiment.

The rest of the paper is organized as follows. First, we provide a brief
problem definition in Sec. 1.1. Section 2 has background information about
SINDy and how it has been combined with autoencoders. Section 3 discusses the
proposed delay SINDy autoencoder, including the architecture and custom loss
functions. The main results from numerical experiments are provided in Sec. 4.
We analyze several canonical nonlinear and chaotic systems, in addition to
video data from a chaotic waterwheel experiment that may be viewed as a Lorenz
analogue. Finally, conclusions and discussions are provided in Sec. 5

###  1.1 Mathematical problem definition

An overarching inverse problem in data-driven dynamical systems is to learn a
model for the evolution of a system, purely from observation data.
Mathematically, this inverse problem can be stated as follows. Let
$\mathbf{x}(t)$, an $m$ dimensional time-dependent state vector defined on the
domain $\mathcal{D}\in\mathbb{R}^{m}\times\mathbb{R}^{+}$, with
$m\in\mathbb{N}$ and $t>0$, be the solution of a nonlinear differential
equation

$\frac{d}{dt}\mathbf{x}(t)=\mathbf{f}(\mathbf{x}(t);\bm{\mu})$ (1)

where $\mathbf{f}$ is a smooth and nonlinear function and $\bm{\mu}$ a vector
of system parameters. Further, let $\mathbf{y}\in\mathbb{R}^{d}$ represent
noisy measurements of the system, given by

$\mathbf{y}(t)=\mathbf{g}(\mathbf{x}(t))+{\bm{\eta}},$ (2)

where ${\bm{\eta}}$ is a noise process. The goal is then to learn an
approximate dynamical system

$\frac{d}{dt}\mathbf{z}(t)=\check{\mathbf{f}}(\mathbf{z}(t);\check{\bm{\mu}})$
(3)

in terms of a new state $\mathbf{z}$, which may either be the measured state
$\mathbf{y}$ or an invertible function of $\mathbf{y}$:

$\mathbf{z}(t)={\bm{\phi}}(\mathbf{y}(t)).$ (4)

If observation data for the full state $\mathbf{x}(t)$ are available, so that
$\mathbf{y}=\mathbf{x}$ and $\mathbf{z}=\mathbf{y}$, an approximation of
$\mathbf{f}$ can be inferred using a variety of machine learning algorithms.
Although neural networks have been shown to provide useful predictive models
[53, 54, 55, 56, 57], they don’t result in interpretable, closed-form
expressions such as (3) that are amenable to analysis. Parsimonious modeling
techniques, such as genetic programming [1, 2] and SINDy [5, 6], are able to
learn analytic, closed-form models that often capture the exact form of the
original dynamics (1) when they are known.

However, in many real world applications, only partial measurements are
available, so the dimension of $\mathbf{y}$ is less than the dimension of
$\mathbf{x}$, i.e. $d<m$. In general, it is not possible to learn $\mathbf{f}$
directly in its original dimension in terms of the measurement vector
$\mathbf{y}$. A direct approach of fitting a recurrent relation for
$\mathbf{y}$, using recurrent neural networks, Bayesian methods, Gaussian
process regression, etc., or discovering the underlying differential equation
using SINDy, do not generalize well if the dynamics of $\mathbf{y}$ are
nonlinear and depend on hidden variables that control its evolution. The
challenge of representing the dynamics in terms of a partial observation
$\mathbf{y}$ of the full-state $\mathbf{x}$ is known as the _closure_ problem,
and it is one of the most fundamental issues in scientific modeling, with a
particularly rich history in turbulence modeling [58, 24, 59, 60].

There are several embedding techniques [61, 62, 22, 50] to enrich the partial
measurement vector, providing a map as in (4), although it remains a challenge
to learn the dynamics (3) in these coordinates. Therefore, in general, it is
necessary to simultaneously learn the coordinate transformation (4) and the
dynamics (3). For simplicity, we will assume scalar measurements $y(t)$,
although the following arguments generalize to multi-dimensional vectors of
partial measurements. Before learning the map (4), we first time-delay embed
the measurement $y(t)$, resulting in the intermediate vector
$\mathbf{h}(t;n,\tau)=[y(t),y(t+\tau),y(t+2\tau),\ldots,y(t+(n-1)\tau)]$,
where $n$ and $\tau$ are the number of delays and the time increments between
successive measurements respectively. Takens’ theorem [7] provides conditions
for when time-delay embedding results in an attractor that is diffeomorphic to
the original system. We then simultaneously learn the map (4) from
$\mathbf{z}={\bm{\phi}}(\mathbf{h})$ and the dynamics (3) using a modification
of the SINDy-autoencoder framework [36]. The goal is an interpretable
differential equation
$\dot{\mathbf{z}}(t)=\check{\mathbf{f}}(\mathbf{z}(t);\check{\bm{\mu}})$ in
the transformed embedded variable. Takens’ embedding theorem also puts a
condition on the embedding dimension to be $n>2m$ for the diffeomorphism to be
guaranteed.

##  2 Background

In this section, we review time-delay embedding, the sparse identification of
nonlinear dynamics (SINDy) [5] algorithm, and the recent combination of SINDy
with deep autoencoders [36].

#### Delay embedding.

When defining the delay embedding variable $\mathbf{h}(t;n,\tau)$ of a
measurement time series $\tilde{\mathbf{y}}$, the proper choice of $\tau$ and
$n$ have been studied extensively [63, 64] (see appendix A). The embedded
variable $\mathbf{h}(t)$ can be assembled into a Hankel matrix

$\mathbf{H}=\begin{bmatrix}y(t_{1})&y(t_{2})&\ldots&y(t_{q})\\\
y(t_{2})&y(t_{3})&\ldots&y(t_{q+1})\\\ \vdots&\vdots&\ddots&\vdots\\\
y(t_{n})&y(t_{n+1})&\ldots&y(t_{n+q+1})\end{bmatrix}=\begin{bmatrix}\mathbf{h}_{1},\mathbf{h}_{2},\ldots,\mathbf{h}_{q}\end{bmatrix}.$
(5)

The corresponding eigen-time-delays are given by the columns of $\mathbf{U}$
from the singular value decomposition (SVD)
$\mathbf{H}=\mathbf{U}\bm{\Sigma}\mathbf{V}^{\top}$, and they are used for
linear system identification [65, 66, 67, 68, 69] and in the Hankel
alternative view of Koopman (HAVOK) algorithm [22]. The first few eigen-time-
delays, whose rank can be chosen based on a Pareto front analysis, are often
sufficient to encode the shape of a given signal, as shown in Fig. 13.

#### SINDy.

The SINDy algorithm frames the problem of discovering the model in (3) as a
sparse generalized linear regression [5]. SINDy assumes a library of candidate
functions $\bm{\theta}_{i}(\mathbf{z})$ that determine possible terms in
$\check{\mathbf{f}}$. Given $q$ discrete-time samples of
$\mathbf{z}(t)\in\mathbb{R}^{m}$ and $\dot{\mathbf{z}}(t)\in\mathbb{R}^{m}$,
these may be cast into data matrices
$\mathbf{Z}=[\mathbf{z}_{1},\mathbf{z}_{2},\ldots,\mathbf{z}_{q}]^{\top}$ and
$\dot{\mathbf{Z}}=[\dot{\mathbf{z}}_{1},\dot{\mathbf{z}}_{2},\ldots,\dot{\mathbf{z}}_{q}]^{\top}$,
where $\mathbf{Z}$, $\dot{\mathbf{Z}}\in\mathbb{R}^{q\times m}$. Assuming $r$
candidate functions, with $r\ll m$, it is possible to build the library data
matrix
$\bm{\Theta}(\mathbf{Z})=[\bm{\theta}_{1}(\mathbf{Z}),\bm{\theta}_{2}(\mathbf{Z}),\ldots,\bm{\theta}_{r}(\mathbf{Z})]\in\mathbb{R}^{m\times
r}$. The SINDy algorithm then minimizes the loss

$\mathcal{L}_{\text{SINDy}}(\bm{\Xi})=\left\|\dot{\mathbf{Z}}-\bm{\Theta}(\mathbf{Z})\bm{\Xi}\right\|^{2}_{2}+\lambda\left\|\bm{\Xi}\right\|^{2}_{0},$
(6)

where the sparse matrix
$\bm{\Xi}=[\bm{\xi}_{1},\bm{\xi}_{2},\ldots,\bm{\xi}_{r}]\in\mathbb{R}^{r\times
m}$ determines which candidate terms are active in the dynamics, so that
$\check{\mathbf{f}}(\mathbf{z})\approx{\bm{\Theta}}(\mathbf{z}){\bm{\Xi}}$.
The sparsity constraint $\left\|\bm{\Xi}\right\|^{2}_{0}$ minimizes the number
of non-zero terms in $\bm{\Xi}$, but it requires an intractable combinatorial
search. Thus, it is either replaced by $\left\|\bm{\Xi}\right\|^{2}_{1}$ or a
sequentially thresholded least-square algorithm for computational efficiency.
A sparse regularizer is motivated by Occam’s razor, which prefers simpler and
more interpretable models. SINDy has been widely applied to model a range of
complex systems [70, 71, 72, 73, 69, 74, 75, 76, 77, 78, 79, 80, 81].

#### SINDy autoencoders.

SINDy is most effective when the choice of the dictionary is informed by the
structure of the problem. When the potential terms of the differential
equations we seek are not known a priori, or if they include integral
operators (e.g. non-local closures [82]), a naive implementation of SINDy
gives sub-optimal results. A particularly versatile solution is to use
autoencoders to transform the coordinates of the input variable before
discovering a differential equation, as in Champion et al. [36]. The resulting
SINDy autoencoder algorithm is summarized in Fig. 3. We will expand on this
algorithm to address the problem of data-driven discovery of high-dimensional
equations from low-dimensional measurements.

Figure 3: SINDy Autoencoder summary algorithm. ${\bm{\phi}}$ and ${\bm{\psi}}$
are encoder and decoder networks respectively. The network Jacobians
$\partial\mathbf{z}/\partial\mathbf{y}$ and
$\partial\bar{\mathbf{y}}/\partial\bar{\mathbf{z}}$ are computed using
automatic differentiation in tensorflow (i.e. tf.jacobian).

SINDy autoencoders simultaneously transform measurement data
$\mathbf{y}(t)\in\mathbb{R}^{p}$ into a lower dimensional variable
$\mathbf{z}(t)\in\mathbb{R}^{m}$ and discover a nonlinear differential
equation $\dot{\mathbf{z}}=\bm{\Theta}(\mathbf{z}){\bm{\Xi}}$ in $\mathbf{z}$.
The main element of the algorithm is an autoencoder that requires a
reconstruction loss

$\mathcal{L}_{\text{recon}}=\left\|{\bm{\psi}}({\bm{\phi}}(\mathbf{y}))-\mathbf{y}\right\|_{2}^{2}.$
(7)

In addition, the method incorporates the SINDy loss given in equation (6) in
two different ways. First, the difference between the output of the SINDy
layer,
$\dot{\bar{\mathbf{z}}}=\bm{\Theta}(\mathbf{z})\bm{\Xi}=\bm{\Theta}({\bm{\phi}}(\mathbf{y}))\bm{\Xi}$,
and that computed by the chain rule
$\dot{\mathbf{z}}=\nabla_{\mathbf{y}}{\bm{\phi}}(\mathbf{y})\dot{\mathbf{y}}$,
is minimized. This SINDy loss in $\dot{\mathbf{z}}$ is given by

$\mathcal{L}_{\dot{\mathbf{z}}}=\left\|\nabla_{\mathbf{y}}{\bm{\phi}}(\mathbf{y})\dot{\mathbf{y}}-\bm{\Theta}({\bm{\phi}}(\mathbf{y}))\bm{\Xi}\right\|_{2}^{2}.$
(8)

Second, a consistency loss ensuring that $\dot{\mathbf{y}}$ can be
reconstructed from $\mathbf{y}$ using the path
$\mathbf{y}\rightarrow\mathbf{z}\rightarrow\dot{\bar{\mathbf{z}}}\rightarrow\dot{\bar{\mathbf{y}}}$
is included. Using the fact that

$\displaystyle\dot{\bar{\mathbf{y}}}$
$\displaystyle=\frac{\partial\bar{\mathbf{y}}}{\partial\bar{\mathbf{z}}}\cdot\frac{\partial\bar{\mathbf{z}}}{\partial\mathbf{t}}$
$\displaystyle=\nabla_{\mathbf{z}}{\bm{\psi}}(\mathbf{z})\cdot\bm{\Theta}(\mathbf{z})\bm{\Xi}$
$\displaystyle=\nabla_{\mathbf{z}}{\bm{\psi}}({\bm{\phi}}(\mathbf{y}))\cdot\bm{\Theta}({\bm{\phi}}(\mathbf{y}))\bm{\Xi},$

the resulting loss in $\dot{\mathbf{y}}$ is given by

$\mathcal{L}_{\dot{\mathbf{y}}}=\left\|\dot{\mathbf{y}}-\nabla_{\mathbf{z}}{\bm{\psi}}({\bm{\phi}}(\mathbf{y}))\bm{\Theta}({\bm{\phi}}(\mathbf{y}))\bm{\Xi}\right\|_{2}^{2}.$
(9)

Finally, a sparsity constraint $\mathcal{L}_{\text{sparse}}=\|\bm{\Xi}\|_{1}$
is included in the loss function to encourage the discovery of a parsimonious
model. While the choice of these losses is sufficient for some high-
dimensional data, the current problem requires more loss terms to account for
delayed inputs.

##  3 Delay SINDy Autoencoders

We now leverage the SINDy autoencoders algorithm to learn coordinates and
dynamics starting from a delay-embedding. The algorithm is extended by adding
losses that take advantage of time-delay embedding properties.

Let $\\{\mathbf{h}_{i}\\}_{i=1}^{q}$ be the embedding of discrete time-series
data $\hat{\mathbf{y}}=[y(\tau),y(2\tau),y(3\tau),\cdots]$, given by the
column vectors of the Hankel matrix in (5). Accordingly, $h_{ij}$ is the
$j$-th entry of the embedding vector $\mathbf{h}_{i}$. Our approach relies on
the following premises:

  1. 1.

Takens’ embedding theorem proves that a diffeomorphic map between the delay-
embedded attractor $\mathbf{h}$ and the original attractor $\mathbf{x}(t)$
exists under certain conditions. The theorem doesn’t specify how to find this
map.

  2. 2.

The universal approximation properties of neural networks [83] can approximate
this diffeomorphism and others where the dynamics are simplified.

The goal is then to design a loss function whose optimal solution is a
predictive model with a sparse analytical form. Ideally, the optimal sparse
differential equation of the latent variable, $\mathbf{z}(t)$, is similar to
that of $\mathbf{x}(t)$ when it is known. Finding the best loss function is
matter of making the hypothesis class as small as possible and as close as
possible to the optimal mappings in such a way to minimize the approximation
and estimation errors. This is illustrated in Fig. 4. We explore these
premises in Sec. 4.

Figure 4: Illustration of the estimation and approximation errors in the
hypothesis class.

In addition to the SINDy autoencoder losses in Sec. 2, we introduce three
losses that take advantage of a priori knowledge of time-delay embedding
properties, as illustrated in Fig. 5.

#### Loss 1.

First, the dimensionality reduction of $\mathbf{h}$ can be performed using
linear techniques, such as computing the SVD of the Hankel matrix in (5), or
using nonlinear techniques, such as a deep autoencoder network. The
diffeomorphism between $\mathbf{h}$ and $\mathbf{x}$ is expected to be highly
nonlinear, which motivates the use of neural network. However, it is possible
to perform an initial linear dimensionality reduction, using the first $p$
dominant modes of the Hankel matrix as an input to the autoencoder. This step
can filter out noise and improve the computational efficiency of the method.
If $n\gg 2m$, then a linear projection of $\mathbf{h}$ on an intermediate
dimension $p$ preserves the majority of the variance in the signal.

Figure 5: Summary of delay SINDy Autoencoders algorithm. The tensor $h_{ij}$
denotes the $i$ row and $j$ column of the Hankel matrix $\mathbf{H}$.

In this case, the SVD projection $\mathbf{h}^{\top}\mathbf{U}_{p}$ and its
inverse are included as part of the encoding and decoding of the time-series
$y(t)$. We define $\mathbf{U}_{p}$ to be the first $p$ columns of
$\mathbf{U}$, and $\bm{\Sigma}_{p}$ the first $p$ rows and columns of
$\bm{\Sigma}$ from the SVD of $\mathbf{H}$. The SINDy losses in
$\dot{\mathbf{z}}$ and $\dot{\mathbf{y}}$, given in the original SINDy
autoencoder algorithm by (8) and (9) respectively, now become

$\displaystyle\mathcal{L}_{\dot{\mathbf{z}}}$
$\displaystyle=\left\|\nabla_{\mathbf{h}}\bm{\phi}(\mathbf{h}^{\top}\mathbf{U}_{p}))\dot{\mathbf{h}}-\bm{\Theta}({\bm{\phi}}(\mathbf{h}^{\top}\mathbf{U}_{p}))\bm{\Xi}\right\|_{2}^{2}$
(10a) $\displaystyle\mathcal{L}_{\dot{\mathbf{h}}}$
$\displaystyle=\left\|\dot{\mathbf{h}}-\nabla_{\mathbf{z}}\left[\bm{\psi}(\bm{\phi}(\mathbf{h}^{\top}\mathbf{U}_{p}))\mathbf{U}_{p}^{\top}\right]\bm{\Theta}(\bm{\phi}(\mathbf{h}^{\top}\mathbf{U}_{p}))\bm{\Xi}\right\|_{2}^{2},$
(10b)

where $\mathbf{y}$ has been replaced by $\mathbf{h}$ for consistency.
Furthermore, the reconstruction loss is obtained from the difference between
the feed-forward prediction $\overline{\mathbf{h}}$ and the input $\mathbf{h}$

$\mathcal{L}_{\text{recon}}=\left\|\mathbf{h}-{\bm{\psi}}({\bm{\phi}}(\mathbf{h}^{\top}\mathbf{U}_{p}))\mathbf{U}_{p}^{\top}\right\|_{2}^{2}.$
(11)

#### Loss 2.

Second, we require the first component of the discovered latent variable
$\mathbf{z}(t)$ to be an exact reconstruction of the input signal $y(t)$. This
significantly constrains the hypothesis class by an $m$-fold reduction of the
number of permutations of the components of $\mathbf{z}$, to $m-1!$
possibilities instead of $m!$. This reduction ultimately makes the
optimization problem more feasible. The first component loss is simply given
by

$\mathcal{L}_{z_{1}}=\left\|h_{i1}-z_{i1}\right\|^{2}_{2},$ (12)

where $\mathbf{h}_{i}$ is the $i$-th delay embedding of $\hat{\mathbf{y}}$ and
$\mathbf{z}_{i}=[z_{i1},z_{i2},\ldots,z_{im}]$ is the corresponding latent
variable transformation.

In addition to favorably constraining the optimization, this loss function
promotes a more physically interpretable solution. The variable $y(t)$
typically has some physical meaning, as it is the quantity that is actually
measured experimentally. Thus, fixing the first component of $\mathbf{z}$ to
be equal to $y(t)$ guarantees that the dynamics in the latent space are
directly relevant for the measured quantity $y$.

#### Loss 3.

Finally, if the constraint in equation (12) is satisfied, the first dimension
of the time integration of $\mathbf{z}_{i}$ using the SINDy model,
$\dot{\mathbf{z}}_{i}=\bm{\Theta}(\mathbf{z}_{i})\bm{\Xi}$, should correspond
to the vector $\mathbf{h}_{i}$. In Fig. 5, we set
$\hat{\mathbf{z}}_{i}(t_{1})\equiv\mathbf{z}_{i}$ as an initial condition
where $t_{1}$ is the initial time, and integrate it using the SINDy model with
a time step $\Delta t=t_{i}-t_{i-1}$ that matches the sampling rate of the
input data $\hat{\mathbf{y}}$. The integration can be expressed as

$\hat{\mathbf{z}}_{i}(t_{i})=\int_{t_{1}}^{t_{i}}\bm{\Theta}({\bm{\phi}}(\mathbf{h}_{i}^{\top}\mathbf{U}_{p}))\bm{\Xi}\,dt.$
(13)

The corresponding SINDy consistency loss is therefore given by

$\mathcal{L}_{\text{cons}}=\sum_{j=1}^{n}\left\|h_{ij}-\left(\int_{t_{1}}^{t_{j}}\bm{\Theta}({\bm{\phi}}(\mathbf{h}_{i}^{\top}\mathbf{U}_{p}))\bm{\Xi}dt\right)_{1}\right\|_{2}^{2},$
(14)

where the subscript $1$ of the parentheses designates the first component of
$\hat{\mathbf{z}}(t_{i})$, and $t_{j}=t_{1}+j\Delta t$ for $j\in[1..n-1]$. The
initial condition $\hat{\mathbf{z}}(t_{1})$ is not included in the loss to
avoid redundancy with the $\mathcal{L}_{z_{1}}$ loss. The integration
procedure is included in the network as a recurrent layer with an RK4
integration scheme.

Finally, all the losses are combined to give

$\mathcal{L}=\mathcal{L}_{\text{recon}}+\lambda_{1}\mathcal{L}_{\dot{\mathbf{h}}}+\lambda_{2}\mathcal{L}_{\dot{\mathbf{z}}}+\lambda_{3}\mathcal{L}_{z_{1}}+\lambda_{4}\mathcal{L}_{\text{cons}}+\lambda_{5}\mathcal{L}_{\text{reg}},$
(15)

with weighting coefficients
$\bm{\lambda}=[\lambda_{1},\lambda_{2},\ldots,\lambda_{5}]$ that are
hyperparameters to be optimized.

The proposed algorithm constrains the hypothesis class to make the inverse
problem computationally feasible and the corresponding optimal solution
predictive of the measurement data. However, the choice of the weights
$\bm{\lambda}$ and the losses to include are not unique. Given more
information about the system, such as conservation of momentum or energy, more
losses can be included to guide the optimization algorithm to find the optimal
transform within the hypothesis class.

##  4 Results

We demonstrate the ability of the proposed deep delay autoencoder network to
discover coordinates and parsimonious models for a number of canonical
dynamical systems. Simultaneously learning the coordinates and dynamics is a
challenging inverse problem. To solve this inverse modeling problem, and to
give insight into potential pitfalls, we break it down into the following sub-
problems, of increasing difficulty:

  1. 1.

Optimize ${\bm{\phi}}$ and ${\bm{\psi}}$ (with and without an intermediate SVD
layer) with known $\mathbf{x}$; i.e. using supervised learning.

  2. 2.

Optimize ${\bm{\phi}}$ and ${\bm{\psi}}$ (with and without an SVD layer) with
known analytical model $\mathbf{f}$, by fixing the coefficient matrix
${\bm{\Xi}}$.

  3. 3.

Optimize ${\bm{\phi}}$, ${\bm{\psi}}$, and ${\bm{\Xi}}$, initializing
$\mathbf{f}$ in the proximity of its original value in a known system of
equations.

  4. 4.

Optimize ${\bm{\phi}}$, ${\bm{\psi}}$, and ${\bm{\Xi}}$ with random
initialization.

The challenge of solving these optimization problems, in increasing order of
difficulty, is due to the size of the hypothesis class which is spanned by the
fitting parameters of ${\bm{\phi}}$, ${\bm{\psi}}$, and ${\bm{\Xi}}$, as
illustrated in Fig. 4. Increasing the complexity of the hypothesis decreases
the approximation error by capturing the nonlinear diffeomophisms
${\bm{\phi}}$ and ${\bm{\psi}}$. In the chaotic systems examples used in this
section, we expect $\mathbf{f}$ to be nonlinear, with at least second order
polynomial features. However, expanding the hypothesis class makes the
optimization task more difficult, thus increasing the estimation error that
arises because of sub-optimal fitting (e.g. local minimum, over-fitting, lack
of data etc.). Accordingly, by assuming some fitting functions to be known
prior to optimization, we slowly increase the size of the hypothesis class to
explore its robustness and generalization.

We apply the full algorithm to the following systems: the Lorenz system, the
Lotka-Volterra equations, the Rössler attractor, and a video of the Lorenz
waterwheel. In the following examples, the first component of $\mathbf{x}$ is
used for the partial measurement $y(t)$. However, not all scalar measurements
$y(t)$ of $\mathbf{x}$ are equally fit to discover the mappings ${\bm{\psi}}$
and ${\bm{\phi}}$. For example, measurement of the third component of the
Lorenz system does not satisfy Takens’ embedding theory, and therefore, the
delay-embedded attractor is not diffeomorphic to the full-state attractor
$\mathbf{x}$. Furthermore, the original full-state system may not be the only
coordinate system that admits a sparse representation of the dynamics, and
indeed there may be even sparser representations, as we will show for the
Lorenz system. Therefore, the maps ${\bm{\phi}}$ and ${\bm{\psi}}$ and the
dynamics $\check{\mathbf{f}}$ may not be unique, which we will discuss in the
following test cases.

The delay embedding dimension $n$ and the time step $\tau$ are chosen to
satisfy the equality $n\tau\approx 0.1$, such that the delay embedding of
$y(t)$ is most unfolded in the embedding phase space (see Appendix A); this
choice is also consistent with that from HAVOK [22]. Along with $n$ and
$\tau$, the hyperparameters over which we optimize ${\bm{\phi}}$,
${\bm{\psi}}$, and $\mathbf{\Xi}$ are the SVD rank dimension $p$, the
dimension of the latent variable $m$, the architecture of the neural network
(number of nodes, number of layers, activation function etc.), and the SINDy
library (polynomial order, trigonometric features, threshold value, etc.). The
high dimensionality of the hyperparameter space adds to the challenge of the
optimization problem, which we address in the current study by running a
hyperparameter search on parallel GPU clusters.

###  4.1 The Lorenz System

As a benchmark problem, we investigate the three dimensional Lorenz system,
which arises in atmospheric modeling and is often used to study chaos and
nonlinear dynamics. The system of equations for $\mathbf{x}(t)$ is given by

$\displaystyle\dot{x}_{1}$ $\displaystyle=\sigma(x_{2}-x_{1})$ (16a)
$\displaystyle\dot{x}_{2}$ $\displaystyle=x_{1}(\rho-x_{3})-x_{2}$ (16b)
$\displaystyle\dot{x}_{3}$ $\displaystyle=x_{1}x_{2}-\beta x_{3},$ (16c)

where one typically assumes $\sigma$, $\rho$, and $\beta$ are positive. The
system is known to be chaotic for $\sigma=10$, $\rho=28$, and $\beta=8/3$. In
this experiment, we use simulation data of the first component $x_{1}$ as the
measurement $y(t)$, which is input to the neural network. The goal is to find
a coordinate transformation from the delay-embedded coordinates $\mathbf{h}$
to a new coordinate system $\mathbf{z}$ where the dynamics have a sparse,
closed-form representation.

Figure 6 shows that sub-problems 1 and 2 can be solved with high accuracy,
when either the full state $\mathbf{x}$ or the dynamics $\mathbf{f}$ are
known, respectively. For small $p=10$, sub-problem 2 has a large error for
large amplitudes of $x_{3}$ where the trajectory switches lobes. This example
also explores performance with and without the use of an intermediate SVD
step.

Figure 6: Learning an autoencoder map from delay embedding to the full-state
attractor, either with a known $\mathbf{x}$ (supervised) or with a known full-
state model $\mathbf{f}$ (known equation). We investigate the performance with
and without using an intermediate SVD layer. The parameters are $n=128$,
$p=10$, and $m=3$, using $3$ hidden layers in both the encoder and decoder.
The computational advantage of using SVD comes at the cost of lower prediction
accuracy.

The diffeomorphic transformations ${\bm{\phi}}$ and ${\bm{\psi}}$ are
approximated in consecutive nonlinear stages through the hidden layers of the
autoencoder. Figure 7 illustrates the evolution of the attractor as a function
of network layers (column-wise) and training epochs (row-wise) when the
analytical model is known and fixed in the SINDy coefficients $\mathbf{\Xi}$.
We use the three dominant modes of every layer for visualization and
initialize the network to map $\mathbf{h}$ to
$\mathbf{v}=[v_{1},v_{2},v_{3}]$, where $\mathbf{v}$ are the three dominant
right singular vector modes of the Hankel matrix $\mathbf{H}$. This
initialization of ${\bm{\phi}}$ and ${\bm{\psi}}$ results in a solution that
is closer to the optimal solution than random initialization, particularly by
making them smooth functions. This is a form of transfer learning in the
context of data-driven modeling.

Figure 7: Evolution of the attractor as a function of epochs and deep network
layers for the known equation case through the encoder. The three linear
dominant modes of every network layer are shown for visualization. Columns
correspond to layers row to selected epochs. At epoch 0, the network is
initialized to map the embedded data $\mathbf{h}$ to the three dominant modes
of $\mathbf{H}$ at the output layer.

Initializations of $\mathbf{\Xi}$ with large coefficient perturbations
(sampled from a normal distribution $\mathcal{N}(\xi_{ij},20)$) around the
original coefficients of $\mathbf{f}$ consistently lead to Lorenz-like chaotic
attractors with a two-lobe structure. However, the convergence of the
algorithm for random initializations of $\mathbf{\Xi}$ is sensitive to
hyperparameters and does not always give rise to a model that has a similar
butterfly structure as the original Lorenz attractor. This sporadic
convergence indicates that the optimizer is prone to over-fitting and getting
stuck in local minima, an issue that is commonly observed in deep learning. We
find that applying a separate SINDy optimization on $\mathbf{z}$ every $k$
epochs (with $k\approx 20-40$) improves sparsity and efficiency significantly,
while nudging the fitting parameters just enough to avoid local minima. For a
given optimal set of weighting coefficients $\bm{\lambda}$, we discover a
sparse system of equations

$\displaystyle\dot{z}_{1}$ $\displaystyle=-16z_{1}+2.5z_{2}z_{3}$ (17a)
$\displaystyle\dot{z}_{2}$ $\displaystyle=8z_{2}-16z_{1}z_{3}$ (17b)
$\displaystyle\dot{z}_{3}$ $\displaystyle=8z_{2}^{2}-1.4z_{3}^{2},$ (17c)

whose corresponding attractor is similar to the original system, as shown in
Fig. 9. The discovered latent representation of the time-delay embedded data
is a almost a scaling and translation transform of the original attractor.
This discrepancy is expected from a method that considers linear transforms
and rotation of a given dynamical systems to be equally valid diffeomorphic
solutions. Note that this Lorenz-like system is more sparse than the original
system of equations, containing 6 instead of 7 terms, but with 4 nonlinear
terms instead of 2. The coefficients of equation 17 also appear to be involve
3 system parameters ($a=8$, $b=2.5$ and $c=1.4$), similar to the original
system with parameters $\rho$, $\beta$ and $\sigma$.

Figure 8: Coefficients evolution as a function of epochs, where the dashed
lines are the coefficients of the original system and the solid lines are the
evolution of the discovered coefficients (left), $z_{1}$ loss and
reconstruction loss as a function of epochs (bottom right), and input
reconstructions as a function of time (top right).

Figure 8 shows the evolution of the coefficients as a function of time, the
autoencoder reconstruction comparing $\mathbf{h}$ and $\bar{\mathbf{h}}$
(optimized through equation (11)), and a comparison between the first
components of $\mathbf{z}$ and $\mathbf{\hat{h}}$ (optimized through equation
(12)). The input $\mathbf{h}$ is nearly perfectly reconstructed at the output
of the autoencoder and in the first dimension of the latent variable, with
loss below $10^{-5}$. The coefficients of the discovered model diverge from
the original coefficients, although with a similar two-lobe attractor
dynamics.

Figure 9: Results for Lorenz, Rössler, and Lotka Volterra systems.

###  4.2 Other systems

We now test the algorithm on the Lotka-Volterra model and the Rössler system,
as well as on a challenging example of a video of a chaotic waterwheel
experiment.

#### Lotka Volterra.

The Lotka-Volterra equations describe the interaction of multiple species with
a predator-prey relationship. In the simplest case, a system of equations that
quantifies the number of a predator $x$ and prey $y$ is given by

$\displaystyle\dot{x}$ $\displaystyle=ax+bxy$ (18a) $\displaystyle\dot{y}$
$\displaystyle=cy+dxy,$ (18b)

where $a$, $b$, $c$, and $d$ are parameters that depend on the death rate and
reproduction rate of both species. In this case, the attractor is a family of
limit cycles that depend on the input parameters and the initial condition.
The discovered models are structurally similar as shown in Fig. 9.

#### Rössler system.

The Rössler system is a three-dimensional system of differential equations,
commonly used for studying chaotic systems, given by

$\displaystyle\dot{x}_{1}$ $\displaystyle=-x_{2}-x_{3}$ (19a)
$\displaystyle\dot{x}_{2}$ $\displaystyle=-x_{1}+ax_{2}$ (19b)
$\displaystyle\dot{x}_{3}$ $\displaystyle=a+x_{3}(x_{1}-b),$ (19c)

where $a$, $b$ are input parameters. The attractor is known for its
intermittently bursting third dimension. Figure 9 shows that the discovered
model is less sparse than the original system; however, the bursting effect in
the third dimension is accurately captured.

Figure 10: Lorenz waterwheel discovery from video [84]. An OpenCV line
detection algorithm is used to find the angle of the spokes and a k-means
algorithm is used for clustering lines that detect the same spoke. The angles
of the lines are post-processed to deduce one angular velocity $\omega$ for
each frame, which is used as an input measurement variable.

#### Chaotic waterwheel.

Finally, we use video footage of the Lorenz waterwheel experiment to discover
a latent representation and its corresponding system of equations, as shown in
Fig. 10. The video is encoded into a time series using the Hough transform to
detect the spokes of the wheel, a k-means clustering algorithm to group lines
that identify the same spoke, which are post-processed to calculate the
angular velocity of the wheel. A latent representation of the Lorenz-like
dynamics of the angular velocity is discovered by the delay SINDy autoencoder,
and third degree polynomial features are required in the model for chaotic
dynamics (see Appendix B). The rough attractor shape and chaotic dynamics are
both reproduced. The evolution of the losses as a function of epochs is shown
in Fig. 11.

Figure 11: Evolution of the losses as a function of epochs for the Lorenz
waterwheel case.

##  5 Discussion

This work has shown that it is possible to discover a sparse, interpretable,
nonlinear dynamical system from incomplete or partial measurements using a
combination of data augmentation through time-delay embedding, nonlinear
dimensionality reduction and coordinate transformation with deep neural
network autoencoders, and model identification with SINDy. Thus, this work
begins to resolve a forty year old open challenge, which began with Takens’
embedding theorem, to explicitly learn a diffeomorphism back from delay
coordinates to a coordinate system where the dynamics may be efficiently
represented. We also highlight the importance of choosing the proper
hypothesis class, as constrained by the relevant loss function. While the
resulting analytical reconstructions do not always match the original system
from which the data was generated, the discovered latent variables and their
underlying models are sparse, sharing key structural features with the
original system when it is known. Surprisingly, we discover Lorenz-like models
that are more sparse than the original Lorenz equations from 1D measurements.
Our models have only six terms, while the original Lorenz equations have seven
terms; however, additional sparsity comes at the cost of a higher number of
quadratic terms. The method’s ability to discover models that are
qualitatively and quantitatively similar to the original system, when the
differential equation is initialized with perturbations around the original
coefficients, suggests that the optimizer is prone to over-fitting and
becoming stuck in local minima, resulting in large estimation errors.
Designing and adding appropriate regularization terms to the loss function to
solve this issue remains an open question and is the subject of further
investigation.

It is unknown whether or not there exists a loss function for which the
original system of equations can be exactly recovered from one-dimensional
measurements and random fitting parameter initialization. There is always the
possibility of adding more physics-informed constraints to the loss function,
such as known problem-specific symmetries, to help the optimizer find the
desired model. Furthermore, the diffeomorphic transform between the delay
embedding and the full-state attractor can be further enforced by using a
bijective neural network, i.e. an invertible neural network, that honors the
one-to-one map by construction. In addition, the SINDy algorithm used in this
study is relatively basic in its use of the $L_{1}$ loss for enforcing
sparsity, and many improvements have been suggested in the past few years [85,
86], which can be integrated into the proposed delay SINDy autoencoder
framework. Particularly, using a probabilistic [87] or a weak [70, 75, 77, 81]
SINDy algorithm would certainly improve on the results obtained from real-
world noisy data such as the Lorenz waterwheel video.

Crucially, this study proposes a general framework for discovering multi-
dimensional models from partial, low-dimensional data, motivating many
theoretical and numerical investigations into the combined use of time-delay
embedding, SINDy, and autoencoders for data-driven modeling. Future work will
involve exploring various deep network architectures, such as convolutional
autoencoders, for discovering physics from high-dimensional video data.

## Acknowledgments

The authors acknowledge support from the Army Research Office (ARO
W911NF-19-1-0045) and the National Science Foundation AI Institute in Dynamic
Systems (Grant No. 2112085).

##  Appendix A Delay embedding hyperparameters

The time delay increment $n\tau$, where $\tau$ is the sample time and $n$ is
the number of delays included in the Hankel matrix, is typically chosen to
make the attractor occupy (or “unfold” in) as much of the embedded phase space
as possible. Methods for choosing $n\tau$ using measures of autocorrelation or
mutual information [88] to maximize the phase-space unfolding have been
studied extensively [15, 63, 89]. Sampling strategies have also been studied
for discovering physics from data [47]. In this study, we choose the embedding
dimension $n$ and the time delay increment $n\tau$ to match that chosen in
HAVOK [22]. A comparison of the unfolding of the attractor based on $n\tau$
alone in figure 12 shows that the optimal unfolding is for $n\tau\approx 0.1$.

\begin{overpic}[trim=0.0pt 0.0pt 0.0pt
40.15pt,clip,width=398.33858pt]{figures/effect_of_tau.pdf}
\put(6.0,28.0){$n\tau=0.01$} \put(32.5,28.0){$n\tau=0.03$}
\put(61.0,28.0){$n\tau=0.10$} \put(87.0,28.0){$n\tau=0.20$} \end{overpic}
Figure 12: The delay time increment $n\tau$ is chosen to optimally unfold the
attractor in the embedding space. In the case of the Lorenz attractor,
$n\tau\approx 0.1$ seems to be the most appropriate choice.

When finding the dominant modes of the Hankel matrix, similar considerations
for the dimension $n$ have to be taken into account. Figure 13 shows that the
first three dominant modes are representative of the majority of the variance
in the data. This helps guide the choice of the latent dimension $m$ when it
is not known a priori. More generally, standard techniques for determining the
number of autoencoder latent variables may be used before adding the SINDy
losses.

Figure 13: The first 3 SVD modes of the Hankel matrix reflect the two-lobe
structure and capture a majority of the variance.

It is reasonable to ask if the the learning of coordinates and identification
of models can be separated into sequential learning steps. Fig. 14 shows the
latent representation $\mathbf{z}$ from an autoencoder without SINDy
constraints. Because the embedding is non-unique, the autoencoder will
typically fail to capture essential symmetries from the original system,
making it difficult or impossible to achieve sparse models, thus motivating
the combined SINDy-autoencoder approach.

Figure 14: Autoencoder latent representation with only a reconstruction loss
(without SINDy). Error is $\varepsilon\sim 7\cdot 10^{-7}$. This coordinate
system is unlikely to yield a sparse model.

##  Appendix B Lorenz waterwheel model

The discovered Lorenz waterwheel attractor shown in Fig. 10 has the same two-
lobe structure as the original Lorenz system. However, the sparsity of the
data, the measurement noise that arises from the line detection algorithm, the
small frame rate of the video (24 frames/second) and its limited time span
make the sparse identification of the differential equation a challenging
task. After optimizing over hyperparameters, a third order polynomial
dictionary (with $n=32$ and a $4$ layer encoder/decoder) is found to recreate
the chaotic dynamics of the Lorenz system. The discovered coefficients are
shown in table 1, where the differential equation is written as

$\dot{z}_{i}(t)=\sum_{j}\xi_{ij}\theta_{j}.$ (20)

A sparser and more generalizing analytical model can most likely be found by
collecting more data and improving the algorithm that maps the video to the
angular velocity. In future studies, we will explore machine learning
techniques (e.g. using a CNN) to incorporate the video to angular velocity
mapping as part of the SINDy autoencoder algorithm.

$\theta$ | $1$ | $z_{1}$ | $z_{2}$ | $z_{3}$ | $z_{1}^{2}$ | $z_{1}z_{2}$ | $z_{1}z_{3}$ | $z_{2}^{2}$ | $z_{1}z_{2}$ | $z_{3}^{2}$ | $z_{1}^{3}$ | $z_{1}^{2}z_{2}$ | $z_{1}^{2}z_{3}$ | $z_{1}z_{2}^{2}$ | $z_{1}z_{2}z_{3}$ | $z_{1}z_{3}^{2}$ | $z_{2}^{3}$ | $z_{2}^{2}z_{3}$ | $z_{2}z_{3}^{2}$ | $z_{3}^{3}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$\dot{z}_{1}$ | -0.12 | 1.70 | 0.04 |  |  | -1.79 |  | 0.16 |  |  |  |  |  |  |  |  | -0.02 | 0.03 | -0.37 |
$\dot{z}_{2}$ | -1.05 | 0.05 | -0.28 |  | 0.12 |  | 1.50 | 0.60 |  |  |  |  |  |  | -0.10 |  |  | 0.18 | -0.33 |
$\dot{z}_{3}$ | 0.08 | -0.12 | 0.30 | 0.07 | 0.08 | -0.29 | -0.22 | 0.35 | -0.55 |  |  | 0.06 |  | -0.07 | 0.20 | 0.03 | -0.25 | 0.28 |  |
Table 1: Discovered coefficients of the Lorenz waterwheel system.

## References

  * [1] Michael Schmidt and Hod Lipson.  Distilling free-form natural laws from experimental data.  science, 324(5923):81–85, 2009.
  * [2] Josh Bongard and Hod Lipson.  Automated reverse engineering of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 104(24):9943–9948, 2007.
  * [3] Michael D Schmidt, Ravishankar R Vallabhajosyula, Jerry W Jenkins, Jonathan E Hood, Abhishek S Soni, John P Wikswo, and Hod Lipson.  Automated refinement and inference of analytical models for metabolic networks.  Physical biology, 8(5):055011, 2011.
  * [4] Bryan C Daniels and Ilya Nemenman.  Automated adaptive inference of phenomenological dynamical models.  Nature communications, 6, 2015.
  * [5] Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the national academy of sciences, 113(15):3932–3937, 2016.
  * [6] Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Data-driven discovery of partial differential equations.  Science Advances, 3(4):e1602614, 2017.
  * [7] Floris Takens.  Detecting strange attractors in turbulence.  In Dynamical systems and turbulence, Warwick 1980, pages 366–381. Springer, 1981.
  * [8] Agustín Somacal, Yamila Barrera, Leonardo Boechi, Matthieu Jonckheere, Vincent Lefieux, Dominique Picard, and Ezequiel Smucler.  Uncovering differential equations from data with hidden variables.  arXiv preprint arXiv:2002.02250, 2020.
  * [9] Steven Atkinson.  Bayesian hidden physics models: Uncertainty quantification for discovery of nonlinear partial differential operators from data.  arXiv preprint arXiv:2006.04228, 2020.
  * [10] H Ribera, S Shirman, AV Nguyen, and NM Mangan.  Model selection of chaotic systems from data with hidden variables using sparse data assimilation.  arXiv preprint arXiv:2105.10068, 2021.
  * [11] J-C Roux, Reuben H Simoyi, and Harry L Swinney.  Observation of a strange attractor.  Physica D: Nonlinear Phenomena, 8(1-2):257–266, 1983.
  * [12] Dr S Broomhead and Gregory P King.  Extracting qualitative dynamics from experimental data.  Physica D: Nonlinear Phenomena, 20(2-3):217–236, 1986.
  * [13] James P Crutchfield and Bruce S McNamara.  Equations of motion from a data series.  Complex systems, 1:417–452, 1987.
  * [14] DS Broomhead and R Jones.  Time-series analysis.  In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, volume 423, pages 103–121. The Royal Society, 1989.
  * [15] Matthew B Kennel and Steven Isabelle.  Method to distinguish possible chaos from colored noise and to determine embedding parameters.  Physical Review A, 46(6):3111, 1992.
  * [16] Ethan R Deyle and George Sugihara.  Generalized theorems for nonlinear state space reconstruction.  PloS one, 6(3):e18295, 2011.
  * [17] George Sugihara, Robert May, Hao Ye, Chih-hao Hsieh, Ethan Deyle, Michael Fogarty, and Stephan Munch.  Detecting causality in complex ecosystems.  Science, 338(6106):496–500, 2012.
  * [18] Dimitrios Giannakis and Andrew J Majda.  Nonlinear laplacian spectral analysis for time series with intermittency and low-frequency variability.  Proceedings of the National Academy of Sciences, 109(7):2222–2227, 2012.
  * [19] Stephen A Billings.  Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains.  John Wiley & Sons, 2013.
  * [20] Hao Ye, Ethan R Deyle, Luis J Gilarranz, and George Sugihara.  Distinguishing time-delayed causal interactions using convergent cross mapping.  Scientific reports, 5(1):1–9, 2015.
  * [21] Hao Ye, Richard J Beamish, Sarah M Glaser, Sue CH Grant, Chih-hao Hsieh, Laura J Richards, Jon T Schnute, and George Sugihara.  Equation-free mechanistic ecosystem forecasting using empirical dynamic modeling.  Proceedings of the National Academy of Sciences, 112(13):E1569–E1576, 2015.
  * [22] Steven L Brunton, Bingni W Brunton, Joshua L Proctor, Eurika Kaiser, and J Nathan Kutz.  Chaos as an intermittently forced linear system.  Nature communications, 8(1):1–9, 2017.
  * [23] J.-C. Loiseau, B. R. Noack, and S. L. Brunton.  Sparse reduced-order modeling: sensor-based dynamics to full-state estimation.  Journal of Fluid Mechanics, 844:459–490, 2018.
  * [24] Shaowu Pan and Karthik Duraisamy.  Data-driven discovery of closure models.  SIAM Journal on Applied Dynamical Systems, 17(4):2381–2413, 2018\.
  * [25] Suddhasattwa Das and Dimitrios Giannakis.  Delay-coordinate maps and the spectra of Koopman operators.  Journal of Statistical Physics, 175(6):1107–1145, 2019.
  * [26] Dimitrios Giannakis.  Data-driven spectral decomposition and forecasting of ergodic dynamical systems.  Applied and Computational Harmonic Analysis, 47(2):338–396, 2019\.
  * [27] Jakob Runge, Sebastian Bathiany, Erik Bollt, Gustau Camps-Valls, Dim Coumou, Ethan Deyle, Clark Glymour, Marlene Kretschmer, Miguel D Mahecha, Jordi Muñoz-Marí, et al.  Inferring causation from time series in earth system sciences.  Nature communications, 10(1):1–13, 2019.
  * [28] Yong Zou, Reik V Donner, Norbert Marwan, Jonathan F Donges, and Jürgen Kurths.  Complex network approaches to nonlinear time series analysis.  Physics Reports, 787:1–97, 2019.
  * [29] Shaowu Pan and Karthik Duraisamy.  On the structure of time-delay embedding in linear models of non-linear dynamical systems.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(7):073135, 2020.
  * [30] Dimitrios Giannakis.  Delay-coordinate maps, coherence, and approximate spectra of evolution operators.  arXiv preprint arXiv:2007.02195, 2020.
  * [31] Seungjoon Lee, Felix Dietrich, George E Karniadakis, and Ioannis G Kevrekidis.  Linking gaussian process regression with data-driven manifold embeddings for nonlinear data fusion.  Interface focus, 9(3):20180083, 2019.
  * [32] Christoph Wehmeyer and Frank Noé.  Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinetics.  The Journal of Chemical Physics, 148(241703):1–9, 2018.
  * [33] Igor Mezić and Andrzej Banaszuk.  Comparison of systems with complex behavior.  Physica D: Nonlinear Phenomena, 197(1-2):101–133, 2004.
  * [34] Yoshihiko Susuki and Igor Mezić.  A prony approximation of koopman mode decomposition.  In Decision and Control (CDC), 2015 IEEE 54th Annual Conference on, pages 7022–7027. IEEE, 2015.
  * [35] Hassan Arbabi and Igor Mezic.  Ergodic theory, dynamic mode decomposition, and computation of spectral properties of the koopman operator.  SIAM Journal on Applied Dynamical Systems, 16(4):2096–2126, 2017\.
  * [36] Kathleen P Champion, Steven L Brunton, and J Nathan Kutz.  Discovery of nonlinear multiscale systems: Sampling strategies and embeddings.  SIAM Journal on Applied Dynamical Systems, 18(1):312–333, 2019\.
  * [37] M. Kamb, E. Kaiser, S. L. Brunton, and J. N. Kutz.  Time-delay observables for Koopman: Theory and applications.  SIAM J. Appl. Dyn. Syst., 19(2):886–917, 2020.
  * [38] Seth M Hirsh, Sara M Ichinaga, Steven L Brunton, J Nathan Kutz, and Bingni W Brunton.  Structured time-delay models for dynamical systems with connections to frenet-serret frame.  arXiv preprint arXiv:2101.08344, 2021.
  * [39] P. J. Schmid.  Dynamic mode decomposition of numerical and experimental data.  Journal of Fluid Mechanics, 656:5–28, August 2010.
  * [40] J Nathan Kutz, Steven L Brunton, Bingni W Brunton, and Joshua L Proctor.  Dynamic mode decomposition: data-driven modeling of complex systems.  SIAM, 2016.
  * [41] Jonathan H Tu, Clarence W Rowley, Dirk M Luchtenburg, Steven L Brunton, and J Nathan Kutz.  On dynamic mode decomposition: Theory and applications.  Journal of Computational Dynamics, 1(2):391, 2014.
  * [42] Bingni W Brunton, Lise A Johnson, Jeffrey G Ojemann, and J Nathan Kutz.  Extracting spatial–temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition.  Journal of neuroscience methods, 258:1–15, 2016.
  * [43] Naoya Takeishi, Yoshinobu Kawahara, and Takehisa Yairi.  Learning koopman invariant subspaces for dynamic mode decomposition.  In Advances in Neural Information Processing Systems, pages 1130–1140, 2017.
  * [44] Bethany Lusch, J Nathan Kutz, and Steven L Brunton.  Deep learning for universal linear embeddings of nonlinear dynamics.  Nature communications, 9(1):1–10, 2018.
  * [45] Enoch Yeung, Soumya Kundu, and Nathan Hodas.  Learning deep neural network representations for koopman operators of nonlinear dynamical systems.  arXiv preprint arXiv:1708.06850, 2017.
  * [46] Andreas Mardt, Luca Pasquali, Hao Wu, and Frank Noé.  VAMPnets: Deep learning of molecular kinetics.  Nature Communications, 9(5), 2018.
  * [47] Kathleen Champion, Bethany Lusch, J Nathan Kutz, and Steven L Brunton.  Data-driven discovery of coordinates and governing equations.  Proceedings of the National Academy of Sciences, 116(45):22445–22451, 2019.
  * [48] Samuel E Otto and Clarence W Rowley.  Linearly-recurrent autoencoder networks for learning dynamics.  SIAM Journal on Applied Dynamical Systems, 18(1):558–593, 2019\.
  * [49] Alec J Linot and Michael D Graham.  Deep learning to discover and predict dynamics on an inertial manifold.  Physical Review E, 101(6):062209, 2020.
  * [50] William Gilpin.  Deep learning of dynamical attractors from time series measurements.  arXiv preprint arXiv:2002.05909, 2020.
  * [51] Kookjin Lee and Kevin T Carlberg.  Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.  Journal of Computational Physics, 404:108973, 2020.
  * [52] Manu Kalia, Steven L Brunton, Hil GE Meijer, Christoph Brune, and J Nathan Kutz.  Learning normal form autoencoders for data-driven discovery of universal, parameter-dependent governing equations.  arXiv preprint arXiv:2106.05102, 2021.
  * [53] Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott.  Model-free prediction of large spatiotemporally chaotic systems from data: a reservoir computing approach.  Physical review letters, 120(2):024102, 2018.
  * [54] Pantelis R Vlachas, Wonmin Byeon, Zhong Y Wan, Themistoklis P Sapsis, and Petros 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(2213):20170844, 2018.
  * [55] Hao Wu, Andreas Mardt, Luca Pasquali, and Frank Noe.  Deep generative markov state models.  32nd Conference on Neural Information Processing Systems (NeurIPS), 2018.
  * [56] Maziar Raissi, Paris Perdikaris, and George 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, 2019.
  * [57] Andreas Mardt, Luca Pasquali, Frank Noé, and Hao Wu.  Deep learning markov and koopman models with physical constraints.  In Mathematical and Scientific Machine Learning, pages 451–475. PMLR, 2020.
  * [58] Julia Ling, Andrew Kurzawski, and Jeremy Templeton.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  Journal of Fluid Mechanics, 807:155–166, 2016.
  * [59] Karthik Duraisamy, Gianluca Iaccarino, and Heng Xiao.  Turbulence modeling in the age of data.  Annual Review of Fluid Mechanics, 51:357–377, 2019.
  * [60] Steven L Brunton, Bernd R Noack, and Petros Koumoutsakos.  Machine learning for fluid mechanics.  Annual Review of Fluid Mechanics, 52:477–508, 2020.
  * [61] Ronald R Coifman, Ioannis G Kevrekidis, Stéphane Lafon, Mauro Maggioni, and Boaz Nadler.  Diffusion maps, reduction coordinates, and low dimensional representation of stochastic systems.  Multiscale Modeling & Simulation, 7(2):842–864, 2008.
  * [62] Or Yair, Ronen Talmon, Ronald R Coifman, and Ioannis G Kevrekidis.  Reconstruction of normal forms by learning informed observation geometries from data.  Proceedings of the National Academy of Sciences, page 201620045, 2017.
  * [63] Hong-guang Ma and Chong-zhao Han.  Selection of embedding dimension and delay time in phase space reconstruction.  Frontiers of Electrical and Electronic Engineering in China, 1(1):111–114, 2006.
  * [64] Matthew B Kennel, Reggie Brown, and Henry DI Abarbanel.  Determining embedding dimension for phase-space reconstruction using a geometrical construction.  Physical review A, 45(6):3403, 1992.
  * [65] J. N. Juang and R. S. Pappa.  An eigensystem realization algorithm for modal parameter identification and model reduction.  Journal of Guidance, Control, and Dynamics, 8(5):620–627, 1985\.
  * [66] J. N. Juang, M. Phan, L. G. Horta, and R. W. Longman.  Identification of observer/Kalman filter Markov parameters: Theory and experiments.  Technical Memorandum 104069, NASA, 1991.
  * [67] M. Phan, J. N. Juang, and R. W. Longman.  Identification of linear-multivariable systems by identification of observers with assigned real eigenvalues.  The Journal of the Astronautical Sciences, 40(2):261–279, 1992\.
  * [68] J. N. Juang.  Applied System Identification.  Prentice Hall PTR, Upper Saddle River, New Jersey, 1994.
  * [69] S. L. Brunton and J. N. Kutz.  Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control.  Cambridge University Press, 2019.
  * [70] Hayden Schaeffer.  Learning partial differential equations via data discovery and sparse optimization.  In Proc. R. Soc. A, volume 473, page 20160446. The Royal Society, 2017.
  * [71] Lorenzo Boninsegna, Feliks Nüske, and Cecilia Clementi.  Sparse learning of stochastic dynamical equations.  The Journal of Chemical Physics, 148(24):241723, 2018.
  * [72] Eurika Kaiser, J Nathan Kutz, and Steven L Brunton.  Sparse identification of nonlinear dynamics for model predictive control in the low-data limit.  Proceedings of the Royal Society of London A, 474(2219), 2018.
  * [73] J.-C. Loiseau and S. L. Brunton.  Constrained sparse Galerkin regression.  Journal of Fluid Mechanics, 838:42–67, 2018.
  * [74] Patrick Gelß, Stefan Klus, Jens Eisert, and Christof Schütte.  Multidimensional approximation of nonlinear dynamical systems.  Journal of Computational and Nonlinear Dynamics, 14(6), 2019.
  * [75] Daniel R Gurevich, Patrick AK Reinbold, and Roman O Grigoriev.  Robust and optimal sparse regression for nonlinear pde models.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(10):103113, 2019.
  * [76] S Beetham and J Capecelatro.  Formulating turbulence closures using sparse regression with embedded form invariance.  Physical Review Fluids, 5(8):084611, 2020.
  * [77] Patrick AK Reinbold, Daniel R Gurevich, and Roman O Grigoriev.  Using noisy or incomplete data to discover models of spatiotemporal dynamics.  Physical Review E, 101(1):010203, 2020.
  * [78] Martin Schmelzer, Richard P Dwight, and Paola Cinnella.  Discovery of algebraic reynolds-stress models using sparse symbolic regression.  Flow, Turbulence and Combustion, 104(2):579–603, 2020.
  * [79] Balachandra Suri, Logan Kageorge, Roman O Grigoriev, and Michael F Schatz.  Capturing turbulent dynamics and statistics in experiments with unstable periodic orbits.  Physical Review Letters, 125(6):064501, 2020.
  * [80] Sarah Beetham, Rodney O Fox, and Jesse Capecelatro.  Sparse identification of multiphase turbulence closures for coupled fluid–particle flows.  Journal of Fluid Mechanics, 914, 2021.
  * [81] Patrick AK Reinbold, Logan M Kageorge, Michael F Schatz, and Roman O Grigoriev.  Robust learning from noisy, incomplete, high-dimensional experimental data via physically constrained symbolic regression.  Nature communications, 12(1):1–8, 2021.
  * [82] Joseph Bakarji and Daniel M Tartakovsky.  Data-driven discovery of coarse-grained equations.  Journal of Computational Physics, 434:110219, 2021.
  * [83] Balázs Csanád Csáji et al.  Approximation with artificial neural networks.  Faculty of Sciences, Etvs Lornd University, Hungary, 24(48):7, 2001\.
  * [84] Harvard Natural Sciences Lecture Demonstrations.  Chaotic waterwheel, 2020.
  * [85] Brian M de Silva, Kathleen Champion, Markus Quade, Jean-Christophe Loiseau, J Nathan Kutz, and Steven L Brunton.  Pysindy: a python package for the sparse identification of nonlinear dynamics from data.  Journal of Open Source Software, 5(49):2104, 2020.
  * [86] Alan A Kaptanoglu, Brian M de Silva, Urban Fasel, Kadierdan Kaheman, Jared L Callaham, Charles B Delahunt, Kathleen Champion, Jean-Christophe Loiseau, J Nathan Kutz, and Steven L Brunton.  Pysindy: A comprehensive python package for robust sparse system identification.  arXiv preprint arXiv:2111.08481, 2021.
  * [87] Urban Fasel, J Nathan Kutz, Bingni W Brunton, and Steven L Brunton.  Ensemble-sindy: Robust sparse model discovery in the low-data, high-noise limit, with active learning and control.  arXiv preprint arXiv:2111.10992, 2021.
  * [88] Andrew M. Fraser.  Information and entropy in strange attractors.  IEEE transactions on Information Theory, 35(2):245–262, 1989.
  * [89] Sebastian Wallot and Dan Mønster.  Calculation of average mutual information (ami) and false-nearest neighbors (fnn) for the estimation of embedding parameters of multidimensional time series in matlab.  Frontiers in psychology, 9:1679, 2018.