# Data-Driven Linearization of Dynamical Systems 111To appear in Nonlinear
Dynamics

George Haller and Bálint Kaszás
Institute for Mechanical Systems
ETH Zürich
Leonhardstrasse 21, 8092 Zürich, Switzerland
Corresponding author. Email<EMAIL_ADDRESS>
Dynamic Mode Decomposition (DMD) and its variants, such as extended DMD
(EDMD), are broadly used to fit simple linear models to dynamical systems
known from observable data. As DMD methods work well in several situations but
perform poorly in others, a clarification of the assumptions under which DMD
is applicable is desirable. Upon closer inspection, existing interpretations
of DMD methods based on the Koopman operator are not quite satisfactory: they
justify DMD under assumptions that hold only with probability zero for generic
observables. Here, we give a justification for DMD as a local, leading-order
reduced model for the dominant system dynamics under conditions that hold with
probability one for generic observables and non-degenerate observational data.
We achieve this for autonomous and for periodically forced systems of finite
or infinite dimensions by constructing linearizing transformations for their
dominant dynamics within attracting slow spectral submanifolds (SSMs). Our
arguments also lead to a new algorithm, data-driven linearization (DDL), which
is a higher-order, systematic linearization of the observable dynamics within
slow SSMs. We show by examples how DDL outperforms DMD and EDMD on numerical
and experimental data.

##  1 Introduction

In recent years, there has been an overwhelming interest in devising linear
models for dynamical systems from experimental or numerical data (see the
recent review by Schmid [51]). This trend was largely started by the Dynamic
Mode Decomposition (DMD), put forward in seminal work by Schmid [50]. The
original exposition of the method has been streamlined by various authors,
most notably by Rowley et al. [49] and Kutz et al. [32].

To describe DMD, we consider an autonomous dynamical system

$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x}),\qquad\mathbf{x}\in\mathcal{\mathbb{R}}^{n},\qquad\mathbf{f}\in
C^{1}\left(\mathcal{\mathbb{R}}^{n}\right),$ (1)

for some $n\in\mathbb{N}^{+}$. Trajectories
$\left\\{\mathbf{x}(t;\mathbf{x}_{0})\right\\}_{t\in\mathbb{R}}$ of this
system evolve from initial conditions $\mathbf{x}_{0}$. The flow map
$\mathbf{F}^{t}\colon\mathcal{\mathbb{R}}^{n}\to\mathcal{\mathbb{R}}^{n}$ is
defined as the mapping taking the initial trajectory positions at time
$t_{0}=0$ to current ones at time $t$, i.e.,

$\mathbf{F}^{t}(\mathbf{x}_{0})=\mathbf{x}(t;\mathbf{x}_{0}).$ (2)

As observations of the full state space variable $\mathbf{x}$ of system (1)
are often not available, one may try to explore the dynamical system (1) by
observing $d$ smooth scalar functions
$\phi_{1}(x),\ldots,\text{$\phi$}_{d}(x)$ along trajectories of the system. We
order these scalar observables into the observable vector

$\boldsymbol{\phi}(x)=\left(\begin{array}[]{c}\phi_{1}(\mathbf{x})\\\
\vdots\\\ \phi_{d}(\mathbf{x})\end{array}\right)\in
C^{1}\left(\mathcal{\mathbb{R}}^{n}\right)$ (3)

The basic idea of DMD is to approximate the observed evolution of
$\boldsymbol{\phi}\left(\mathbf{F}^{t}(\mathbf{x}_{0})\right)$ of the
dynamical system with the closest fitting autonomous linear dynamical system

$\dot{\boldsymbol{\phi}}=\mathbf{L}\boldsymbol{\phi},\qquad\mathbf{L}\in\mathbb{R}^{d\times
d},$ (4)

based on available trajectory observations.

This is a challenging objective for multiple reasons. First, the original
dynamical system (1) is generally nonlinear whose dynamics cannot be well
approximated by a single linear system on a sizable open domain. For instance,
one may have several isolated, coexisting attracting or repelling stationary
states (such as periodic orbits, limit cycles or quasiperiodic tori), which
linear systems cannot have. Second, it is unclear why the dynamics of $d$
observables should be governed by a self-contained autonomous dynamical system
induced by the original system (1), whose dimension is $n$. Third, the result
of fitting system (4) to observable data will clearly depend on the initial
conditions used, the number and the functional form of the observables chosen,
as well as on the objective function used in minimizing the fit.

Despite these challenges, we may proceed to find an appropriately defined
closest linear system (4) based on available observable data. We assume that
for some fixed time step $\Delta t$, discrete observations of $m$ initial
conditions, $\mathbf{x}^{1}(t_{0}),\ldots,\mathbf{x}^{m}(t_{0})$, and their
images $\mathbf{F}^{t_{0}+\Delta
t}(\mathbf{x}^{1}(t_{0})),\ldots,\mathbf{F}^{t_{0}+\Delta
t}\left(\mathbf{x}^{m}(t_{0})\right)$, under the sampled flow map
$\mathbf{F}^{t_{0}+\Delta t}$ are available in the data matrices

$\boldsymbol{\Phi}=\left[\boldsymbol{\phi}\left(\mathbf{x}^{1}(t_{0})\right),\ldots,\boldsymbol{\phi}\left(\mathbf{x}^{m}(t_{0})\right)\right],\qquad\hat{\boldsymbol{\Phi}}=\left[\boldsymbol{\phi}\left(\mathbf{F}^{t_{0}+\Delta
t}(\mathbf{x}^{1}(t_{0}))\right),\ldots,\boldsymbol{\phi}\left(\mathbf{F}^{t_{0}+\Delta
t}\left(\mathbf{x}^{m}(t_{0})\right)\right)\right],$ (5)

respectively. We seek the best fitting linear system of the form (4) for which

$\hat{\boldsymbol{\Phi}}\approx\boldsymbol{\mathcal{D}}\boldsymbol{\Phi},\quad\boldsymbol{\mathcal{D}}=e^{\mathbf{L}\Delta
t},$ (6)

holds. The eigenvalues of such a $\boldsymbol{\mathcal{D}}$ are usually called
DMD eigenvalues, and their corresponding eigenvectors are called the DMD
modes.

Various norms can be chosen with respect to which the difference of
$\hat{\boldsymbol{\Phi}}$ and
$\mathcal{\boldsymbol{\mathcal{D}}}\boldsymbol{\Phi}$ is to be minimized. The
most straightforward choice is the Euclidean matrix norm
$\left|\,\cdot\,\right|$, which leads to the minimization principle

$\boldsymbol{\mathcal{D}}^{*}=\underset{\boldsymbol{\mathcal{D}}\in\mathbb{R}^{d\times
d}}{\mathrm{argmin}\left|\hat{\boldsymbol{\Phi}}-\boldsymbol{\mathcal{D}}\boldsymbol{\Phi}\right|^{2}}.$
(7)

An explicit solution to this problem is given by

$\boldsymbol{\mathcal{D}}=\left(\hat{\boldsymbol{\Phi}}\boldsymbol{\Phi}^{\mathrm{T}}\right)\left(\hat{\boldsymbol{\Phi}}\boldsymbol{\Phi}^{\mathrm{T}}\right)^{\dagger},$
(8)

with the dagger referring to the pseudo-inverse of a matrix (see, e.g., Kutz
et al. [32] for details). We note that the original formulation of Schmid [50]
is for discrete dynamical processes and assumes observations of a single
trajectory (see also Rowley et al. [49]).

Among several later variants of DMD surveyed by Schmid [51], the most broadly
used one is the Extended Dynamic Mode Decomposition (EDMD) of Williams et al.
[58]. This procedure seeks the best-fitting linear dynamics for an a priori
unknown set of functions $\mathbf{K}(\boldsymbol{\phi}(\mathbf{x}))$ of
$\boldsymbol{\phi}(\mathbf{x})$, rather than for
$\boldsymbol{\phi}(\mathbf{x})$ itself. In practice, one often chooses
$\mathbf{K}$ as an $N(d,k)$-dimensional vector of $d$-variate scalar monomials
of order $k$ or less, where $N(d,k)=\left(\begin{array}[]{c}d+k\\\
k\end{array}\right)$ is the total number of all such monomials. The underlying
assumption of EDMD is that a self-contained linear dynamical system of the
form

$\frac{d}{dt}\mathbf{K}\left(\boldsymbol{\phi}(\mathbf{F}^{t}(\mathbf{x}_{0}))\right)=\mathbf{L}\mathbf{K}\left(\boldsymbol{\phi}(\mathbf{F}^{t}(\mathbf{x}_{0}))\right)$
(9)

can be obtained on the feature space $\mathbb{R}^{N(d,k)}$ by optimally
selecting $\mathbf{L}\in\mathbb{R}^{N(d,k)\times N(d,k)}$. For physical
systems, the $N(d,k)$-dimensional ODE in eq. (9) defined on the feature space
$\mathbb{R}^{N(d,k)}$ can be substantially higher-dimensional than the
$d$-dimensional ODE (4). In fact, $N(d,k)$ may be substantially higher than
the dimension $n$ of the phase space $\mathbb{R}^{n}$ of the original
nonlinear system (1).

Once the function library used in EDMD is fixed, one again seeks to choose
$\mathbf{L}$ so that

$\mathbf{K}(\hat{\boldsymbol{\Phi}})\approx\boldsymbol{\mathcal{D}}\mathbf{K}(\boldsymbol{\Phi}^{\mathrm{}}),\quad\boldsymbol{\mathcal{D}}=e^{\mathbf{L}\Delta
t}.$

This again leads to a linear optimization problem that can be solved using
linear algebra tools. For higher-dimensional systems, a kernel-based version
of EDMD was developed by Williams et al. [59]. This method computes inner
products necessary for EDMD implicitly, without requiring an explicit
representation of (polynomial) basis functions in the space of observables. As
a result, kernel-based EDMD operates at computational costs comparable to
those of the original DMD.

##  2 Prior justifications for DMD methods

Available justifications for DMD (see Rowley et al. [49]) and EDMD (see
Williams et al. [58]) are based on the Koopman operator, whose basics we
review in Appendix A.1 for completeness. The argument starts with the
observation that special observables falling in invariant subspaces of this
operator in the space of all observables obey linear dynamics. Consequently,
DMD should recover the Koopman operator restricted to this subspace if the
observables are taken from such a subspace.

In this sense, DMD is viewed as an approximate, continuous immersion of a
nonlinear system into an infinite dimensional linear dynamical system. While
such an immersion is not possible for typical nonlinear systems with multiple
limit sets (see Liu et al. [38, 39]), one still hopes that this approximate
immersion is attainable via DMD or EDMD for nonlinear systems with a single
attracting steady state that satisfies appropriate nondegeneracy conditions
(see Kvalheim and Arathoon [33]). In that case, unlike classic local
linearization near fixed points, the linearization via DMD or EDMD is argued
to be non-local, as it covers the full domain of definition of Koopman
eigenfunctions spanning the underlying Koopman-invariant subspace.

However, Koopman eigenfunctions, whose existence, domain of definition and
exact form are a priori unknown for general systems, are notoriously
difficult–if not impossible–to determine accurately from data. More
importantly, even if Koopman-invariant subspaces of the observable space were
known, any countable set of generically chosen observables would lie outside
those subspaces with probability one. As a consequence, DMD eigenvectors
(which are generally argued to be approximations of Koopman eigenfunctions and
can be used to compute Koopman modes222Assuming that each coordinate component
$x_{j}$ of the full phase space vector $\mathbf{x}$ of system (1) falls in a
span of the same set of Koopman eigenfunctions
$\left\\{\phi_{1}(\mathbf{x}),\ldots,\phi_{N}(\mathbf{x})\right\\}$, one
defines the $j^{th}$ Koopman mode associated with $x_{j}$ as the vector
$\mathbf{v}_{j}=\left(v_{1j},\dots,v_{Nj}\right)^{\mathrm{T}}$ for which
$x_{j}=\sum_{i=1}^{N}v_{ij}\phi_{i}(x)$ (see, e.g., Williams et al. [58]).)
would also lie outside Koopman-invariant subspaces, given that such
eigenvectors are just linear combinations of the available observables.
Consequently, practically observed data sets would fall under the realm of
Koopman-based explanation for DMD with probability zero. This is equally true
for EDMD, whose flexibility in choosing the function set
$\mathbf{K}(\boldsymbol{\phi}(\mathbf{x}))$ of observables also introduces
further user-dependent heuristics beyond the dimension $d$ of the DMD.

One may still hope that by enlarging the dimension $d$ of observables in DMD
and enlarging the function library $\mathbf{K}(\boldsymbol{\phi}(\mathbf{x}))$
in EDMD, the optimization involved in these methods brings DMD and EDMD
eigenvectors closer and closer to Koopman eigenfunctions. The required
enlargements, however, may mean hundreds or thousand of dimensions even for
dynamical system governed by simple, low-dimensional ODEs (Williams et al.
[59]). These enlargements succeed in fitting linear systems closely to sets of
observer trajectories, but they also unavoidably lead to overfits that give
unreliable predictions for initial conditions not used in the training of DMD
or EDMD. Indeed, the resulting large linear systems can perform substantially
worse in prediction than much lower dimensional linear or nonlinear models
obtained from other data-driven techniques (see, e.g., Alora et al. [1]).

Similar issues arise in justifying the kernel-based EDMD of Williams et al.
[59] based on the Koopman operator. Additionally, the choice of the kernel
function that represents the inner product of the now implicitly defined
polynomial basis functions remains heuristic and problem-dependent. Again, the
accuracy of the procedure is not guaranteed, as available observer data is
generically not in a Koopman eigenspace. As Williams et al. [59] write: “Like
most existing data-driven methods, there is no guarantee that the kernel
approach will produce accurate approximations of even the leading eigenvalues
and modes, but it often appears to produce useful sets of modes in practice if
the kernel and truncation level of the pseudoinverse are chosen properly.”

Finally, a lesser known limitation of the Koopman-based approach to DMD is the
limited domain in the phase space over which Koopman eigenfunctions (and hence
their corresponding invariant subspaces) are defined in the observable space.
Specifically, all $C^{\infty}$ Koopman eigenfunctions necessarily blow up near
basin boundaries of attracting and repelling fixed points, periodic orbits and
more general invariant sets (see Proposition 1 of our Appendix A.4 for a
precise statement and proof).

Expansions of observables in terms of such blowing-up eigenfunctions have even
smaller domains of convergence, as was shown explicitly in a simple example by
Page and Kerswell [44]. This is a fundamental obstruction to the often
envisioned concept of global linear models built of different Koopman
eigenfunctions over multiple domains of attraction (see, e.g., Williams et al.
[58], p. 1309). While it is broadly known that such models would be
discontinuous along basin boundaries (Liu et al. [38, 39], Kvalheim and
Arathoon [33]), it is rarely noted that such models would also generally blow
up at those boundaries and hence would become unmanageable even before
reaching the boundaries.

For these reasons, an alternative mathematical foundation for DMD is
desirable. Ideally, such an approach should be defined on an equal or lower
dimensional space, rather than on higher or even infinite-dimensional spaces,
as suggested by the Koopman-based view on DMD. This should help in avoiding
overfitting and computational difficulties. Additionally, an ideal treatment
of DMD should also provide specific non-degeneracy conditions on the
underlying dynamical system, on the available observables, and on the specific
data to be used in the DMD procedure.

In this paper, we develop a treatment of DMD that satisfies these
requirements. This enables to derive conditions for DMD to approximate the
dominant linearized observable dynamics near hyperbolic fixed points and
periodic orbits of finite- and infinite-dimensional dynamical systems.

Our approach to DMD also leads to a refinement of DMD which we call data-
driven linearization (or DDL, for short). DDL effectively carries out exact
local linearization via nonlinear coordinate changes on a lower-dimensional
attracting invariant manifold (spectral submanifold) of the dynamical system.
We illustrate the increased accuracy and domain of validity of DDL models
relative to those obtained from DMD and EDMD on examples of autonomous and
forced dynamical systems.

##  3 A simple justification for the DMD algorithm

Here we give an alternative interpretation of DMD and EDMD as approximate
models for a dynamical system known through a set of observables. The main
idea (to be made more precise shortly) is to view DMD executed on $d$
observables $\phi_{1},\ldots,\phi_{d}$ as a model reduction tool that captures
the leading-order dynamics of $n\geq d$ phase space variables along a
$d$-dimensional slow manifold in terms of $\phi_{1},\ldots,\phi_{d}$.

Such manifolds arise as slow spectral submanifolds (SSMs) under weak non-
degeneracy assumptions on the linearized spectrum at stable hyperbolic fixed
points of the $n$-dimensional dynamical system (see Cabré et al. [10], Haller
and Ponsioen [21], Cenedese et al. [12]). Specifically, a slow SSM
$\mathcal{W}\left(E\right)$ is tangent to the real eigenspace $E$ spanned by
the $d$ slowest decaying linearized modes at the fixed point. If $m$ sample
trajectories $\left\\{\mathbf{x}^{j}(t)\right\\}_{j=1}^{m}$ are released from
a set of initial conditions $\left\\{\mathbf{x}^{j}(0)\right\\}_{j=1}^{m}$at
time $t=0$, then due to their fast decay along the remaining fast spectral
subspace $F$ , the $\mathbf{x}^{j}(t)$ trajectories will become exponentially
close to $\mathcal{W}\left(E\right)$ by some time $t=t_{0}\geq 0$, and closely
synchronize with its internal dynamics, as seen in Fig. 1.

Figure 1: The geometric meaning of DMD performed on $d$ observables,
$\phi_{1}(\mathbf{x}),\ldots,\phi_{d}(\mathbf{x})$, after initial fast
transients die out at an exponential rate in the data. DMD then identifies the
leading-order (linear) dynamics on a $d$-dimensional attracting spectral
submanifold (SSM) $\mathcal{W}\left(E\right)$ tangent to the $d$-dimensional
slow spectral subspace $E$. These linearized dynamics can be expressed in
terms of the SSM-restricted observables
$\boldsymbol{\varphi}=\boldsymbol{\phi}|_{\mathcal{W}\left(E\right)}$. Also
shown is the spectral subspace $F$ of the faster decaying modes and its
associated nonlinear continuation, the fast spectral subspace
$\mathcal{W}\left(F\right)$.

If $\mathcal{W}\left(E\right)$ admits a non-degenerate parametrization with
respect to the observables $\phi_{1},\ldots,\phi_{d}$, then one can pass to
these observables as new coordinates in which to approximate the leading
order, linearized dynamics inside $\mathcal{W}\left(E\right)$. Specifically,
DMD executed over time the interval $\left[t_{0},t_{0}+\Delta t\right]$
provides the closest linear fit to the reduced dynamics of
$\mathcal{W}\left(E\right)$ from the available observable histories
$\phi_{1}\left(\mathbf{x}^{j}(t)\right),\ldots,\phi_{d}\left(\mathbf{x}^{j}(t)\right)$
for $t\in\left[t_{0},t_{0}+\Delta t\right]$, as illustrated in Fig. 1. This
fit is a close representation of the actual linearized dynamics on the SSM
$\mathcal{W}\left(E\right)$ if the trajectory data
$\left\\{\mathbf{x}^{j}(t)\right\\}_{j=1}^{m}$ is sufficiently diverse for
$t\geq t_{0}$. The resulting DMD model with coefficient matrix
$\mathcal{\boldsymbol{\mathcal{D}}}$ will then be smoothly conjugate to the
linearized reduced dynamics on $\mathcal{W}\left(E\right)$ with an error of
the order of the distance of$\left\\{\mathbf{x}^{j}(t)\right\\}_{j=1}^{m}$
from $\mathcal{W}\left(E\right)$ between the times $t_{0}$ and $t_{0}+\Delta
t$.

The slow SSM $\mathcal{W}\left(E\right)$ may also contain unstable modes in
applications. Similar results hold for that case as well, provided that we
select a small enough $\Delta t$, ensuring that the trajectories
$\left\\{\mathbf{x}^{j}(0)\right\\}_{j=1}^{m}$are not ejected from the
vicinity of the fixed point origin for $t\in\left[t_{0},t_{0}+\Delta
t\right]$. As DMD will show sensitivity with respect to the choice of $t_{0}$
and $\Delta t$ in this case, we will not discuss the justification of DMD near
unstable fixed points beyond Remarks 1 and 7.

In the following sections, we make this basic idea more precise both for
finite and infinite-dimensional dynamical systems. This approach also reveals
explicit, previously undocumented non-degeneracy conditions on the underlying
dynamical system, on the available observable functions and on the specific
data set used, under which DMD should give meaningful results.

###  3.1 Justification of DMD for continuous dynamical systems

We start by assuming that the observed dynamics take place in a domain
containing a fixed point, which is assumed to be at the origin, without loss
of generality, i.e.,

$\mathbf{f}(\mathbf{0})=\mathbf{0}.$ (10)

We can then rewrite the dynamical system (1) in the more specific form

$\dot{\mathbf{x}}=\mathbf{Ax}+\tilde{\mathbf{f}}(\mathbf{x}),\qquad\mathbf{x}\in\mathcal{\mathbb{R}}^{n},\qquad\mathbf{A}=D\mathbf{f}(\mathbf{0}),\qquad\tilde{\mathbf{f}}(\mathbf{x})=o\left(\left|\mathbf{x}\right|\right),$
(11)

where the classic $o\left(\left|\mathbf{x}\right|\right)$ notation refers to
the fact that
$\lim_{\mathbf{x}\to\mathbf{0}}\left[\left|\tilde{\mathbf{f}}(\mathbf{x})\right|/\left|\mathbf{x}\right|\right]=0$.

Most expositions of DMD methods and their variants do not state assumption
(10) explicitly and hence may appear less restrictive than our treatment here.
However, all of them implicitly assume the existence of such a fixed point, as
all of them end up returning homogeneous linear ODEs or mappings with a fixed
point at the origin. Indeed, all known applications of these methods that
produce reasonable accuracy target the dynamics of ODEs or discrete maps near
their fixed points.

Assumption (10) can be replaced with the existence of a limit cycle in the
original system (1), in which case the first return map (or Poincaré map)
defined near the limit cycle will have a fixed point. We give a separate
treatment on justifying DMD as a linearization for such a Poincaré map in
Section 3.2.

We do not advocate, however, the often used procedure of applying DMD to fit a
linear system to the flow map (rather than the Poincaré map) near a stable
limit cycle. Such a fit only produces the desired limiting periodic behavior
if one or more of the DMD eigenvalues are artificially constrained to be on
the complex unit circle by the user of DMD. This renders the DMD model both
structurally unstable and conceptually inaccurate for prediction. Indeed, the
model will approximate the originally observed limit cycle and convergence to
it only within a measure zero, cylindrical set of its phase space. Outside
this set, all trajectories of the DMD model will converge to some other member
of the infinite family of periodic orbits or invariant tori within the center
subspace corresponding to the unitary eigenvalues. These periodic orbits or
tori have a continuous range of locations and amplitudes, and hence represent
spurious asymptotic behaviors that are not seen in the original dynamical
system (1).

In the special case of $d=n$ and for the special observable
$\boldsymbol{\phi}(x)=\mathbf{x}$, the near-linear form of eq. (11) motivates
the DMD procedure because a linear approximation to the system near
$\mathbf{x}=\mathbf{0}$ seems feasible. It is a priori unclear, however, to
what extent the nonlinearities distort the linear dynamics and how DMD would
account for that. Additionally, in a data-driven analysis, choosing the full
phase space variable $\mathbf{x}$ as the observable
$\boldsymbol{\phi}(\mathbf{x})$ is generally unrealistic. For these reasons, a
mathematical justification of DMD requires further assumptions, as we discuss
next.

Let $\lambda_{1},\ldots,\lambda_{n}\in\mathbb{C}$ denote the eigenvalues of
$A$ and let $\mathbf{e}_{1},\ldots,\mathbf{e}_{n}\in\mathbb{C}^{n}$ denote the
corresponding generalized eigenvectors. We assume that at least one of the
modes of the linearized system at $\mathbf{x}=\mathbf{0}$ decays exponentially
and there is at least one other mode that decays slower or even grows. More
specifically, for some positive integer $d<n$, we assume that the spectrum of
$A$ can be partitioned as

$\mathrm{Re}\text{$\lambda_{n}$}\leq\ldots\leq\mathrm{Re}\text{$\lambda_{d+1}$
}<\mathrm{Re}\text{$\lambda_{d}$
}\leq\ldots\leq\mathrm{Re}\text{$\lambda_{1}<0$. }$ (12)

This guarantees the existence of a $d$-dimensional, normally attracting _slow
spectral subspace_

$E=\mathrm{span\,\left\\{\mathrm{Re}\,\text{$\mathbf{e}_{1}$,\,}\mathrm{Im\,}\text{$\mathbf{e}_{1}$,}\ldots,\mathrm{Re\,}\text{$\mathbf{e}_{d}$,\,}\mathrm{Im}\,\text{$\mathbf{e}_{d}$}\right\\}}$
(13)

for the linearized dynamics, with linear decay rate towards $E$ strictly
dominating all decay rates inside $E$. We also define the (real) spectral
subspace of faster decaying linear modes:

$F=\mathrm{span\,\left\\{\mathrm{Re}\,\text{$\mathbf{e}_{d+1}$,\,}\mathrm{Im\,}\text{$\mathbf{e}_{d+1}$,}\ldots,\mathrm{Re\,}\text{$\mathbf{e}_{n}$,\,}\mathrm{Im}\,\text{$\mathbf{e}_{n}$}\right\\}}.$
(14)

We will also use matrices containing the left and right eigenvectors of the
operator $\mathbf{A}$ and of its restrictions, $\mathbf{A}|_{E}$ and
$\mathbf{A}|_{F}$, to its spectral subspaces $E$ and $F$, respectively.
Specifically, we let

$\displaystyle\mathbf{T}$
$\displaystyle=\left[\mathbf{T}_{E},\mathbf{T}_{F}\right],\qquad\mathbf{P}=\left[\begin{array}[]{c}\mathbf{P}_{E}\\\
\mathbf{P}_{F}\end{array}\right],$ (17)

where the columns of $\mathbf{T}_{E}\in\mathbb{R}^{n\times d}$ are the real
and imaginary parts of the generalized right eigenvectors of $\mathbf{A}|_{E}$
and the columns of $\mathbf{T}_{F}\in\mathbb{R}^{n\times(n-d)}$ are defined
analogously for $\mathbf{A}|_{F}.$ Similarly, the rows of
$\mathbf{P}_{E}\in\mathbb{R}^{d\times n}$ are the real and imaginary parts of
the generalized left eigenvectors of $\mathbf{A}|_{E}$ and the rows of
$\mathbf{P}_{F}\in\mathbb{R}^{n\times(n-d)}$ are defined analogously
for$\mathbf{A}|_{F}.$ Under assumption (12), $F$ is always a fast spectral
subspace, containing all trajectories of the linearized system that decay
faster to the origin as any trajectory in $E.$

We will use the notation

$\mathbf{X}=\left[\mathbf{x}^{1}(t_{0}),\ldots,\mathbf{x}^{m}(t_{0})\right],\quad\hat{\mathbf{X}}=\left[\mathbf{F}^{t_{0}+\Delta
t}(\mathbf{x}^{1}(t_{0})),\ldots,\mathbf{F}^{t_{0}+\Delta
t}\left(\mathbf{x}^{m}(t_{0})\right)\right]$ (18)

for trajectory data in the underlying dynamical system ((1)) on which the
observable data matrices $\boldsymbol{\Phi}$ and $\hat{\boldsymbol{\Phi}}$
defined in (5) are defined. In truly data-driven applications, the matrices
$\mathbf{X}$ and $\hat{\mathbf{X}}$ are not known. We will nevertheless use
them to make precise statements about a required dominance of the slow linear
modes of $E$ in the available data. Such a dominance will arise for generic
initial conditions if one selects the initial conditions
$\mathbf{x}^{1}(t_{0}),\ldots,\mathbf{x}^{m}(t_{0})$ after initial fast
transients along $F$ have died out. This can be practically achieved by
initializing $\mathbf{x}^{1}(t_{0}),\ldots,\mathbf{x}^{m}(t_{0})$ after a
linear spectral analysis of the observable data matrix $\boldsymbol{\Phi}$
only returns only $d$ dominant frequencies.

We now state a theorem that provides a general justification for the DMD
procedure under explicit nondegeneracy conditions and with specific error
bounds. Specifically, we give a minimal set of conditions under which DMD can
be justified as an approximate, leading-order, $d$-dimensional reduced-order
model for an nonlinear system of dimension $n\geq d$ near its fixed point.
Based on relevance for applications, we only state Theorem 1 for stable
hyperbolic fixed points, but discuss subsequently in Remark 1 its extension to
unstable fixed points.

######  Theorem 1.

_[__Justification of DMD for ODEs with stable hyperbolic fixed points]__Assume
that_

  * _(A1)_

_The $\mathbf{x}=\mathbf{0}$ is a stable hyperbolic fixed point of system
_(11) with a spectral gap, i.e., the spectrum
$\mathrm{Spect}\left[\mathbf{A}\right]$ satisfies eq. (12).

  * _(A2)_

_$\mathbf{f}\in C^{2}$ in a neighborhood of the origin._

  * _(A3)_

_For some integer $d\in\left[1,n\right]$, a $d$-dimensional observable
function $\boldsymbol{\phi}\in C^{2}$ and the $d$-dimensional slow spectral
subspace $E$ of the hyperbolic fixed point $\mathbf{x}=\mathbf{0}$ of system
(_11 _) satisfy the non-degeneracy condition_

$\mathrm{rank\,}\left[D\boldsymbol{\phi}\left(\mathbf{0}\right)|_{E}\right]=d.$
(19)

  * _(A4)_

_The data matrices $\boldsymbol{\Phi}$ and $\hat{\boldsymbol{\Phi}}$ are non-
degenerate and the initial conditions in $\mathbf{X}$ and $\hat{\mathbf{X}}$
have been selected after fast transients from the modes outside $E$ have
largely died out, i.e.,_

$\mathrm{rank_{\mathrm{row}}\,}\boldsymbol{\Phi}=d,\qquad\left|\mathbf{P}_{F}^{\mathrm{}}\mathbf{X}\right|,\left|\mathbf{P}_{F}^{\mathrm{}}\hat{\mathbf{X}}\right|\leq\left|\mathbf{P}_{E}^{\mathrm{}}\mathbf{X}\right|^{1+\beta},$
(20)

_for some $\beta\in(0,1]$._

_Then the DMD computed from $\boldsymbol{\Phi}$ and $\hat{\boldsymbol{\Phi}}$
yields a matrix $\boldsymbol{\mathcal{D}}$ that is locally topologically
conjugate with order
$\mathcal{O}\left(\left|\mathbf{P}_{E}^{\mathrm{}}\mathbf{X}\right|^{\beta}\right)$
error to the linearized dynamics on a $d$-dimensional, slow attracting
spectral submanifold_ $\mathcal{W}(E)\in C^{1}$ tangent to $E$ at
$\mathbf{x}=\mathbf{0}$_. Specifically, we have_

$\boldsymbol{\mathcal{D}}=D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}e^{\mathbf{\Lambda}_{E}\Delta
t}\left(D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\right)^{-1}+\mathcal{O}\left(\left|\mathbf{P}_{E}^{\mathrm{}}\mathbf{X}\right|^{\beta}\right).$
(21)

###### Proof.

Under assumptions (A1) and (A2), any trajectory in a neighborhood of the
origin in the nonlinear system (11) converges at an exponential rate
$e^{\mathrm{Re}\text{$\lambda_{d+1}$ }t}$ to a $d$-dimensional attracting
spectral submanifold $\mathcal{W}(E)$ tangent to a $d$-dimensional attracting
slow spectral subspace $E$ of the linearized system at the origin. This
follows from the $C^{1}$ linearization theorem of Hartman [23], which is
applicable to $C^{2}$ dynamical systems with a stable hyperbolic fixed point.
Under assumption (A3), the $d$-dimensional observable function
$\boldsymbol{\phi}(\mathbf{x})$ restricted to $\mathcal{W}(E)$ can be used to
parametrize $\mathcal{W}(E)$ near the origin, and hence a $d$-dimensional,
self-contained nonlinear dynamical system can be written down for the
restricted observable
$\boldsymbol{\varphi}=\boldsymbol{\phi}|_{\mathcal{W}(E)}$ along
$\mathcal{W}(E)$. Under the first assumption in (A4), the available
observational data matrices $\boldsymbol{\Phi}$ and $\hat{\boldsymbol{\Phi}}$
are rich enough to characterize the reduced dynamics on $\mathcal{W}(E)$.
Under the second assumption in (A4), transients from the faster modes outside
$E$ have largely died out before the selection of the initial conditions in
$\mathbf{X}$, so that the linear part of the dynamics on $\mathcal{W}(E)$ can
be approximately inferred from $\boldsymbol{\Phi}$ and
$\hat{\boldsymbol{\Phi}}$. In that case, up to an error proportional to the
distance of the training data from $\mathcal{W}(E)$, the matrix
$\boldsymbol{\mathcal{D}}\in\mathbb{R}^{d\times d}$ returned by DMD is similar
to the time-$\Delta t$ flow map of the linearized flow of the underlying
dynamical system restricted to $\mathcal{W}(E)$. This linearized flow then
acts as a local reduced-order model with which nearby trajectory observations
synchronize exponentially fast in the observable space. We give a more
detailed proof of the theorem in Appendix (B). ∎

######  Remark 1.

In Theorem 1, we can replace assumption (A1) with

$\mathrm{Re}\text{$\lambda_{n}$}\leq\ldots\leq\mathrm{Re}\text{$\lambda_{d+1}$
}<0,\qquad\mathrm{Re}\text{$\lambda_{d+1}$ }<\mathrm{Re}\text{$\lambda_{d}$
}\leq\ldots\leq\mathrm{Re}\text{$\lambda_{1}$,
}\quad\mathrm{Re}\text{$\lambda_{j}\neq 0,\quad j=1,\ldots,n.$ }$ (22)

This means that the $\mathbf{x}=\mathbf{0}$ fixed point is only assumed
hyperbolic with a spectral gap and $\mathbf{A}$ has an attracting
$d$-dimensional spectral subspace $E$ that possibly contains some
instabilities, i.e., eigenvalues with positive real parts. Then the statements
of Theorem 1 still hold, but $\mathcal{W}(E)$ will be only be guaranteed
$C^{1}$ at $\mathbf{x}=\mathbf{0}$ and Hölder-continuous at other points near
the fixed point. This follows by replacing the linearization theorem of
Hartman [23] with that of van Strien [56], which still enables us to use eq.
(87) in the proof. Therefore, slow subspaces $E$ containing a mixture of
stable and unstable modes can also be allowed, as long as $F$ contains only
fast modes consistent with the splitting assume in eq. (12). In that case,
however, the time $t_{0}+\Delta t$ must be chosen carefully to ensure that
$\left|\mathbf{P}_{F}^{\mathrm{}}\hat{\mathbf{X}}\right|\leq\left|\mathbf{P}_{E}^{\mathrm{}}\mathbf{X}\right|^{1+\beta}$
still holds, i.e., the data used in DMD still samples a neighborhood of the
origin.

######  Remark 2.

In related work, Bollt et al. [7] construct the transformation relating a pair
of conjugate dynamical systems based on a limited set of matching Koopman
eigenfunctions, which are either known explicitly or constructed from EDMD
with dictionary learning (EDMD-DL; see Li et al. [36]). In principle, this
could be used to construct linearizing transformation as well. However, even
when the eigenfunctions are approximated from data, the approach assumes that
the linearized system, as well as a linearized trajectory and its preimage
under the linearization, are available. As these assumptions are not satisfied
in practice, only very simple and low-dimensional analytic examples are
treated by Bollt et al. [7].

In Appendix (B), Remarks 4-5 summarize technical points on the application and
possible further extensions of Theorem 1. In practice, Theorem 1 provides
previously unspecified non-degeneracy conditions on the linear part of the
dynamical system to be analyzed via DMD (assumption (A1)), on the regularity
of the nonlinear part of the system (assumption (A2)), on the type of
observables available for the analysis (assumption (A3)) and on the specific
observable data used in the analysis (assumption (A4)). The latter assumption
requires that there have to be at least as many independent observations in
time as observables. This specifically excludes the popular use of tall
$\boldsymbol{\Phi}$ observable data matrices which provide more free
parameters to pattern-match observational data but will also lead to an
overfit that diminishes the predictive power of the DMD model on initial
conditions not used in its training.

To illustrate these points, we demonstrate the necessity of assumptions
(A2)-(A4) of Theorem 1 in Appendix C on simple examples.

###  3.2 Justification of DMD for discrete and for time-periodic continuous
dynamical systems

The linearization results we have applied to deduce Theorem 1 are equally
valid for discrete dynamical systems defined by iterated mappings. Such
mappings are of the form

$\mathbf{x}_{n+1}=\mathbf{f}(\mathbf{x}_{n})=\mathbf{A}\mathbf{x}_{n}+\tilde{\mathbf{f}}(\mathbf{x}_{n}),\qquad\mathbf{x}_{j}\in\mathbb{R}^{n},\qquad\mathbf{A}\in\mathbb{R}^{n\times
n},\qquad\tilde{\mathbf{f}}(\mathbf{x})=o\left(\left|\mathbf{x}\right|\right).$
(23)

We will use a similar ordering for the eigenvalues of $\mathbf{A}$ as in the
continuous time case:

$\left|\text{$\lambda_{n}$}\right|\leq\ldots\leq\left|\text{$\lambda_{d+1}$
}\right|<\left|\text{$\lambda_{d}$
}\right|\leq\ldots\leq\left|\text{$\lambda_{1}$ }\right|<1.$ (24)

As in the continuous time case, we will use the observable data matrices

$\boldsymbol{\Phi}=\boldsymbol{\phi}\left(\mathbf{X}\right),\qquad\hat{\boldsymbol{\Phi}}=\boldsymbol{\phi}\left(\mathbf{f}\left(\mathbf{X}\right)\right),$
(25)

with the initial conditions for the map $\mathbf{f}$ stored in $\mathbf{X}$.

With these ingredients, we need only minor modifications in the assumptions of
the theorems that account for the usual differences between the spectrum of an
ODE and a map.

######  Theorem 2.

_[__Justification of DMD for maps with stable hyperbolic fixed points_ _]
Assume that_

  * (A1)

_$\mathbf{x}=\mathbf{0}$ is a stable hyperbolic fixed point of system _(23),
i.e., assumption (24) holds.

  * _(A2)_

_In eq. (_ 23),_ $\tilde{\mathbf{f}}\in C^{2}$ holds in a neighborhood of the
origin._

  * _(A3)_

_For some integer $d\in\left[1,n\right]$, a $d$-dimensional observable
function $\boldsymbol{\phi}\in C^{2}$ and the $d$-dimensional slow spectral
subspace $E$ of the hyperbolic fixed point $\mathbf{x}=\mathbf{0}$ of system
(_23 _) satisfy the non-degeneracy condition_

$\mathrm{rank\,}\left[D\boldsymbol{\phi}\left(\mathbf{0}\right)|_{E}\right]=d.$
(26)

  * _(A4)_

_The data matrices $\boldsymbol{\Phi}$ and $\hat{\boldsymbol{\Phi}}$ collected
from iterations of system (_23 _) are non-degenerate and are dominated by data
near $E$, i.e.,_

$\mathrm{rank_{\mathrm{row}}\,}\boldsymbol{\Phi}=d,\qquad\left|\mathbf{P}_{F}^{\mathrm{}}\mathbf{X}\right|,\left|\mathbf{P}_{F}^{\mathrm{}}\hat{\mathbf{X}}\right|\leq\left|\mathbf{P}_{E}^{\mathrm{}}\mathbf{X}\right|^{1+\beta},$
(27)

_for some $\beta\in(0,1)$._

_Then the DMD computed from $\boldsymbol{\Phi}$ and $\hat{\boldsymbol{\Phi}}$
yields a matrix $\boldsymbol{\mathcal{D}}$ that is locally topologically
conjugate with order
$\mathcal{O}\left(\left|\mathbf{P}_{E}^{\mathrm{}}\mathbf{X}\right|^{\beta}\right)$
error to the linearized dynamics on a $d$-dimensional attracting spectral
submanifold _$\mathcal{W}(E)\in C^{1}$_tangent to $E$ at
$\mathbf{x}=\mathbf{0}$. Specifically, we have_

$\boldsymbol{\mathcal{D}}=D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\mathbf{\Lambda}_{E}\left(D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\right)^{-1}+\mathcal{O}\left(\left|\mathbf{P}_{E}^{\mathrm{}}\mathbf{X}\right|^{\beta}\right).$
(28)

_The spectral submanifold $\mathcal{W}(E)$ and its reduced dynamics are of
class $C^{1}$ at the origin, and at least Hölder continuous in a neighborhood
of the origin._

###### Proof.

The proof is identical to the proof of Theorem 1 but uses the discrete version
of the linearization result by Hartman [23] for stable hyperbolic fixed points
of maps. ∎

Theorem 1 can be immediately applied to justify DMD as a linearization tool
for period-one maps (or Poincaré maps) of time-periodic, non-autonomous
dynamical systems near their periodic orbits. This requires the data matrices
_$\boldsymbol{\Phi}$ _and _ $\hat{\boldsymbol{\Phi}}$_ to contain trajectories
of such a Poincaré map. Remark 1 on the treatment of slow spectral subspaces
$E$ containing possible instabilities also applies here under the modified
assumption

$\left|\text{$\lambda_{n}$}\right|\leq\ldots\leq\left|\text{$\lambda_{d+1}$
}\right|<1,\qquad\left|\text{$\lambda_{d+1}$
}\right|<\left|\text{$\lambda_{d}$
}\right|\leq\ldots\leq\left|\text{$\lambda_{1}$
}\right|,\quad\left|\text{$\lambda_{j}$ }\right|\neq 1,\quad j=1,\ldots,n,$
(29)

which only requires the fixed point to be hyperbolic and $\mathbf{A}$ to have
a $d$-dimensional normally attracting subspace.

###  3.3 Justification of DMD for infinite-dimensional dynamical systems

Most data sets of interest arguably arise from infinite-dimensional dynamical
systems of fluids and solids. Examples include experimental or numerical data
describing fluid motion, continuum vibrations, climate dynamics or salinity
distribution in the ocean. In the absence of external forcing, these problems
are governed by systems of autonomous nonlinear partial differential equations
that can often be viewed as evolutionary differential equations in a form
similar to eq. (1), but defined on an appropriate infinite-dimensional Banach
space. Accordingly, time-sampled solutions of these equations can be viewed as
iterated mappings of the form (23) but defined on Banach spaces.

Our approach to justifying DMD generally carries over to this infinite-
dimensional setting, as long as the observable vector
$\boldsymbol{\phi}(\mathbf{x})$ remains finite-dimensional, and both the
Banach space and the discrete or continuous dynamical system defined on it
satisfy appropriate regularity conditions. These regularity conditions tend to
be technical, but when they are satisfied, they do guarantee the extension of
Theorems 1-2 to Banach spaces. This offers a justification to use DMD to
obtain an approximate finite-dimensional linear model for the dynamics of the
underlying continuum system on a finite-dimensional attracting slow manifold
(or inertial manifold) in the neighborhood of a non-degenerate stationary
solution.

To avoid major technicalities, we only state here a generalized version of
Theorem 1 to justify the use of DMD for observables defined on Banach spaces
for a discrete evolutionary process with a stable hyperbolic stationary state.
We consider mappings of the form

$\mathbf{x}_{n+1}=\mathbf{f}(\mathbf{x}_{n})=\mathbf{A}\mathbf{x}_{n}+\tilde{\mathbf{f}}(\mathbf{x}_{n}),\qquad\mathbf{x}_{j}\in\mathcal{B},\qquad\tilde{\mathbf{f}}\colon
U\subset\mathcal{B}\to\mathcal{B},\qquad\tilde{\mathbf{f}}(\mathbf{0})=\mathbf{0}\in
U,$ (30)

where $\mathcal{B}$ is a Banach space, $U$ is an open set in $\mathcal{B}$,
and $\mathbf{A}\colon\mathcal{B}\to\mathcal{B}$ is an invertible linear
operator that is bounded in the norm defined on $\mathcal{B}$. The function
$\mathbf{f}$ can be here the time-sampled version of an infinite-dimensional
flow map of an autonomous evolutionary PDE or the Poincaré map of a time-
periodic evolutionary PDE. We assume that for some
$\alpha\in\left(0,1\right)$, $\tilde{\mathbf{f}}\in C^{1,\alpha}(U)$ holds,
i.e., $\tilde{\mathbf{f}}$ is (Fréchet-) differentiable in $U$ and its
derivative, $D\tilde{\mathbf{f}}$, is Hölder-continuous in $\mathbf{x}\in U$
with Hölder exponent $\alpha$.

The _spectral radius_ of $A$ is defined as

$\rho(\mathbf{A})=\lim_{k\to\infty}\left|\mathbf{A}^{k}\right|^{\frac{1}{k}}.$

We recall that in the special case $\mathcal{B}=\mathbb{R}^{n}$ treated in
Section 3.2, we have $\rho(\mathbf{A})=\max_{1\leq j\leq
n}\left|\text{$\lambda_{j}$}\right|$. For some $\alpha\in\left(0,1\right)$,
the linear operator $\mathbf{A}$ is called _$\alpha$ -contracting_ if

$\rho(\mathbf{A})^{1+\alpha}\rho\left(\mathbf{A}^{-1}\right)<1,$ (31)

which can only hold if $\rho(\mathbf{A})<1$ (see Newhouse [41]). Therefore, in
the simple case of $\mathcal{B}=\mathbb{R}^{n}$, $A$ is _$\alpha$
-contracting_ if it is a contraction (i.e., all its eigenvalues are less than
one in norm) and

$\left|\lambda_{1}\right|^{1+\alpha}<\left|\lambda_{n}\right|,$

showing that the spectrum of $\mathbf{A}$ is confined to an annulus of outer
radius $\left|\lambda_{1}\right|<1$ and inner radius
$\left|\lambda_{1}\right|^{1+\alpha}$. We can now state our main result on the
justification of DMD for infinite-dimensional discrete dynamical systems.

######  Theorem 3.

_[__Justification of DMD for infinite-dimensional maps with stable hyperbolic
fixed points_ _] Assume that_

  * (A1)

For some $\alpha\in\left(0,1\right)$, the linear operator $\mathbf{A}$ is
_$\alpha$ -contracting _(and hence the $x=0$ fixed point of system (30) is
linearly stable).

  * _(A2)_

_In eq. (_ 30),__ $\tilde{\mathbf{f}}\in C^{1,\alpha}(U)$_holds in a $U$
neighborhood of the origin._

  * _(A3)_

_For some integer $d\in\mathbb{N}^{+}$, there is a splitting
$\mathcal{B}=E\oplus F$ of $\mathcal{B}$ into two $A$-invariant subspaces
$E,F\subset\mathcal{B}$ such that $E$ is $d$-dimensional and slow, i.e.,_

$\rho\left(\mathbf{A}|_{E}\right)<\frac{1}{\rho\left(\mathbf{A}^{-1}|_{F}\right)}$
(32)

_Furthermore, a $d$-dimensional observable function $\boldsymbol{\phi}\in
C^{2}$ satisfies the non-degeneracy condition_

$\mathrm{rank\,}\left[D\boldsymbol{\phi}\left(\mathbf{0}\right)|_{E}\right]=d.$
(33)

  * _(A4)_

_The data matrices $\boldsymbol{\Phi}$ and $\hat{\boldsymbol{\Phi}}$ collected
from iterations of system (_30 _) are non-degenerate and are dominated by data
near $E$, i.e.,_

$\mathrm{rank_{\mathrm{row}}\,}\boldsymbol{\Phi}=d,\qquad\left|\mathbf{P}_{F}^{\mathrm{}}\mathbf{X}\right|,\left|\mathbf{P}_{F}^{\mathrm{}}\hat{\mathbf{X}}\right|\leq\left|\mathbf{P}_{E}^{\mathrm{}}\mathbf{X}\right|^{1+\beta},$
(34)

_for some $\beta\in(0,1)$._

_Then the DMD computed from $\boldsymbol{\Phi}$ and $\hat{\boldsymbol{\Phi}}$
yields a matrix $\boldsymbol{\mathcal{D}}$ that is locally topologically
conjugate with order
$\mathcal{O}\left(\left|\mathbf{P}_{E}^{\mathrm{}}\mathbf{X}\right|^{\beta}\right)$
error to the linearized dynamics on a $d$-dimensional attracting spectral
submanifold $\mathcal{W}(E)$ tangent to $E$ at $\mathbf{x}=\mathbf{0}$.
Specifically, we have_

$\boldsymbol{\mathcal{D}}=D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\mathbf{\Lambda}_{E}\left(D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\right)^{-1}+\mathcal{O}\left(\left|\mathbf{P}_{E}^{\mathrm{}}\mathbf{X}\right|^{\beta}\right).$
(35)

_The spectral submanifold $\mathcal{W}(E)$ and its reduced dynamics are of
class $C^{1}$ in a neighborhood of the origin._

###### Proof.

The proof follows the steps in the proof of Theorem 2 but uses an infinite-
dimensional linearization result, Theorem 3.1 of Newhouse [41], for stable
hyperbolic fixed points of maps on Banach spaces. Specifically, if
$\mathbf{A}$ is $\alpha$-contracting, then Newhouse [41] shows the existence
of a near-identity linearizing transformation
$\mathbf{x}=\mathbf{y}+\mathbf{h}(\mathbf{y})$ for the discrete dynamical
system (30) such that $\mathbf{h}\in C^{1,\alpha}(B)$ holds on a small enough
ball $B\subset U$ centered at $\mathbf{x}=\mathbf{0}$. Using this
linearization theorem instead of its finite-dimensional version from Hartman
[23], we can follow the same steps as in the proof of Theorem 2 to conclude
the statement of the theorem. ∎

In Appendix B, Remarks 6-7 summarize technical remarks on possible further
extensions of Theorem 3.

##  4 Data-Driven Linearization (DDL)

###  4.1 Theoretical foundation for DDL

Based on the results of the previous section, we now refine the first-order
approximation to the linearized dynamics yielded by DMD near a hyperbolic
fixed point. Specifically, we construct the specific nonlinear coordinate
change that linearizes the restricted dynamics on the attracting spectral
submanifold $\mathcal{W}(E)$ illustrated in Fig. 1. This classic notion of
linearization on $\mathcal{W}(E)$ yields a $d$-dimensional linear reduced
model, which can be of significantly lower dimension than the original
$n$-dimensional nonlinear system. This is to be contrasted with the broadly
pursued Koopman embedding approach (see, e.g., Rowley et al. [49], Budišić et
al. [8], Mezić [40]), which seeks to immerse nonlinear systems into linear
systems of dimensions substantially higher (or even infinite) relative to $n$.

The following result gives the theoretical basis for our subsequent data-
driven linearization (DDL) algorithm. We will use the notation $C^{a}$ for the
class of real analytic functions.

######  Theorem 4.

_[__DDL principle for ODEs with a stable hyperbolic fixed points_ _] Assume
that the $\mathbf{x}=\mathbf{0}$ stable hyperbolic fixed point of system (11)
and the spectrum of $\mathbf{A}$ has a spectral gap as in eq. (12). Assume
further that for some $r\in\mathbb{N}^{+}\cup\left\\{\infty,a\right\\}$, the
following conditions are satisfied:_

  * (B1)

_$\mathbf{f}_{2}\in C^{r}$ and the nonresonance conditions _

$\lambda_{k}\neq\sum_{j=1}^{n}m_{j}\lambda_{j},\qquad j,k=1,\ldots,n,\qquad
2\leq\sum_{j=1}^{n}m_{j}\leq Q\leq r,\quad
Q:=\mathrm{Int}\left[\frac{\max_{i}\left|\mathrm{Re}\,\lambda_{i}\right|}{\min_{i}\left|\mathrm{Re}\,\lambda_{i}\right|}\right]+1,$
(36)

_hold for the eigenvalues of $\mathbf{A}$._

  * (B2)

_For some integer $d\in\left[1,n\right]$, a $d$-dimensional observable
function $\boldsymbol{\phi}\in C^{r}$ and the $d$-dimensional slow spectral
subspace $E$ of the stable fixed point $\mathbf{x}=\mathbf{0}$ of system (_11
_) satisfy the non-degeneracy condition._

$\mathrm{rank\,}\left[D\boldsymbol{\phi}\left(\mathbf{0}\right)|_{E}\right]=d.$
(37)

_Then the following hold:_

_(i) On the unique $d$-dimensional attracting spectral submanifold
$\mathcal{W}(E)\in C^{r}$ tangent to $E$ at $x=0$, the reduced observable
vector $\boldsymbol{\varphi}=\boldsymbol{\phi}|_{\mathcal{W}(E)}$ can be used
to describe the reduced dynamics as_

$\dot{\boldsymbol{\varphi}}=\mathbf{B}\boldsymbol{\varphi}+\mathbf{q}\left(\boldsymbol{\varphi}\right),\qquad\mathbf{B}=D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\mathbf{\Lambda}_{E}\left(D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\right)^{-1},\quad\mathbf{q}\left(\boldsymbol{\varphi}\right)=\mathcal{O}\left(\left|\boldsymbol{\varphi}\right|^{2}\right).$
(38)

_(ii) There exists a unique, $C^{r}$ change of coordinates_

$\boldsymbol{\varphi}=\boldsymbol{\kappa}(\boldsymbol{\gamma})=\boldsymbol{\gamma}+\boldsymbol{\ell}(\boldsymbol{\gamma}),$
(39)

_that transforms the reduced dynamics on $\mathcal{W}(E)$ to its
linearization_

$\dot{\boldsymbol{\gamma}}=\mathbf{B}\boldsymbol{\gamma}$ (40)

_inside the domain of attraction of $\mathbf{x}=\mathbf{0}$ within the
spectral submanifold $\mathcal{W}(E).$ _

_(iii) The transformation (_ 39 _) satisfies the $d$-dimensional system of
nonlinear PDEs_

$D_{\boldsymbol{\gamma}}\boldsymbol{\ell}(\boldsymbol{\gamma})\mathbf{B}\boldsymbol{\gamma}=\mathbf{B}\boldsymbol{\ell}(\boldsymbol{\gamma})+\mathbf{q}\left(\boldsymbol{\gamma}+\boldsymbol{\ell}(\boldsymbol{\gamma})\right).$
(41)

_If $r\in\mathbb{N}^{+}\cup\left\\{\infty\right\\}$, solutions of this PDE can
locally be approximated as_

$\boldsymbol{\ell}(\boldsymbol{\gamma})=\sum_{\left|\mathbf{k}\right|=2}^{r}\mathbf{l}_{\mathbf{k}}\gamma^{\mathbf{k}}+o\left(\left|\boldsymbol{\gamma}\right|^{r}\right),\quad\mathbf{k}\in\mathbb{N}^{d},\quad\mathbf{l}_{\mathbf{k}}\in\mathbb{R}^{d},\quad\boldsymbol{\gamma}^{\mathbf{k}}:=\gamma_{1}^{k_{1}}\cdots\gamma_{d}^{k_{d}}.$
(42)

_If $r=a$, then the local approximation (_42) _can be refined to a convergent
Taylor series_

$\boldsymbol{\ell}(\boldsymbol{\gamma})=\sum_{\left|\mathbf{k}\right|=2}^{\infty}\mathbf{l}_{\mathbf{k}}\gamma^{\mathbf{k}},\quad\mathbf{k}\in\mathbb{N}^{d},\quad\mathbf{l}_{\mathbf{k}}\in\mathbb{R}^{d},\quad\boldsymbol{\gamma}^{\mathbf{k}}:=\gamma_{1}^{k_{1}}\cdots\gamma_{d}^{k_{d}}$
(43)

_in a neighborhood of the origin. In either case, the coefficients
$\mathbf{l}_{\mathbf{k}}$ can be determined by substituting the expansion for
$\boldsymbol{\ell}(\boldsymbol{\gamma})$ into the PDE (41), equating
coefficients of equal monomials $\boldsymbol{\gamma}^{\mathbf{k}}$ and solving
the corresponding recursive sequence of $d$-dimensional linear algebraic
equations for increasing $\left|\mathbf{k}\right|$._

###### Proof.

The proof builds on the existence of the $d$-dimensional spectral submanifold
$\mathcal{W}\left(E\right)$ guaranteed by Theorem 1. For a $C^{r}$ dynamical
system with $r\in\mathbb{N}^{+}\cup\left\\{\infty,a\right\\}$,
$\mathcal{W}\left(E\right)$ is also $C^{r}$ smooth based on the linearization
theorems of Poincaré [46] and Sternberg [52], as long as the nonresonance
condition (36) holds. Condition (37) then ensures that
$\mathcal{W}\left(E\right)$ can be parametrized locally by the restricted
observable vector $\boldsymbol{\varphi}$ and hence its reduced dynamics can be
written as a nonlinear ODE for $\boldsymbol{\varphi}$. This ODE can again be
linearized by a near-identity coordinate change (39) using the appropriate
linearization theorem of the two cited above. The result is the restricted
linear system (40) to which the dynamics is $C^{r}$ conjugate within the whole
domain of attraction of the $\boldsymbol{\varphi}=\mathbf{0}$ fixed point
inside $\mathcal{W}\left(E\right)$. The invariance PDE (41) can be obtained by
substituting the linearizing transformation (39) into the reduced dynamics on
$\mathcal{W}\left(E\right)$. This PDE can then be solved via a Taylor
expansion up to order $r$. We give more a more detailed proof in Appendix D. ∎

Note that $1:1$ resonances are not excluded by the condition (36), and hence
repeated eigenvalues arising from symmetries in physical systems are still
amenable to DDL. Also of note is that the non-resonance conditions (36) do not
exclude frequency-type resonances among imaginary parts of oscillatory
eigenvalues. Rather, they exclude simultaneous resonances of the same type
between the real and the imaginary parts of the eigenvalues. Such resonances
will be absent in data generated by generic oscillatory systems.

Assuming hyperbolicity is essential for Theorem ‘4 to hold, since in this case
the linearization is the same as transforming the dynamics to the Poincaré-
normal form. For a non-hyperbolic fixed point, this normal form transformation
results in nonlinear dynamics on the center manifold. This would, however,
only arise in highly non-generic systems, precisely tuned to be at
criticality. Since this is unlikely to happen in experimentally observed or
numerically simulated systems, the hyperbolicity assumption is not
restrictive.

Finally, under the conditions of Theorem 3, the DDL results of Theorem 4 also
apply to data from infinite-dimensional dynamical systems, such as the fluid
sloshing experiments we will analyze using DDL in Section 5.4. In practice,
the most restrictive condition of Theorem 3 is (A1), which requires the
solution operator to have a spectrum uniformly bounded away from zero. Such
uniform boundedness is formally violated in important classes of infinite-
dimensional evolution equations, presenting a technical challenge for the
direct applications of SSM results to certain delay-differential equations
(see Szaksz [54]) and partial differential equations (see, e.g., Kogelbauer
and Haller [30] and Buza [9]). However, this challenge only concerns rigorous
conclusions on the existence and smoothness of a finite-dimensional,
attracting SSM. If the existence of such an SSM is convincingly established
from an alternative mathematical theory (as is Buza [9]) or inferred from data
(as in Szaksz [54]), then the DDL algorithm based on Theorem 4 can be used to
obtain a data-driven linearization of the dynamics on that SSM.

###  4.2 DDL vs. EDMD

Here we examine whether there is a possible relationship between DDL and the
extended DMD (or EDMD) algorithm of Williams et al. [58]. For simplicity, we
assume analyticity for the dynamical system ($r=a$) and hence we can write the
inverse of the linearizing transformation (42) behind the DDL algorithm as a
convergent Taylor expansion of the form

$\boldsymbol{\gamma}=\boldsymbol{\kappa}^{-1}\left(\boldsymbol{\varphi}\right)=\boldsymbol{\varphi}+\sum_{\left|\mathbf{k}\right|=2}^{\infty}\mathbf{q}_{\mathbf{k}}\boldsymbol{\varphi}^{\mathbf{k}}.$
(44)

We then differentiate this equation in time to obtain from the linearized
equation (40) a $d$-dimensional system of equations

$\dot{\boldsymbol{\varphi}}+\sum_{\left|\mathbf{k}\right|=2}^{\infty}\mathbf{q}_{\mathbf{k}}\frac{d}{dt}\boldsymbol{\varphi}^{\mathbf{k}}=\mathbf{B}\boldsymbol{\varphi}+\sum_{\left|\mathbf{k}\right|=2}^{\infty}\mathbf{B}\mathbf{q}_{\mathbf{k}}\boldsymbol{\varphi}^{\mathbf{k}}$

that the restricted observable $\boldsymbol{\varphi}$ and its monomials
$\boldsymbol{\varphi}^{\mathbf{k}}$ must satisfy. This last equation can be
rewritten as a $d$-dimensional autonomous system of linear system of ODEs,

$\left[\begin{array}[]{cc}\mathbf{I}_{d\times
d}&\mathbf{Q}_{2}\end{array}\right]\frac{d}{dt}\left[\begin{array}[]{c}\boldsymbol{\varphi}\\\
\mathbf{K}_{\geq
2}\left(\boldsymbol{\varphi}\right)\end{array}\right]=\left[\begin{array}[]{cc}\mathbf{B}&\mathbf{B}\mathbf{Q}_{2}\end{array}\right]\left[\begin{array}[]{c}\boldsymbol{\varphi}\\\
\boldsymbol{K}_{\geq 2}\left(\boldsymbol{\varphi}\right)\end{array}\right],$
(45)

for the reduced observable $\boldsymbol{\varphi}$ and the infinite-dimensional
vector $\mathbf{K}_{\geq 2}\left(\varphi\right)$ of all nonlinear monomials of
$\boldsymbol{\varphi}$. Here $\mathbf{I}_{d\times d}$ denotes the
$d$-dimensional identity matrix and $\mathbf{Q}_{2}$ contains all coefficients
$\mathbf{q}_{\mathbf{k}}$ as column vectors starting from order
$\left|\mathbf{k}\right|=2$.

If we truncate the infinite-dimensional vector of monomials $\mathbf{K}_{\geq
2}\left(\boldsymbol{\varphi}\right)$ to the vector
$\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)$ of nonlinear monomials
up to order $k$, then eq. (45) becomes

$\left[\begin{array}[]{cc}\mathbf{I}_{d\times
d}&\mathbf{Q}_{2}^{k}\end{array}\right]\frac{d}{dt}\left[\begin{array}[]{c}\boldsymbol{\varphi}\\\
\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)\end{array}\right]=\left[\begin{array}[]{cc}\mathbf{B}&\mathbf{B}\mathbf{Q}_{2}^{k}\end{array}\right]\left[\begin{array}[]{c}\boldsymbol{\varphi}\\\
\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)\end{array}\right].$ (46)

This is a $d$-dimensional implicit system of linear ODEs for the dependent
variable vector
$\left(\boldsymbol{\varphi},\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)\right)$
whose dimension is always larger than $d$. Consequently, the operator
$\left[\begin{array}[]{cc}\mathbf{I}_{d\times
d}&\mathbf{Q}_{2}^{k}\end{array}\right]$ is never invertible and hence,
contrary to the assumption of EDMD, there is no well-defined linear system of
ODEs that governs the evolution of an observable vector and the monomials of
its components.

The above conclusion remains unchanged even if one attempts to optimize with
respect to the choice of the coefficients $\mathbf{q}_{\mathbf{k}}$ in the
matrix $\mathbf{Q}_{2}^{k}$.

###  4.3 Implementation and applications of DDL

####  4.3.1 Basic implementation of DDL for model reduction and linearization

Theorem 4 allows us to define a numerical procedure to construct a linearizing
transformation on the $d$-dimensional attracting slow manifold _
$\mathcal{W}(E)$ _systematically from data. From eq. (46), the matrices
$\mathbf{B}$ and $\mathbf{Q}$ are to be determined, given a set of observed
trajectories. In line with the notation used in Section 4.2, let the data
matrix $\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)$ contain monomials
(from order $2$ to order $k)$ of the observable vector $\boldsymbol{\varphi}$
and let $\hat{\mathbf{K}_{2}^{k}}\left(\boldsymbol{\varphi}\right)$ contain
denote the evaluation of $\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)$
time $\Delta t$ later. Passing to the discrete version of the invariance
equation (46), we obtain

$\left[\begin{array}[]{cc}\mathbf{I}&\mathbf{Q}\end{array}\right]\left[\begin{array}[]{c}\hat{\boldsymbol{\varphi}}\\\
\hat{\mathbf{K}_{2}^{k}}\left(\boldsymbol{\varphi}\right)\end{array}\right]=\left[\begin{array}[]{cc}\boldsymbol{\mathcal{B}}&\boldsymbol{\mathcal{B}}\mathbf{Q}\end{array}\right]\left[\begin{array}[]{c}\boldsymbol{\varphi}\\\
\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)\end{array}\right]$

for some matrices $\mathbf{Q}\in\mathbb{R}^{d\times N(d,k)-d}$ and
$\boldsymbol{\mathcal{B}}=e^{\mathbf{B}\Delta t}\in\mathbb{R}^{d\times d}$.
Moreover, the inverse transformation of the linearization on the SSM
$\mathcal{W}\left(E\right)$ is well-defined, and hence with an appropriate
matrix $\mathbf{Q}^{inv}\in\mathbb{R}^{d\times N(d,k)-d}$, we can write

$\left[\begin{array}[]{cc}\mathbf{I}&\mathbf{Q}^{inv}\end{array}\right]\left[\begin{array}[]{c}\boldsymbol{\varphi}+\mathbf{Q}\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)\\\
\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}+\mathbf{Q}\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)\right)\end{array}\right]=\boldsymbol{\varphi}.$

This allows us to define the cost functions

$\displaystyle\mathcal{L}^{(1)}(\mathbf{Q},\boldsymbol{\mathcal{B}})$
$\displaystyle=\left|\left[\begin{array}[]{cc}\mathbf{I}&\mathbf{Q}\end{array}\right]\left[\begin{array}[]{c}\hat{\boldsymbol{\varphi}}\\\
\hat{\mathbf{K}_{2}^{k}}\left(\boldsymbol{\varphi}\right)\end{array}\right]-\left[\begin{array}[]{cc}\boldsymbol{\mathcal{B}}&\boldsymbol{\mathcal{B}}\mathbf{Q}\end{array}\right]\left[\begin{array}[]{c}\boldsymbol{\varphi}\\\
\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)\end{array}\right]\right|^{2},$
(53) $\displaystyle\mathcal{L}^{(2)}\left(\mathbf{Q},\mathbf{Q}^{inv}\right)$
$\displaystyle=\left|\mathbf{Q}\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)+\mathbf{Q}^{inv}\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}+\mathbf{Q}\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)\right)\right|^{2},$

where $\mathcal{L}^{(1)}$ measures the invariance error along the observed
trajectories and $\mathcal{L}^{(2)}$ measures the error due to the computation
of the inverse. We aim to jointly minimize $\mathcal{L}^{(1)}$ and
$\mathcal{L}^{(2)}$. To this end, we define the combined cost function

$\mathcal{L}_{\nu}\left(\mathbf{Q},\mathbf{Q}^{inv},\boldsymbol{\mathcal{B}}\right)=\mathcal{L}^{(1)}(\mathbf{Q},\boldsymbol{\mathcal{B}})+\nu\mathcal{L}^{(2)}\left(\mathbf{Q},\mathbf{Q}^{inv}\right),$
(54)

for some $\nu\geq 0$. In our examples, we choose $\nu=1$, which puts the same
weight on both terms in the cost function (54). Minimizers of
$\mathcal{L}_{\nu}$ provide optimal solutions to the DDL principle and can be
written as

$(\mathbf{Q}^{\star},\mathbf{Q}^{inv,\star}\boldsymbol{\mathcal{B}}^{\star})=\underset{\mathbf{Q},\mathbf{Q}^{inv},\boldsymbol{\mathcal{B}}}{\mathrm{argmin}\,\,}\mathcal{L_{\nu}}(\mathbf{Q},\mathbf{Q}^{inv},\boldsymbol{\mathcal{B}}),$
(55)

or, equivalently, as solutions of the system of equations

$\displaystyle\frac{\partial\mathcal{L}_{\nu}}{\partial Q_{ij}}$
$\displaystyle=0\quad i=1,...,d,j=1,...,N(d,k)-d,$ (56)
$\displaystyle\frac{\partial\mathcal{L}_{\nu}}{\partial Q_{ij}^{inv}}$
$\displaystyle=0,\quad i=1,...,d,j=1,...,N(d,k)-d,$ (57)
$\displaystyle\frac{\partial\mathcal{L}_{\nu}}{\partial\mathcal{B}_{ij}}$
$\displaystyle=0\quad i,j=1,...,d.$

The optimal solution (55) does not necessarily coincide with the Taylor-
coefficients of the linearizing transformation (43). Instead of giving the
best local approximation,
$(\mathbf{Q}^{\star},\mathbf{Q}^{inv,\star}\boldsymbol{\mathcal{B}}^{\star})$
approximates the linearizing transformation and the linear dynamics in a
least-squares sense over the domain of the training data. This means that DDL
is not hindered by the convergence properties of the analytic linearization.
Note that for $d=1$, one can estimate the radius of convergence of (43), for
example, by constructing the Domb-Skyes plot (see Domb and Sykes [15]), or by
finding the radius of the circle in the complex plane onto which the roots of
the truncated expansion accumulate under increasing orders of truncation (see
Jentzsch [25], Ponsioen et al. [47]). For $d>1$, such analysis is more
difficult, since multivariate Taylor-series have more complicated domains of
convergence. In our numerical examples, we estimate the domain of convergence
of such analytic linearizations as the domain on which
$\boldsymbol{\kappa}\circ\boldsymbol{\kappa}^{-1}=\mathbf{I}$ holds to a good
approximation. As we will see, this domain of convergence may be substantially
smaller the domain of validity of transformations determined in a fully data-
driven way.

Since the cost function (53) is not convex, the optimization problem (55) has
to be solved iteratively starting from an initial guess
$(\mathbf{Q}_{0},\mathbf{Q}_{0}^{inv},\boldsymbol{\mathcal{B}}_{0})$. For the
examples presented in the paper, we use the Levenberg–Marquardt algorithm (see
Bates and Watts [4]), but other nonlinear optimization methods, such as
gradient descent or Adam (see Kingma and Ba [29]) could also be used. For our
implementation, which is available from the repository Kaszás et al. [28], we
used the Scipy and Pytorch libraries of Virtanen et al. [57], Paszke et al.
[45]. In summary, we will use the following Algorithm 1 in our examples for
model reduction via DDL.

Data:

$d$: Dimension of the slow SSM $\mathcal{W(E)}$.

$k$: Maximal polynomial order.

$\boldsymbol{\varphi},\hat{\boldsymbol{\varphi}}\in\mathbb{R}^{d\times
n_{samples}}$: Data matrices of the reduced coordinates and their forward-
shifted images. The input data may need to be truncated, in order to ensure
that it lies sufficiently close to the SSM. This needs to be checked, e.g., by
time-frequency analysis (see Cenedese et al. [12]).

$\nu\geq 0$: Weight parameter in the cost function

$tol$: Tolerance value for the cost function.

Result:

$\mathbf{Q}\in\mathbb{R}^{d\times N(d,k)-d}$: Coefficients of the
transformation
$\boldsymbol{\varphi}=\boldsymbol{\kappa}(\boldsymbol{\gamma})$.

$\mathbf{Q}^{inv}\in\mathbb{R}^{d\times N(d,k)-d}$: Coefficients of the
inverse transformation
$\boldsymbol{\gamma}=\boldsymbol{\kappa}^{-1}(\boldsymbol{\varphi})$.

$\boldsymbol{\mathcal{B}}$: The linear dynamics.

1 Compute the monomials $\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)$,
$\hat{\mathbf{K}_{2}^{k}}\left(\boldsymbol{\varphi}\right)\in\mathbb{R}^{N(d,k)-d\times
n_{samples}}$

2Choose an initial guess, either randomly or inferred from DMD

3
$(\mathbf{Q},\mathbf{Q}^{inv},\boldsymbol{\mathcal{B}})\leftarrow(\mathbf{Q}_{0},\mathbf{Q}_{0}^{inv},\boldsymbol{\mathcal{B}}_{0})$

4while
_$\mathcal{L}_{\nu}(\mathbf{Q},\mathbf{Q}^{inv},\boldsymbol{\mathcal{B}})
>tol$_ do

5 Compute the monomials
$\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}+\mathbf{Q}\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)\right)$

6 Evaluate
$\mathcal{L}_{\nu}\left(\mathbf{Q},\mathbf{Q}^{inv},\boldsymbol{\mathcal{B}}\right)=\mathcal{L}^{(1)}(\mathbf{Q},\boldsymbol{\mathcal{B}})+\nu\mathcal{L}^{(2)}\left(\mathbf{Q},\mathbf{Q}^{inv}\right)$

7 Perform Optimization step (e.g., Levenberg–Marquardt or gradient descent)

8
$(\mathbf{Q},\mathbf{Q}^{inv},\boldsymbol{\mathcal{B}})\leftarrow(\mathbf{Q}_{new},\mathbf{Q}_{new}^{inv},\boldsymbol{\mathcal{B}}_{new})$

9 end while

Algorithm 1 Model reduction with DDL

######  Remark 3.

The expressions (53)-(54) define one of the possible choices for the cost
function. With $\nu=0$, (54) simply corresponds to a one-step-ahead prediction
with the linearized dynamics. Alternatively, a multi-step prediction can also
be enforced. For a training trajectory $\boldsymbol{\varphi}(t),$ the
invariance

$\left[\begin{array}[]{cc}\mathbf{I}&\mathbf{Q}\end{array}\right]\left[\begin{array}[]{c}\boldsymbol{\varphi}\\\
\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}\right)\end{array}\right]=\left[\begin{array}[]{cc}\boldsymbol{\mathcal{B}}^{1:m}&\boldsymbol{\mathcal{B}}^{1:m}\mathbf{Q}\end{array}\right]\left[\begin{array}[]{c}\boldsymbol{\varphi}(\mathbf{0})\\\
\mathbf{K}_{2}^{k}\left(\boldsymbol{\varphi}(\mathbf{0})\right)\end{array}\right]$

could be required, where $\boldsymbol{\mathcal{B}}^{1:m}$ is a tensor composed
of powers of the linear map $\boldsymbol{\mathcal{B}}$. Optimizing over the
entire trajectory is, however, more costly than simply minimizing (53), and we
found no noticeable improvement in accuracy in our numerical examples.

####  4.3.2 Relationship with DMD implementations

Note that setting $\mathbf{Q}=\mathbf{Q}^{inv}=\mathbf{0}$ in the optimization
problem (55) turns the problem into DMD. In this case, the usual DMD algorithm
surveyed in the Introduction returns

$\mathcal{\boldsymbol{\mathcal{B}}}_{0}=\underset{\boldsymbol{\mathcal{B}}}{\mathrm{argmin}\,\,}\underset{}{\mathcal{L}_{0}(\mathbf{0},\mathbf{0},\boldsymbol{\mathcal{B}})},$

which is a good initial guess for the non-convex optimization problem (55).
More importantly, since Theorem 4 guarantees the existence of a near-identity
linearizing transformation, we expect that the true minimizer is close to the
DMD-solution. Therefore, we may explicitly expand the cost function (53)
around the DMD solution as

$\displaystyle\mathcal{L_{\nu}}(\mathbf{Q},\mathbf{Q}^{inv},\mathbf{\mathcal{B}})$
$\displaystyle=\mathcal{L}_{\nu}(\mathbf{0},\mathbf{0},\boldsymbol{\mathcal{B}}_{0})+D\mathcal{L}_{(\mathbf{0},\mathbf{0},\mathcal{\boldsymbol{\mathcal{B}}}_{0})}\cdot\left(\begin{array}[]{c}\mathbf{Q}\\\
\mathbf{Q}^{inv}\\\
\mathcal{\boldsymbol{\mathcal{B}}}-\boldsymbol{\mathcal{B}}_{0}\end{array}\right)$
$\displaystyle+\frac{1}{2}\left[D^{2}\mathcal{L}_{(\mathbf{0},\mathbf{0},\mathcal{\boldsymbol{\mathcal{B}}}_{0})}\cdot\left(\begin{array}[]{c}\mathbf{Q}\\\
\mathbf{Q}^{inv}\\\
\boldsymbol{\mathcal{B}}-\boldsymbol{\mathcal{B}}_{0}\end{array}\right)\right]\cdot\left(\begin{array}[]{c}\mathbf{Q}\\\
\mathbf{Q}^{inv}\\\
\boldsymbol{\mathcal{B}}-\boldsymbol{\mathcal{B}}_{0}\end{array}\right)$
$\displaystyle+\mathcal{}\left(\left|\mathbf{Q}\right|^{3},\left|\mathbf{Q}^{inv}\right|^{3},\left|\boldsymbol{\mathcal{B}}-\boldsymbol{\mathcal{B}}_{0}\right|^{3}\right),$

where $D\mathcal{L}_{(\mathbf{0},\mathbf{0},\boldsymbol{\mathcal{B}}_{0})}$
and $D^{2}\mathcal{L}_{(\mathbf{0},\mathbf{0},\boldsymbol{\mathcal{B}}_{0})}$
are the Jacobian and the Hessian of the cost function evaluated at the DMD
solution, respectively. Since the Jacobian is nonsingular at the DMD solution,
the minimum of the quadratic approximation of the cost function satisfies the
linear equation

$-D\mathcal{L}_{(\mathbf{0},\mathbf{0},\mathcal{\boldsymbol{\mathcal{B}}}_{0})}=D^{2}\mathcal{L}_{(\mathbf{0},\mathbf{0},\boldsymbol{\mathcal{B}}_{0})}\left(\begin{array}[]{c}\mathbf{Q}\\\
\mathbf{Q}^{inv}\\\
\boldsymbol{\mathcal{B}}-\boldsymbol{\mathcal{B}}_{0}\end{array}\right).$ (58)

This serves as the first-order correction to the DMD-solution in the DDL
procedure. The equation (58) is explicitly solvable and is equivalent to
performing a single Levenberg–Marquardt step on the non-convex cost function
(53), with the DMD solution
$(\mathbf{0},\mathbf{0},\mathcal{\boldsymbol{\mathcal{B}}}_{0})$ serving as an
initial guess.

Minimization of (53) leads to a non-convex optimization problem. Besides
computing the leading-order approximation (58), a possible workaround to this
challenge is to carry out the linearization in two steps. First, one can fit a
polynomial map to the reduced dynamics by linear regression. Then, if the
reduced-dynamics is non-resonant, it can be analytically linearized. Axås et
al. [3] follow this approach to automatically find the extended normal form
style reduced dynamics on SSMs using the implementation of SSMTool by Jain et
al. [24]. Although this procedure does convert the DDL principle into a convex
problem, the drawback is that the linearization is obtained as a Taylor-
expansion, with possibly limited convergence properties.

####  4.3.3 Using DDL to construct spectral foliations

The mathematical foundation of SSM-reduced modeling is that any trajectory
converging to a slow SSM is guaranteed to synchronize up to an exponentially
decaying error with one of the trajectories on the SSM. This follows from the
general theory of invariant foliations by Fenichel [18], when applied to the
$d$-dimensional normally hyperbolic invariant manifold
$\mathcal{W}(E)$.333More precisely, Fenichel’s foliation results become
applicable after the wormhole construct in Proposition B1 of Eldering et al.
[17] is applied to extend $\mathcal{W}(E)$ smoothly into a compact normally
attracting invariant manifold without boundary. This is needed because
Fenichel’s results only apply to compact normally attracting invariant
manifolds with an empty or overflowing boundary, whereas the boundary of
$\mathcal{W}(E)$ is originally inflowing. The main result of the theory is
that off-SSM initial conditions synchronizing with the same on-SSM trajectory
turn out to form a class $C^{r-1}$ smooth, $\left(n-d\right)$-dimensional
manifold, denoted $\mathcal{F}_{\mathbf{p}}$, which intersects
$\mathcal{W}(E)$ in a unique point $\mathbf{p}\in\mathcal{W}(E)$. The manifold
$\mathcal{F}_{\mathbf{p}}$ is called the stable fiber emanating from the base
point $\mathbf{p}$. Fenichel proves that any off-SSM trajectory
$\mathbf{x}(t;\mathbf{x}_{0})$ with initial condition
$\mathbf{x}_{0}\in\text{$\mathcal{F}_{\mathbf{p}_{0}}$ }$ converges to the
specific on-SSM trajectory $\mathbf{p}(t;\mathbf{p}_{0})\in\mathcal{W}(E)$
with initial condition $\mathbf{p}_{0}\in\mathcal{W}(E)$ faster than any other
nearby trajectory might converge to $\mathbf{p}(t;\mathbf{p}_{0}).$ Recently,
Szalai [55] studied this foliation in more detail under the name “invariant
spectral foliation”, discussed its uniqueness in an appropriate smoothness
class and proposed its use in model reduction.

To predict the evolution of a specific, off-SSM initial condition
$\mathbf{x}_{0}$ up to time $t$ from an SSM-based model, we first need to
relate that initial condition to the base point $\mathbf{p}_{0}$ of the stable
fiber $\mathcal{F}_{\mathbf{p}_{0}}$. Next, we need to run the SSM-based
reduced model up to time $t$ to obtain $\mathbf{p}(t;\mathbf{p}_{0})$. Based
on the exponentially fast convergence of the full solution
$\mathbf{x}(t;\mathbf{x}_{0})$ to the SSM-reduced solution
$\mathbf{p}(t;\mathbf{p}_{0}),$ we obtain an accurate longer-term prediction
for $\mathbf{x}(t;\mathbf{x}_{0})$ using this procedure. Such a longer-term
prediction is helpful, for instance, when we wish to predict steady states,
such as fixed points and limit cycles, from the SSM-reduced dynamics.

Constructing this spectral foliation directly from data, however, is
challenging for nonlinear systems. Indeed, one would need a very large number
of initial conditions that cover uniformly a whole open neighborhood of the
fixed point in the phase space. For example, while one or two training
trajectories are generally sufficient to infer accurate SSM-reduced models
even for very high-dimensional systems (see e.g., Cenedese et al. [12, 13],
Axås et al. [3]), thousands of uniformly distributed initial conditions in a
whole open set of a fixed point are required to infer accurate spectral
foliation-based models even for low-dimensional systems (see Szalai [55]). The
latter number and distribution of initial conditions is unrealistic to acquire
in a truly data-driven setting.

To avoid constructing the full foliation, one may simply project an initial
condition $x_{0}$ orthogonally to an observed spectral submanifold
$\mathcal{W}(E)$ to obtain $\mathbf{p}_{0}$, but this may result in large
errors if $E$ and $F$ are not orthogonal. In that case, $\mathcal{W}(E)$ may
divert substantially from $E$ (see [48, 42, 43] for a discussion of the
limitations of this projection for general invariant manifolds).

A better solution is to project $\mathbf{x}_{0}$ orthogonally to the slow
spectral subspace $E$ over which $\mathcal{W}(E)$ is a graph in an (often
large) neighborhood of the fixed point. This approach assumes that $E$ and $F$
are nearly normal and $\mathcal{W}(E)$ is nearly flat. As the latter is
typically the case for delay-embedded observables (Axås and Haller [2]),
orthogonal projection onto $E$ has been the choice so far in data-driven SSM-
based reduction via the SSMLearn algorithm (Cenedese et al. [11]). This
approach has produced highly accurate reduced-order models in a number of
examples (see Cenedese et al. [12, 13], Axås et al. [3]). There are
nevertheless examples in which the linear part of the dynamical system is
significantly non-normal and hence $E$ and $F$ are not close to being
orthogonal (see Bettini et al. [6]).

Near hyperbolic fixed points, the use of DDL eliminates the need to construct
involved nonlinear spectral foliations. Indeed, let us assume that the slow
spectral subspace $E$ in Theorem 4 can be decomposed into a direct sum
$E=E_{1}\oplus E_{2}$, where $E_{1}$ denotes the slowest spectral subspace
with dim $E_{1}=d_{1}$ and $E_{2}$ denotes the second-slowest spectral
subspace with dim $E_{2}=d_{2}$, as sketched in Fig. 2.

Figure 2: (a): The linearized phase space geometry governed by the slow
spectral subspace $E=E_{1}\oplus E_{2}$ and the slow invariant foliation
within $E$. (b): Phase space geometry in the original coordinates.

Reducing the dynamics to the SSM $\mathcal{W}(E)$ is accurate for transient
times given by the decay rate of $E_{2}.$ This initial reduction can be done
simply by a normal projection onto $E$. Inside $E$, one can simply locate
spectral foliations of the DDL-linearized systems explicitly and map them back
to the original nonlinear system under the DDL transformation
$\boldsymbol{\kappa}\left(\boldsymbol{\gamma}\right)$. The unique class
$C^{a}$ foliation of a linear system within $E$ is the family of stable fibers
forming the affine space

$\mathcal{F}_{\mathbf{p}}=\mathbf{p}+E_{2},$

where $\mathbf{p}\in E_{1}$. The trajectories started inside
$\mathcal{F}_{\mathbf{p}}$ all synchronize with $\mathbf{p}\in E_{1}$. The
linear projection $\mathbf{P}_{E_{2}}$ onto $E_{1}$ along directions parallel
to $E_{2}$, when applied to an initial condition
$\mathbf{y}\in\mathcal{F}_{\mathbf{p}}$, returns the base point

$\mathbf{P}_{E_{2}}\mathbf{y}=\mathbf{p}.$ (59)

In the nonlinear system (1), the leaves of the smooth foliation within
$\mathcal{W}(E)$

$\mathcal{F}_{\boldsymbol{\kappa}\left(\mathbf{p}\right)}^{0}=\boldsymbol{\kappa}\left(\mathcal{F}_{\mathbf{p}}\right)\subset\mathcal{W}(E),$

where $\boldsymbol{\kappa}\left(\mathbf{p}\right)\in\mathcal{W}(E_{1})$ is the
image of $\mathbf{p}$ under the mapping $\boldsymbol{\kappa}$ defined in (39).
The SSM $\mathcal{W}(E)$ can then be parametrized via the foliation

$\mathcal{W}(E)=\bigcup_{\mathbf{q}\in\mathcal{W}(E_{1})}\mathcal{F}_{\mathbf{q}}^{0}.$

####  4.3.4 Using DDL to predict nonlinear forced response from unforced data

We now discuss how DDL performed near the fixed point of an autonomous
dynamical system can be used to predict nonlinear forced response under
additional weak periodic forcing in the domain of DDL. The addition of such
small forcing is frequent in structural vibration problems in which the
unforced structure (e.g., a beam or disk) is rigid enough to react with small
displacements under practically relevant excitation levels (see, e.g.,
Cenedese et al. [12, 13] for specific examples).

We append system (11) with a small, time-periodic forcing term
$\varepsilon\mathbf{F}(\mathbf{x},t)$ to obtain the system

$\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\tilde{\mathbf{f}}(\mathbf{x})+\varepsilon\mathbf{F}(\mathbf{x},t),\qquad\mathbf{x}\in\mathcal{\mathbb{R}}^{n},\qquad\mathbf{A}=D\mathbf{f}(0),\qquad\tilde{\mathbf{f}}(\mathbf{x})=\mathcal{O}\left(\left|\mathbf{x}\right|^{2}\right),\qquad
0\leq\epsilon\ll 1,$ (60)

with $\mathbf{F}(\mathbf{x},t)=\mathbf{F}(\mathbf{x},t+T)$ for some period
$T>0$. If the conditions of Theorem 4 hold for the system (60) for
$\varepsilon=0$_,_ then, for $\epsilon>0$ small enough, exists a unique
$d$-dimensional, $T$-periodic, attracting spectral submanifold
$\mathcal{W}_{\varepsilon}(E,t)\in C^{r}$ of a locally unique attracting
$T$-periodic orbit $\mathbf{x}_{\epsilon}(t)$ perturbing from
$\mathbf{x}=\mathbf{0}$ (see, e.g., Cabré et al. [10], Haller and Ponsioen
[21]). The manifold $\mathcal{W}_{\varepsilon}(E,t)$ is $O(\varepsilon)$
$C^{1}$-close to $\mathcal{W}_{0}(E,t)\equiv\mathcal{W}(E)$ and hence its
reduced dynamics can be parametrized using the reduced observable vector
$\boldsymbol{\varphi}=\boldsymbol{\phi}|_{\mathcal{W}(E)}$ in the form __

$\displaystyle\dot{\boldsymbol{\varphi}}$
$\displaystyle=\mathbf{B}\boldsymbol{\varphi}+\mathbf{q}\left(\boldsymbol{\varphi}\right)+\varepsilon\hat{\mathbf{F}}(\boldsymbol{\varphi},t),\qquad\mathbf{B}=D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\boldsymbol{\Lambda}_{E}\left(D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\right)^{-1},\quad\mathbf{q}\left(\boldsymbol{\varphi}\right)=\mathcal{O}\left(\left|\boldsymbol{\varphi}\right|^{2}\right),$
(61)
$\displaystyle\hat{\mathbf{F}}(\boldsymbol{\varphi},t)=\left(D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\right)^{-1}\left(\mathbf{I}+D\mathbf{h}\left(\boldsymbol{\varphi},\mathbf{0}\right)\right){}^{-1}\mathbf{F}(\mathbf{0},t)+\mathcal{O}\left(\varepsilon\left|\boldsymbol{\varphi}\right|^{2}\right),$

where we have relegated the details of this calculation to Appendix E.__

Then the unique, $C^{r}$ change of coordinates,

$\boldsymbol{\varphi}=\boldsymbol{\kappa}\left(\boldsymbol{\gamma}\right)=\boldsymbol{\gamma}+\boldsymbol{\ell}(\boldsymbol{\gamma}),$
(62)

guaranteed by statement (iii) of Theorem 4 transforms the reduced dynamics
(61) to its final form

$\dot{\boldsymbol{\gamma}}=\mathbf{B}\boldsymbol{\gamma}+\varepsilon(\mathbf{I}+D\boldsymbol{\ell}(\boldsymbol{\gamma}))^{-1}\hat{\mathbf{F}}(\mathbf{0},t).$
(63)

The transformation is valid on trajectories of (60) as long as they remain in
the domain of definition of the coordinate change (62).

Note that eq. (63) is a weakly perturbed, time-periodic nonlinear system. The
matrix $\mathbf{B}$ and the nonlinear terms
$\boldsymbol{\ell}(\boldsymbol{\gamma})$ can be determined using data from the
unforced ($\varepsilon=0$) system. As a result, nonlinear time-periodic forced
response can be predicted _solely from unforced data_ by applying numerical
continuation to system (63) for $\varepsilon>0$. This is not expected to be as
accurate as SSM-based forced response prediction (see, e.g., Cenedese et al.
[12, 13], Axås et al. [3], Axås and Haller [2]), but nevertheless offers a way
to make predictions for non-linearizable forced response based solely on DDL
performed on unforced data. These predictions are valid for forced
trajectories that stay in the domain of convergence of DDL carried out on the
unforced system. We will illustrate such predictions using actual experimental
data from fluid sloshing in Section 5.4.

Setting $\boldsymbol{\ell}(\boldsymbol{\gamma})=0$ in formula (63) enables us
to carry out a forced-response prediction based on DMD as well. Such a
prediction will be fundamentally linear with respect to the forcing and can
only be reasonably accurate for very small forcing amplitudes, as we will
indeed see in examples. There is no systematic way to model the addition of
non-autonomous forcing in the EDMD procedure, and hence EDMD will not be
included in our forced response prediction comparisons.

We also note, that one might be tempted to solve an approximate version of
(63) by assuming

$\varepsilon(\mathbf{I}+D\boldsymbol{\ell}(\boldsymbol{\gamma}))^{-1}\approx\varepsilon\mathbf{I}.$
(64)

This assumption simplifies the computation of the forced response of the
nonlinear system (63) to those of a simple linear system. Although the forced
response computed using this approximate DDL method turns out to be more
accurate than DMD on our example, we do not recommend this approach. This is
because neglecting the nonlinear effects of the coordinate change in (63) is,
in general, inconsistent with $\boldsymbol{\ell}(\boldsymbol{\gamma})\neq 0$.
We give more detail on this approximation in Appendix F of the Supplementary
Information.

##  5 Examples

In this section, we compare the DMD, EDMD and DDL algorithms on specific
examples. When applicable, we also compute the exact analytic linearization of
the dynamical system near its fixed point as a benchmark. On a slow SSM
$\mathcal{W}\left(E\right)$, an observer trajectory $\boldsymbol{\varphi}(t)$,
starting from a select initial condition $\boldsymbol{\varphi}(0)$, will be
tracked as the image of the linearized reduced observer trajectory
$\boldsymbol{\gamma}(t)$ under the linearizing transformation (62):

$\boldsymbol{\varphi}(t)=\boldsymbol{\kappa}\left(e^{\mathbf{B}t}\boldsymbol{\gamma}(0)\right)=e^{\mathbf{B}t}\boldsymbol{\gamma}(0)+\boldsymbol{\ell}\left(e^{\mathbf{B}t}\boldsymbol{\gamma}(0)\right),\qquad\boldsymbol{\gamma}(0)=\boldsymbol{\kappa}^{-1}\left(\boldsymbol{\varphi}(0)\right).$
(65)

When model reduction has also taken place, i.e., when the observable vector
$\boldsymbol{\varphi}$ is not defined on the full phase space, we will
nevertheless provide a prediction in the full phase space via the
parametrization of the slow SSM.

By Theorem 1, DMD can be interpreted as setting
$\boldsymbol{\ell}(\boldsymbol{\gamma})\equiv 0$ in (65) and finding the
linear operator $\mathbf{B}$ as a best fit from the available data. In
contrast, DDL finds the linear operator $\mathbf{B}$, the transformation
$\boldsymbol{\varphi}=\boldsymbol{\gamma}+\boldsymbol{\ell}(\boldsymbol{\gamma})$,
and its inverse simultaneously. As we explained in Section 4.2, EDMD cannot
quite be interpreted in terms of the linearizing transformation (65) as it is
an attempt to immerse the dynamics into a higher dimensional space. For our
EDMD tests, we will use monomials of the observable vector
$\boldsymbol{\varphi}$.

###  5.1 1D nonlinear system with two isolated fixed points

Consider the one-dimensional ODE obtained as the radial component of the
Stuart-Landau equation, i.e.,

$\dot{r}=\mu r-r^{3},$

which can be rescaled to

$\dot{R}=R-R^{3}.$ (66)

For $R\geq 0$, the system has a repelling fixed point at $R=0$ and an
attracting one at $R=1$. Page and Kerswell [44] show that local expansion of
observables in terms of the Koopman eigenfunctions computed near each fixed
point are possible, but the expansions at the two fixed points are not
compatible with each other and both diverge at $R=\sqrt{2}/2$. This is a
consequence of the more general result that the Koopman eigenfunctions
themselves inevitably blow up near basin boundaries (see our Proposition 1 in
Appendix A of the Supplementary Information). Both DMD and EDMD can
nevertheless be computed from data, even for a trajectory crossing the turning
point at $R=\sqrt{2}/2$, but the resulting models cannot have any connection
to the Koopman operator.

In each comparison performed on system (66), we generate a single trajectory
in the domain of attraction of the $R=1$ fixed point and use it as training
data for DMD, EDMD and DDL. In each subplot of Fig. 3, the single training
trajectory starts from the intersection of the red horizontal line “IC of
training trajectory” with the black curve labeled “True trajectory” at time
$t=0$. We then also generate a new test trajectory (black) with its initial
condition denoted with a black dot over the line $t=0$. We place this initial
condition slightly outside the domain of linearization for system (66) (under
the grey line labelled “Turning point”). We use DMD, order-$k=5$ EDMD, and DDL
trained on a single training trajectory to make predictions for the black
testing trajectory (not used in the training).

Figure 3a shows DDL to be the most accurate of the three methods when applied
to forward-time ($t\geq 0$) segments of the test trajectory. If we try to
predict the backward-time ($t<0$) segment of the same trajectory as it leaves
the training domain, DDL diverges immediately upwards, whereas DMD and EDMD
diverge more gradually downwards. As we increase the training domain in Fig
3b, DDL continues to be the most accurate in both forward and backward time
until it reaches the domain of its training range in backward time. At that
point, it diverges quickly upwards, while DMD and EDMD diverge more slowly
downwards.

Importantly, increasing the approximation order for DDL first to $k=10$ then
to $k=18$ (see Figs. 3c-d), makes DDL predictions more and more accurate in
backward time inside the training domain. At the same time, the same increase
in order makes EDMD less and less accurate inside the same domain. This is not
surprising for EDMD because it seeks to approximate the dynamics within a
Koopman-invariant subspaces for increasing $k$, and Koopman mode expansions
blow up at the “Turning point line”, as shown both analytically and
numerically by Page and Kerswell [44]. Interestingly, however, EDMD becomes
less accurate even within the domain of linearization under increasing $k$.
This is clearly visible in Fig. 3d which shows spurious, growing oscillations
in the EDMD predictions close to the $R=1$ fixed point.

In summary, of the three methods tested, DDL makes the most accurate
predictions in forward time. This remains true in backward time as longs as
the trajectory remains in the training range used for the three methods, even
if this range is larger than the theoretical domain of linearization. Inside
the training range, an increase of the order $k$ of the monomials used
increases the accuracy of DDL but introduces growing errors in EDMD.

Figure 3: Predictions of DMD, EDMD and DDL on trajectories of (66). (a)
Training trajectory starts inside the domain of convergence of the
linearization, i.e. $R(0)=0.8$ (see Page and Kerswell [44]). For both DDL and
EDMD the order of the monomials used is $k=5$. (b) Same for a different
training trajectory with $R(0)=0.4$ and $k=5$ (c) Same for $R(0)=0.1$ and
$k=10$ (d) Same for $R(0)=0.1$ and $k=18$.

###  5.2 3D linear system studied via nonlinear observables

Wu et al. [60] studied the ability of DMD to recover a 3D linear system based
on the time history of three nonlinear observables evaluated on the
trajectories of the system. To define the linear system, they use a block-
diagonal matrix $\boldsymbol{\Lambda}$ and a basis transformation matrix
$\mathbf{R}$ of the form

$\boldsymbol{\Lambda}=\left(\begin{array}[]{ccc}a&-b&0\\\ b&a&0\\\
0&0&c\end{array}\right),\quad
a,b,c\in\mathbb{R},\qquad\mathbf{R}=\left(\begin{array}[]{ccc}1&0&\sin\theta_{1}\cos\theta_{2}\\\
0&1&\sin\theta_{1}\sin\theta_{2}\\\ 0&0&\cos\theta_{2}\end{array}\right),$
(67)

to define the linear discrete dynamical system

$\mathbf{x}(n+1)=\left(\mathbf{R}\boldsymbol{\Lambda}\mathbf{R}^{-1}\right)\mathbf{x}(n).$
(68)

The linear change of coordinates $\mathbf{R}$ rotates the real eigenspace of
$\boldsymbol{\Lambda}$ corresponding to the eigenvalue $c$ and hence
introduces non-normality in system (68). This system is then assumed to be
observed via a 3D nonlinear observable vector

$\mathbf{y}(\mathbf{x})=\left(\begin{array}[]{c}x_{1}+0.1\left(x_{1}^{2}+x_{2}x_{3}\right)\\\
x_{2}+0.1\left(x_{2}^{2}+x_{1}x_{3}\right)\\\
x_{3}+0.1\left(x_{3}^{2}+x_{1}x_{2}\right)\end{array}\right).$ (69)

Ideally, DMD should closely approximate the linear dynamics of system (68)
because the observable function defined in eq. (69) is close to the identity
and has only weak nonlinearities. Wu et al. [60] find, however, that this
system poses a challenge for DMD, which produced inaccurate predictions for
the spectrum of $\mathbf{R}\boldsymbol{\Lambda}\mathbf{R}^{-1}$.

Following one of the parameter settings of Wu et al. [60], we set
$a=0.45\sqrt{3}$, $b=0.5$, $c=0.6$, $\theta_{1}=1.5$, and $\theta_{2}=0$. We
initialize three training trajectories with $\left\|\mathbf{x}(0)\right\|<1,$
each containing $100$ iterations of system (68). We then compute the
predictions of a $5^{th}$ order DDL model and compare to those of DMD and EDMD
on a separate test trajectory not used in training these three methods. The
predictions and the spectrum obtained from the three methods are shown in Fig.
4.

Figure 4: Predictions by DMD, EDMD, and DDL on the discrete dynamical system
(67),(69).(a) Predicted and true $y_{1}-$components of a test trajectory. (b)
Spectra identified by DMD, EDMD, and DDL superimposed on the true spectrum
(marked by crosses). The dashed line represents the unit circle.

The predictions of DMD and EDMD can only be considered accurate for very low
amplitude oscillations, while DDL returns accurate predictions throughout the
whole trajectory. This example consists of linear dynamics and monomial
observables of the state, and hence should be an ideal test case for EDMD.
Yet, EDMD is inaccurate in identifying the spectrum of system (68). Indeed, as
seen in Fig. 4b, a number of spurious eigenvalues arise from EDMD, both real
and complex, and only one can be considered to be close to a true eigenvalue.
DMD performs clearly better but it is still markedly less accurate than DDL.
These inaccuracies in the predictions of EDMD and DMD spectra are also
reflected by considerable errors in their predictions for trajectories, as
seen in Fig. 4a. In contrast, DDL produces the most accurate prediction for
the test trajectory.

###  5.3 Damped and periodically forced Duffing equation

We consider the damped and forced Duffing equation

$\displaystyle\dot{x}$ $\displaystyle=y,$ $\displaystyle\dot{y}$
$\displaystyle=x-x^{3}-dy+\varepsilon\cos\Omega t,$ (70)

with damping coefficient $d=0.0141$, forcing frequency $\Omega$ and forcing
amplitude $\varepsilon$. We perform a change of coordinates
$\left(x,y\right)\mapsto\boldsymbol{\varphi}=\left(\varphi_{1},\varphi_{2}\right)$
that moves the stable focus at $(x,y)=(1,0)$ to the origin and makes the
linear part block-diagonal. The resulting system is of the form

$\dot{\boldsymbol{\varphi}}=\mathbf{A}\boldsymbol{\varphi}+\mathbf{f}\left(\boldsymbol{\varphi}\right)+\varepsilon\hat{\mathbf{F}}(t),\qquad\mathbf{f}\left(\boldsymbol{\varphi}\right)=\mathcal{O}\left(\left|\boldsymbol{\varphi}\right|^{2}\right),$
(71)

where

$\mathbf{A}=\begin{pmatrix}-\alpha&-\omega\\\
\omega&-\alpha\end{pmatrix},\qquad\omega=1.4142,\quad\alpha=0.00707,$ (72)

and $\hat{\mathbf{F}}(t)$ is the transformed image of the physical forcing
vector in (70). We first consider the unforced system with $\varepsilon=0$. In
this case, the 2D slow SSM of the fixed point coincides with the phase space
$\mathbb{R}^{2}$ and hence no further model reduction is possible. However,
since the non-resonance conditions (36) hold for the linear part (72), the
system is analytically linearizable near the origin. The linearizing
transformation and its inverse can both be computed from eq. (71), as outlined
in eq. (43). For reference, we carry out this linearization analytically up to
order $k=9$. The Taylor series of the linearization is estimated to converge
for $\left|\boldsymbol{\varphi}\right|<R_{crit}\approx 0.15$. The details of
the calculation can be found in the repository [28].

We now compare the analytic linearization results it to DMD, EDMD and DDL,
with all three trained on the same three trajectories, launched both inside
and outside the domain of convergence of the analytic linearization. The
polynomial order of approximation is $k=5$ for both the EDMD and the DDL
algorithms. The performance of the various methods is compared in Fig. 5.
Close to the fixed point, in the domain of convergence of the analytic
linearization, all three methods perform well. Moving away from the fixed
point, the analytic linearization is no longer possible. Both DMD and EDMD
perform worse, while DDL continues to accurately linearize the system even
outside the domain of convergence of the analytic linearization.

Figure 5: Comparison of the time evolution of the linearized trajectories
(red) and the full trajectories of the nonlinear system (71) (blue), (71).
(a)-(d): Analytic linearization, DMD, EDMD and DDL models trained and
evaluated on trajectories inside the domain of convergence. (e)-(h): Same as
(a)-(d) but outside the domain of convergence of the analytic linearization.

Using formula (63) and our DDL-based model, we can also predict the response
of system (71) for the forcing term of the form

$\varepsilon\hat{\mathbf{F}}(t)=\varepsilon\left(\begin{array}[]{c}-0.006\\\
1.225\end{array}\right)\cos\Omega t,$ (73)

without using any data from the forced system. As the forced DDL model (63) is
nonlinear, it can capture non-linearizable phenomena such as coexisting of
stable and unstable periodic orbits arising under the forcing. We can also
make a forced response prediction from DMD simply by setting
$D\boldsymbol{\ell}(\boldsymbol{\gamma})=\boldsymbol{0}$ in eq. (63). As an
inhomogeneous linear system of ODEs, however, this forced DMD model cannot
predict coexisting stable and unstable periodic orbits.

In Fig. 6, we compare the forced predictions of the analytic linearization,
DMD, and DDL to those computed from the nonlinear system directly via the
continuation software COCO of Dankowicz and Schilder [14]. Since the forced
and linearized systems are also nonlinear, we use the same continuation
software to determine the stable and unstable branches of periodic orbits.

As expected, the analytic linearization is accurate while the forced response
is inside the domain of convergence but deteriorates quickly for larger
amplitudes. DMD gives good predictions for the peaks of the forced response
diagrams, but cannot account for any of the nonlinear softening behavior,
i.e., the overhangs in the curves that signal multiple coexisting periodic
responses at the same forcing frequency. In contrast, while the DDL model of
order $k=5$ starts becoming inaccurate for peak prediction at larger
amplitudes outside the domain of analytic linearization, it continues to
capture accurately the overhangs arising from non-linearizable forced response
away from the peaks. Notably, DDL even identifies the unstable branches (in
dashed lines) of the periodic response accurately. For completeness, we also
show results of approximate DDL, by assuming (64) in Appendix F.

Figure 6: Periodic response of the Duffing oscillator under the forcing (71).
The three distinctly colored forced response curves correspond to
$\varepsilon=0.001,0.002,0.0028$. (a) Actual nonlinear forced response from
numerical continuation (blue) and prediction for it from analytic
linearization (purple) (b) Forced response predictions from DDL (red) (c)
Forced response predictions from DMD (green). The training data for panels (b)
and (c) is the same unforced trajectory data set as the one used in Fig. 5.

###  5.4 Water sloshing experiment in a tank

In this section, we analyze experimental data generated by Bäuerlein and Avila
[5] for forced and unforced fluid sloshing in a tank. Previous studies of this
data set used nonlinear SSM-reduction to predict forced response ( Cenedese et
al. [12], Axås et al. [3], Axås and Haller [2]). Here we will use DMD and DDL
to extract and compare linear reduced-order models from unforced trajectory
data, then use them to predict and verify forced response curves obtained from
forced trajectory data. Neither DMD nor DDL is expected to outperform the
fully nonlinear approach of SSM reduction, so we will only compare them
against each other.

The tank in the experiments is mounted on a platform that is displaced
sinusoidally in time with various forcing amplitudes and frequencies (Fig.
7(a)). To train DMD and DDL, we use unforced sloshing data obtained by
freezing the movement of the tank near a resonance and recording the ensuing
decaying oscillations of the water surface with a camera under they die out.
The resulting videos serve as input data to our analysis. Specifically, the
horizontal position of the center of mass of the fluid is extracted from each
video frame tracked and used as the single scalar observable.

During such a resonance decay experiment, the system approaches its stable
unforced equilibrium via oscillations that are dominated by a single mode. In
terms of the phase space geometry, this means an approach to a stable fixed
point along its 2D slow SSM $\mathcal{W}\left(E\right)$ tangent to the slowest
2D real eigenspace $E$ . As we only have a single observable from the videos,
we use delay embedding to generate a larger observable space that can
accommodate the 2D manifold $\mathcal{W}\left(E\right)$. As discussed by
Cenedese et al. [12], we need an at least 5D observable space for this purpose
by the Takens embedding theorem. In this space, $\mathcal{W}\left(E\right)$
turns out to be nearly flat for short delays (see Axås and Haller [2]), which
allows us to use a linear approximation for its parametrization. The reduced
coordinates on $\mathcal{W}\left(E\right)\approx E$ can then be identified via
a singular value decomposition of the data after one removes initial
transients from the experimental data. The end of the transients can be
identified as a point beyond which a frequency analysis of the data shows only
one dominant frequency, the imaginary part of the eigenvalue corresponding to
$E$.

All this analysis has been carried out using the publicly available SSMLearn
package Cenedese et al. [11]. With $\mathcal{W}\left(E\right)$ identified, we
use the DDL method with order $k=5$ to find the linearizing transformation and
the linearized dynamics on $\mathcal{W}\left(E\right)$. In Fig 7(b) we show
the prediction of the model on a decaying trajectory reserved for testing. The
displacements are reported as percentage values, with respect to the depth of
the tank. In Fig. 7c-d, we show predictions from the DDL model for the forced
response, compared with the experimentally observed response. Since the exact
forcing function is unknown, we follow the calibration procedure outlined by
Cenedese et al. [12] to find an equivalent forcing amplitude in the reduced-
order model.

We present data for three forcing amplitudes. The DDL predictions are accurate
up to $0.17\%$ amplitude, even capturing the softening trend. The largest-
amplitude forcing resulted in response significantly outside the range of the
training data; in this range, we were unable to find the converged forced
response from DDL. We also show the corresponding DMD-predictions In Fig.
7c-d. Although the linear response can formally be evaluated for any forcing
amplitude, DMD shows no trace of the softening trend, and is even inaccurate
for low forcing amplitudes.

Figure 7: (a) Schematic representation of the experimental setup (adopted from
Cenedese et al. [12]). (b) Prediction of the decay of a test trajectory with
order-$k=5$ DDL. (c) Prediction of the forced response from DMD. (d)
Prediction of the forced response from DDL. Blue shading indicates the domain,
in which training data for DMD and DDL was available.

###  5.5 Model reduction and foliation in a nonlinear oscillator chain

As a final example, we consider the dynamics of a chain of nonlinear
oscillators, which has been analyzed in the SSMLearn package Cenedese et al.
[11]. Denoting the positions of the oscillators as $q_{i}$ for $i=1,...,5$, we
assume that the springs and dampers are linear, except for the first
oscillator. The non-dimensionalized equations of motion can be written as

$\mathbf{M\ddot{\mathbf{q}}}+\mathbf{C}\dot{\mathbf{q}}+\mathbf{K}\mathbf{q}+\mathbf{f}(\mathbf{q},\dot{\mathbf{q}})=\mathbf{0},$
(74)

where $\mathbf{M}=\mathbf{I}$; the springs have the same linear stiffness
$k=1$ which is encoded in $\mathbf{K}$ via nearest-neighbor coupling. The
damping is assumed to be proportional, i.e., we specifically set
$\mathbf{C}=0.002\mathbf{M}+0.005\mathbf{K}$.

Three numerically generated training trajectories show decay to the
$\mathbf{q}=\mathbf{0}$ fixed point, as expected from the damped nature of the
linear part of the system. In this example, we also seek to capture some of
the transients, which motivates us to select the slow SSM $\mathcal{W}(E)$ to
be 4D, tangent to the spectral subspace $E=E_{1}\oplus E_{2}$ spanned by the
the slowest mode ($E_{1})$ and the second slowest mode ($E_{2}$). As the mode
corresponding to $E_{2}$ does disappear over time from the decaying signal,
there is no resonance between the eigenvalues and hence Theorem 4 is
applicable. As already noted, numerical from a generic physical system
described by eq. (74) will be free from resonances. An exception is a $1\colon
1$ resonance arising from a perfect symmetry, but this resonance is not
excluded by Theorem 4 and has is amenable to DDL.

Within the 4D SSM $\mathcal{W}(E)$, we also demonstrate how to optimally
reduce the dynamics to its 2D slowest SSM $\mathcal{W}(E_{1})$. As explained
in Section 4.3.3, to find the trajectory in $\mathcal{W}(E_{1})$ with which a
given trajectory close to $\mathcal{W}(E)$ ultimately synchronizes, we need to
project along a point $\mathbf{q}_{0}$ of the full trajectory $\mathbf{q}(t)$
first onto $\mathcal{W}(E)$ orthogonally to obtain a point
$\mathbf{q}_{0}^{4D}\in\mathcal{W}(E)$. We then need to identify the stable
fiber $\mathcal{F}_{\mathbf{q}_{0}^{2D}}$ in $\mathcal{W}(E)$ for which
$\mathbf{q}_{0}^{4D}\in\mathcal{F}_{\mathbf{q}_{0}^{2D}}$ holds. Finally, one
has to project along $\mathcal{F}_{\mathbf{q}_{0}^{2D}}$ to locate its base
point $\mathbf{q}_{0}^{2D}\in\mathcal{W}(E_{1})$. The trajectory through
$\mathbf{q}_{0}^{2D}$ in $\mathcal{W}(E_{1})$ will then be the one with which
the full trajectory $\mathbf{q}(t)$ will synchronize faster than with any
other trajectory. As noted in Section 4.3.3, computing the full nonlinear
stable foliation

$\mathcal{W}(E)=\bigcup_{\mathbf{q}_{0}^{2D}\in\mathcal{W}(E_{1})}\mathcal{F}_{\mathbf{q}_{0}^{2D}}^{0}$
(75)

of $\mathcal{W}(E)$ is simple in the linearized coordinates, in which it can
be achieved via a linear projection along the faster eigenspace $E_{2}$.

We use a third-order polynomial approximation for $\mathcal{W}(E)$ based on
the three training trajectories. The polynomials depend on the reduced
coordinates we introduce along $E$ using a singular value decomposition of the
trajectory data. These reduced coordinates are shown in Fig. 8, where we show
a representative training trajectory, the 2D slow SSM $E_{1}$, as well as the
foliation (75) computed from DDL for this specific problem.

Figure 8: Reduced coordinates of the 4D SSM $\mathcal{W}(E)$ for the
oscillator chain. The slow 2D SSM $\mathcal{W}(E_{1})$, a typical trajectory
and its projection to the slow SSM along the fibers
$\mathcal{F}_{\mathbf{q}_{0}^{2D}}^{0}$ are also shown. Panel (a) shows the
linearized coordinates computed from DDL, and (b) shows their image under the
inverse of the linearizing transformation. The order of approximation used in
DDL is 3.

We also evaluate the DDL-based predictions on $\mathcal{W}(E)$ and
$\mathcal{W}(E_{1})$ by comparing them to predictions from DMD and EDMD.
Performing DMD and EDMD with the data first projected to $E$ can be
interpreted as finding the linear approximation to the dynamics in
$\mathcal{W}(E)$. Similarly, performing DMD and EDMD with the data first
projected to $E_{1}$ can be interpreted as finding the linear approximation to
the dynamics in $\mathcal{W}(E_{1})$. These are to be contrasted with
performing DDL that finds the linearized reduced dynamics within
$\mathcal{W}(E)$, which in turn contains the linearized reduced dynamics
within $\mathcal{W}(E_{1})$. Figure 9 shows that DMD and EDMD both perform
similarly to DDL on $\mathcal{W}(E)$. However, the 2D DMD and EDMD results
obtained for $\mathcal{W}(E_{1})$ are noticeably less accurate than the DDL
results.

Figure 9: Predictions of (a) DMD (b) EDMD and (c) DDL models on a test
trajectory of the oscillator chain. The order of approximation for DDL and
EDMD is $k=3$ .

##  6 Conclusions

We have given a new mathematical justification for the broadly used DMD
procedure to eliminate the shortcomings of prior proposed justifications.
Specifically, we have shown that under specific non-degeneracy conditions on
the $n$-dimensional dynamical system, on $d\leq n$ observable functions
defined for that system, and on the actual data from these observables, DMD
gives a leading-order approximation to the observable dynamics on an
attracting $d$-dimensional spectral submanifold (SSM) of the system.

This result covers both discrete and continuous dynamical systems even for
$n=\infty$. Our Theorem 1 only makes explicit non-degeneracy assumptions on
the observables which will hold with probability one in practical
applications. This is to be contrasted with prior approaches to DMD and its
variants based on the Koopman operator, whose assumptions on the observables
fail with probability one on generic observables.

Our approach also yields a systematic procedure that gradually refines the
leading-order DMD approximation of the reduced observable dynamics on SSMs to
higher orders. This procedure, which we call data-driven linearization (DDL),
builds a nonlinear coordinate transformation under which the observable
becomes linear on the attracting SSM. We have shown on several examples how
DDL indeed outperforms DMD and extended DMD (EDMD), as expected. In addition
to this performance increase, DDL also enables a prediction of truly nonlinear
forced response from unforced data within its training range. Although we have
only illustrated this for periodically forced water sloshing experiments in a
tank, recent results on aperiodically time-dependent SSMs by Haller and
Kaundinya [20] allow us to predict more general forced response using DDL
trained on unforced observable data.

Despite all these advantages, DDL (as any linearization method) remains
applicable only in parts of the phase space where the dynamics are
linearizable. Yet SSMs continue to exist across basin boundaries and hence are
able to carry characteristically nonlinear dynamics with multiple coexisting
attractors. For such nonlinearizable dynamics, data-driven nonlinear SSM-
reduction algorithms, such as SSMLearn and fastSSM, are preferable and have
been showing high accuracy and predictive ability in a growing number of
physical settings (see, e.g., Cenedese et al. [12, 13], Axås et al. [3],
Kaszás et al. [27], Alora et al. [1], Kaszás and Haller [26], Liu et al.
[37]).

Funding This work was supported by the Swiss National Science Foundation.

Data Availability All data and codes used in this work will be made available
after publication.

Competing Interests The authors declare that they have no conflict of
interest.

## References

  * Alora et al. [2023] J.I. Alora, M. Cenedese, E. Schmerling, G. Haller, and M. Pavone.  Practical deployment of spectral submanifold reduction for optimal control of high-dimensional systems.  _IFAC PapersOnLine_ , 56-2:4074–4081, 2023.
  * Axås and Haller [2023] J. Axås and G. Haller.  Model reduction for nonlinearizable dynamics via delay-embedded spectral submanifolds.  _Nonlinear Dyn._ , 2023.
  * Axås et al. [2022] J. Axås, M. Cenedese, and G. Haller.  Fast data-driven model reduction for nonlinear dynamical systems.  (https://doi.org/10.1007/s11071-022-08014-0), 2022.
  * Bates and Watts [1988] D. M. Bates and D. G. Watts.  _Nonlinear Regression Analysis and Its Applications_.  Wiley, 1988.
  * Bäuerlein and Avila [2021] B. Bäuerlein and K. Avila.  Phase lag predicts nonlinear response maxima in liquid-sloshing experiments.  _J. Fluid Mech._ , 925, 2021.
  * Bettini et al. [2024] L. Bettini, B. Kaszás, Bernhard Zybach, Jürg Dual, and George Haller.  Model reduction to spectral submanifolds via oblique projection.  _Preprint_ , 2024.
  * Bollt et al. [2018] E. M. Bollt, Q. Li, F. Dietrich, and I. Kevrekidis.  On matching, and even rectifying, dynamical systems through Koopman operator eigenfunctions.  _SIAM J. Appl. Dyn. Sys_ , 17(2):1925–1960, 2018\.
  * Budišić et al. [2012] M. Budišić, R. Mohr, and I. Mezić.  Applied Koopmanism.  _Chaos_ , 22:047510, 2012.
  * Buza [2024] G. Buza.  Spectral submanifolds of the Navier–Stokes equations.  _SIAM J. Appl. Dyn. Sys._ , 23(2):1052–1089, 2024.
  * Cabré et al. [2003] X. Cabré, E. Fontich, and R. de la Llave.  The parameterization method for invariant manifolds i: Manifolds associated to non-resonant subspaces.  _Indiana Univ. Math. J._ , 52(2):283–328, 2003\.
  * Cenedese et al. [2021] M. Cenedese, J. Axås, and G. Haller.  SSMLearn.  _https://github.com/haller-group/SSMLearn_ , 2021.
  * Cenedese et al. [2022a] 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.  _Nature Comm._ , 198, 2022a.
  * Cenedese et al. [2022b] M. Cenedese, J. Axås, H. Yang, M. Eriten, and G. Haller.  Data-driven nonlinear model reduction to spectral submanifolds in mechanical systems.  _Phil. Trans. R. Soc. A_ , 380:20210194, 2022b.
  * Dankowicz and Schilder [2013] H. Dankowicz and F. Schilder.  _Recipes for Continuation_.  SIAM, 2013.
  * Domb and Sykes [1957] C. Domb and M. F. Sykes.  On the susceptibility of a ferromagnetic above the curie point.  _Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences_ , 240(1221):214–228, 1957\.
  * Elbialy [2001] S.M. Elbialy.  Local contractions of Banach spaces and spectral gap conditions.  _J. Func. Anal._ , 182:108–150, 2001.
  * Eldering et al. [2018] J. Eldering, M Kvalheim, and S. Revzen.  Global linearization and fiber bundle structure of invariant manifolds.  _Nonlinearity_ , 31:4202–4245, 2018.
  * Fenichel [1971] N. Fenichel.  Persistence and smoothness of invariant manifolds for flows.  _Indiana Univ. Math. J._ , 21(3):193–226, 1971\.
  * Guckenheimer and Holmes [1983] J. Guckenheimer and P. Holmes.  _Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields_.  Springer, New York, 1983.
  * Haller and Kaundinya [2024] G. Haller and R. Kaundinya.  Nonlinear model reduction to temporally aperiodic spectral submanifolds.  _Chaos_ , Chaos (to appear);arXiv:2404.05355, 2024.
  * 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(3):1493–1534, 2016.
  * 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(6):063138, 2023.
  * Hartman [1960] P. Hartman.  On local homeomorphisms of Euclidean spaces.  _Bol. Soc. Mat. Mexicana_ , 5:220–241, 1960.
  * Jain et al. [2023] S. Jain, T. Thurner, M. Li, and G. Haller.  SSMTool 2.3: Computation of invariant manifolds and their reduced dynamics in high-dimensional mechanics problems.  _https://doi.org/10.5281/zenodo.4614201_ , pages 1417–1450, 2023\.
  * Jentzsch [1916] Robert Jentzsch.  Untersuchungen zur theorie der folgen analytischer funktionen.  _Acta Mathematica_ , 41:219–251, 1916.
  * Kaszás and Haller [2024] B. Kaszás and G. Haller.  Capturing the edge of chaos as a spectral submanifold in pipe flows.  _J. Fluid. Mech._ , 979:A48, 2024.
  * Kaszás et al. [2022] B. Kaszás, M. Cenedese, and G. Haller.  Dynamics-based machine learning of transitions in couette flow.  _Phys. Rev. Fluids._ , 7:L082402, 2022.
  * Kaszás et al. [2024] B. Kaszás, , and G. Haller.  Data-Driven Linearization: Numerical Implementation.  _[Link to the repository will be inserted after publication]_ , 2024.
  * Kingma and Ba [2015] Diederik Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In _International Conference on Learning Representations (ICLR)_ , 2015.
  * Kogelbauer and Haller [2018] F. Kogelbauer and G. Haller.  Rigorous model reduction for a damped-forced nonlinear beam model: An infinite-dimensional analysis.  _J. Nonlin. Sci._ , 28:1109–1150, 2018.
  * Koopman [1931] B.O. Koopman.  Hamiltonian systems and transformation Hilbert space.  _Proc.Natl. Acad. Sci. USA_ , 7:315–318, 1931.
  * Kutz et al. [2016] J.N. Kutz, S.L. Brunton, B.W. Brunton, and J.L. Proctor.  _Dynamic Mode Decomposition_.  SIAM, Philadelphia, PA, 2016.
  * Kvalheim and Arathoon [2023] M.D. Kvalheim and P. Arathoon.  Linearizability of flows by embeddings.  _arXiv:2305.18288_ , 2023.
  * Kvalheim and Revzen [2021] M.D. Kvalheim and S. Revzen.  Existence and uniqueness of global koopman eigenfunctions for stable fixed points and periodic orbits.  _Physica D_ , 425:132959, 2021.
  * Lan and Mezic [2013] Y. Lan and I. Mezic.  Linearization in the large of nonlinear systems and Koopman operator spectrum.  _Physica D_ , 242:42–53, 2013.
  * Li et al. [2017] Q. Li, F. Dietrich, E. M. Bollt, and I. G. Kevrekidis.  Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the Koopman operator.  _Chaos_ , 27(10):103111, 10 2017.
  * Liu et al. [2022] A. Liu, J. Axås, and G. Haller.  Data-driven modeling and forecasting of chaotic dynamics on inertial manifolds constructed as spectral submanifolds.  _Chaos_ , 34:033140, 2022.
  * Liu et al. [2023] Z. Liu, N. Ozay, and E.D. Sontag.  On the non-existence of immersions for systems with multiple omega-limit sets.  _IFAC-PapersOnLine_ , 56(2):60–64, 2023.
  * Liu et al. [2024] Z. Liu, N. Ozay, and E.D. Sontag.  Properties of immersions for systems with multiple limit sets with implications to learning Koopman embeddings.  _arXiv:2312.17045_ , pages 1–14, 2024.
  * Mezić [2013] I. Mezić.  Analysis of fluid flows via spectral properties of the Koopman operator.  _Ann. Rev. Fluid Mech._ , 45(1):357–378, 2013\.
  * Newhouse [2017] S. E. Newhouse.  On a differentiable linearization theorem of Philip Hartman.  _Contemporary. Math._ , 692:209–262, 2017.
  * Otto et al. [2022] S. E. Otto, A. Padovan, and C. W. Rowley.  Optimizing oblique projections for nonlinear systems using trajectories.  _SIAM J. Sci. Comput._ , 44(3):A1681–A1702, 2022.
  * Otto et al. [2023] S. E. Otto, A. Padovan, and C. W. Rowley.  Model reduction for nonlinear systems by balanced truncation of state and gradient covariance.  _SIAM J. Sci. Comput._ , 45(5):A2325–A2355, 2023.
  * Page and Kerswell [2019] J. Page and R.R. Kerswell.  Koopman mode expansions between simple invariant solutions.  _J. Fluid Mech._ , 879:1–27, 2019.
  * Paszke et al. [2017] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer.  Automatic differentiation in pytorch.  In _NIPS-W_ , 2017.
  * Poincaré [1892] H. Poincaré.  _Les Méthodes Nouvelles de la Mécanique Céleste_.  Gauthier-Villars et Fils, Paris, 1892.
  * Ponsioen et al. [2018] S. Ponsioen, T. Pedergnana, and G. Haller.  Automated computation of autonomous spectral submanifolds for nonlinear modal analysis.  _J. Sound Vibr._ , 420:269–295, 2018.
  * Rowley [2005] C. W. Rowley.  Model reduction for fluids, using balanced proper orthogonal decomposition.  _Int. J. Bifurcat. Chaos_ , 15(03):997–1013, 2005\.
  * Rowley et al. [2009] C.W. Rowley, I. Mezić, S. Bagheri, P. Schlachter, and D.S. Henningson.  Spectral analysis of nonlinear flows.  _J. Fluid Mech._ , 641:115–127, 2009.
  * Schmid [2010] P.J. Schmid.  Dynamic mode decomposition of numerical and experimental data.  _J. Fluid Mech._ , 656:5–28, 2010.
  * Schmid [2022] P.J. Schmid.  Dynamic mode decomposition and its variants.  _Ann. Rev. Fluid Mech._ , 54:225–254, 2022.
  * Sternberg [1957] S. Sternberg.  Local contractions and theorem of Poincaré.  _Am. J. Math._ , 79(4):809–824, 1957.
  * Sternberg [1958] S. Sternberg.  On the structure of local homeomorphisms of Euclidean n-space, II.  _Am. J. Math._ , 80(3):623–631, 1958.
  * Szaksz [2024] B. Szaksz.  _The Stabilizing and Destabilizing Effects of Time Delays in Nonlinear Dynamical Systems_.  Ph.D. Thesis, Budapest University of Technology and Economics, 2024\.
  * Szalai [2020] R. Szalai.  Invariant spectral foliations with applications to model order reduction and synthesis.  _Nonlin. Dyn._ , 101:2645–2669, 2020.
  * van Strien [1990] S. van Strien.  Smooth linearization of hyperbolic fixed points without resonance conditions.  _J. Diff. Eqs._ , 85(1):66–90, 1990.
  * Virtanen et al. [2020] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17:261–272, 2020.
  * Williams et al. [2015a] M.O. Williams, I.G. Kevrekidis, and C.W. Rowley.  A data–driven approximation of the Koopman operator: Extending dynamic mode decomposition.  _J. Nonlin. Sci._ , 9:1307–1346, 2015a.
  * Williams et al. [2015b] M.O. Williams, C.W. Rowley, and I.G. Kevrekidis.  A kernel-based method for data-driven Koopman spectral analysis.  _J. Comp. Dyn._ , 2:247–265, 2015b.
  * Wu et al. [2021] Z. Wu, S. L. Brunton, and S. Revzen.  Challenges in dynamic mode decomposition.  _J. Royal Soc. Interface_ , 18(185):20210686, 2021.

Supplementary Information for

_Data-Driven Linearization of Dynamical Systems_

George Haller444Corresponding author. Email<EMAIL_ADDRESS>and Bálint
Kaszás
Institute for Mechanical Systems
ETH Zürich
Leonhardstrasse 21, 8092 Zürich, Switzerland

##  Appendix A DMD and the Koopman operator

###  A.1 Observable dynamics and the Koopman operator

A classic approach, introduced first by Koopman [31] and revived recently by
multiple authors (see, e.g., by Budišić et al. [8], Mezić [40], Kutz et al.
[32] and the references cited therein), describes observable evolution via the
Koopman operator
$\boldsymbol{\mathcal{K}}^{t}\colon\mathcal{G}^{d}\to\mathcal{G}^{d}$, defined
on a Banach space $\mathcal{G}^{d}$ of $d$-dimensional observable functions as
the pull-back operation on observables under the flow map $\mathbf{F}^{t}$ of
system (1). Specifically,

$\boldsymbol{\mathcal{K}}^{t}\left[\boldsymbol{\phi}\right]\,(\mathbf{x}_{0}):=\boldsymbol{\phi}\left(\mathbf{F}^{t}(\mathbf{x}_{0})\right),$
(76)

or, in more compact notation,

$\boldsymbol{\mathcal{K}}^{t}\left[\boldsymbol{\phi}\right]=\boldsymbol{\phi}\circ\mathbf{F}^{t},\quad\boldsymbol{\phi}\in\mathcal{G}^{d}.$
(77)

As $\mathcal{G}^{d}$ is a complete vector space, linear combinations of $\ell$
observables $\phi_{1},\ldots,\phi_{\ell}\in\mathcal{G}^{d}$ are also
observables in $\mathcal{G}^{d}$ and hence $\boldsymbol{\mathcal{K}}^{t}$ can
be evaluated on them. Such an evaluation gives

$\displaystyle\boldsymbol{\mathcal{K}}^{t}\left(c_{1}\phi_{1}+\ldots+c_{1}\phi_{\ell}\right)(\mathbf{x}_{0})$
$\displaystyle=c_{1}\boldsymbol{\mathcal{K}}^{t}\left[\phi_{1}\right](\mathbf{x}_{0})+\ldots+c_{\ell}\boldsymbol{\mathcal{K}}^{t}\left[\phi_{\ell}\right](\mathbf{x}_{0})$
(78)

_This last equation shows that the mapping $\boldsymbol{\mathcal{K}}^{t}$ is
linear with respect to the choice of observables while the system’s initial
condition $\mathbf{x}_{0}$ is kept fixed. In contrast, DMD seeks to find an
approximate mapping that is linear with respect to the choice of the system’s
initial conditions while the observable is kept fixed._ Indeed, in
applications of DMD, eq. (6) is posed for a fixed observable vector evaluated
on various initial conditions. EDMD is also trained on system trajectories
launched from different initial conditions, observed under the same general
function of an initially fixed set of observables.

On a broader note, the linearity of the Koopman operator with respect to
changes in the choice of an observable, as seen in eq. (78), does not imply
linear dynamics for observations of the system under a fixed observable. For
instance, let $\phi_{1}(\mathbf{x})$ denote the observation of the
displacement, measured in meters, of the end point of a nonlinear beam whose
current state in its phase space is denoted by $\mathbf{x}$. If, instead, we
want to observe the position of the endpoint in millimeters, then we switch
the observable $\phi_{2}(\mathbf{x})=10^{3}\phi_{1}(\mathbf{x}),$ and hence
for the observed initial state $\mathbf{x}_{0}$ of the beam, we have the
change

$\phi_{2}(\mathbf{x}_{0})=c\phi_{1}(\mathbf{x}_{0}),\qquad c=10^{3},$

in the initial observation. It is then no surprise that this change in the
initial observable will transform to the same change in the current
observations of the endpoint at time $t$, i.e.,

$\boldsymbol{\mathcal{K}}^{t}\left[\phi_{2}\right](\mathbf{x}_{0})=c\boldsymbol{\mathcal{K}}^{t}\left[\phi_{1}\right](\mathbf{x}_{0}).$

This simple fact, however, just reflects that the final results at time $t$
change by the same multiple as their initial conditions if we decide to change
the physical units in our measurements. Clearly, this fact does not imply any
linear dynamics for the displacement of the endpoint of the beam. Indeed,
$\boldsymbol{\mathcal{K}}^{t}\left[\phi_{2}\right](\mathbf{x}_{0})$ may well
be a chaotic signal.

Another example would be tracking the total population of a continent and
denoting the current state of the population by $\mathbf{x}$. Let
$\phi_{1}(\mathbf{x})$ denote specifically the observable returning the
population of country $A$ and $\phi_{2}(\mathbf{x})$ denote the observable
returning the population of country $B$ within the same continent. If we
decide now to jointly observe the population of these two countries, then we
are in effect passing to a third observable at at initial state
$\mathbf{x}_{0}$ to obtain

$\phi_{3}(\mathbf{x}_{0})=\phi_{2}(b_{0})+\phi_{1}(\mathbf{x}_{0}),$

for the initial observations. To obtain the total population of countries $A$
and $B$ at the current state $\mathbf{x}$ of the total population, we can
simply write

$\boldsymbol{\mathcal{K}}^{t}\left[\phi_{3}\right](\mathbf{x}_{0})=\boldsymbol{\mathcal{K}}^{t}\left[\phi_{1}\right](\mathbf{x}_{0})+\boldsymbol{\mathcal{K}}^{t}\left[\phi_{2}\right](\mathbf{x}_{0}).$

Again, this fact by itself does not say anything about the linearity or
nonlinearity of the population growth in either country $A$ or country $B$.

Finally, _the Koopman operator is as much nonlinear as it is linear with
respect to the changes in the observables_ as it commutes with any nonlinear
single or multi-variable function of observables, as longs as that function is
also an observable. Indeed, for any function

$\mathbf{g}\colon\left(\mathcal{G}^{d}\right)^{\ell}\to\mathcal{G}^{d},$

the observable
$\boldsymbol{\psi}(\mathbf{x}_{0})=\mathbf{g}\left(\phi_{1}(\mathbf{x}_{0}),\ldots,\phi_{\ell}(\mathbf{x}_{0})\right)$
satisfies

$\boldsymbol{\mathcal{K}}^{t}\left[\mathbf{g}\left(\phi_{1},\ldots,\phi_{\ell}\right)\right]=\mathbf{g}\left(\boldsymbol{\mathcal{K}}^{t}\left[\phi_{1}\right],\ldots,\boldsymbol{\mathcal{K}}^{t}\left[\phi_{\ell}\right]\right).$

To illustrate this general property of the Koopman operator, let $\mathbf{x}$
denote the instantaneous velocity field of a 2D Navier-Stokes flow with
density $\rho$ and let $\phi_{ij}(\mathbf{x})=\left|v_{ij}\right|$ be the
specific velocity magnitude along the $(i,j)$ location of a fixed spatial
grid. If we now want to track the local kinetic energy of the flow at $(i,j)$,
we do not have to reformulate and re-solve the Navier–Stokes equation for
$\phi_{\mathrm{kin}}(\mathbf{x})=\frac{1}{2}\rho\phi_{ij}^{2}(\mathbf{x})=\left|x_{ij}\right|^{2}$.
Rather, we can use our already available observations of
$\phi_{ij}(\mathbf{x})$ of the velocity field starting from an initial
velocity field $\mathbf{x}_{0}$ to obtain
$\boldsymbol{\mathcal{K}}^{t}\left[\phi_{\mathrm{kin}}\right](\mathbf{x}_{0})=\frac{1}{2}\rho\left[\boldsymbol{\mathcal{K}}^{t}\left[\phi_{ij}\right](\mathbf{x}_{0})\right]^{2}.$
This, however, does not imply a statement that “the dynamics of the 2D
Navier–Stokes equation are quadratic in the space of observables”.

###  A.2 Differential equation for the dynamics of observables under the
Koopman operator

We can directly verify that

$\boldsymbol{\mathcal{K}}^{t_{1}+t_{2}}\boldsymbol{\phi}=\mathcal{\boldsymbol{\mathcal{K}}}^{t_{2}}\boldsymbol{\mathcal{K}}^{t_{1}}\boldsymbol{\phi}=\boldsymbol{\mathcal{K}}^{t_{1}}\mathcal{\boldsymbol{\mathcal{K}}}^{t_{2}}\boldsymbol{\phi},\quad\left(\boldsymbol{\mathcal{K}}^{t}\right)^{-1}\boldsymbol{\phi}=\boldsymbol{\mathcal{K}}^{-t}\boldsymbol{\phi},$

i.e., $\boldsymbol{\mathcal{K}}^{t}$ defines a flow map on the observable
space $\mathcal{G}^{d}$, as long as system (1) generates a flow
$\mathbf{F}^{t}$ on the phase space $\mathcal{P}$. To describe the dynamics
generated by this flow map, let us select an arbitrary scalar observable
$\boldsymbol{\phi}\in\mathcal{G}^{d}$ and differentiate the defining relation
(76) with respect to $t$ to obtain __

$\frac{d}{dt}\left\\{\text{$\mathcal{\boldsymbol{\mathcal{K}}}^{t}\left[\boldsymbol{\phi}\right]$}(\boldsymbol{x}_{0})\right\\}=D\boldsymbol{\phi}\left(\mathbf{F}^{t}(x_{0})\right)\mathbf{f}\left(\mathbf{F}^{t}(x_{0})\right).$
(79)

Since we have