0.8pt 11affiliationtext: Department of Mathematics, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA 22affiliationtext: Department of
Mechanical Engineering, University of Chile, Beauchef 851, Santiago, Chile
33affiliationtext: Graduate Aerospace Laboratories, California Institute of
Technology, Pasadena CA 91125, USA 44affiliationtext: Department of Applied
Mathematics, University of Washington, Seattle, WA 98195, USA
55affiliationtext: Department of Mechanical Engineering, University of
Washington, Seattle, WA 98195, USA

# Physics-informed dynamic mode decomposition (piDMD)

Peter J. Baddoo∗ Benjamin Herrmann  Beverley J. McKeon  J. Nathan Kutz  Steven
L. Brunton

###### Abstract

In this work, we demonstrate how physical principles – such as symmetries,
invariances, and conservation laws – can be integrated into the dynamic mode
decomposition (DMD). DMD is a widely-used data analysis technique that
extracts low-rank modal structures and dynamics from high-dimensional
measurements. However, DMD frequently produces models that are sensitive to
noise, fail to generalize outside the training data, and violate basic
physical laws. Our physics-informed DMD (piDMD) optimization, which may be
formulated as a Procrustes problem, restricts the family of admissible models
to a matrix manifold that respects the physical structure of the system. We
focus on five fundamental physical principles – conservation, self-
adjointness, localization, causality, and shift-invariance – and derive
several closed-form solutions and efficient algorithms for the corresponding
piDMD optimizations. With fewer degrees of freedom, piDMD models are less
prone to overfitting, require less training data, and are often less
computationally expensive to build than standard DMD models. We demonstrate
piDMD on a range of challenging problems in the physical sciences, including
energy-preserving fluid flow, travelling-wave systems, the Schrödinger
equation, solute advection-diffusion, a system with causal dynamics, and
three-dimensional transitional channel flow. In each case, piDMD significantly
outperforms standard DMD in metrics such as spectral identification, state
prediction, and estimation of optimal forcings and responses.

††∗ Corresponding author (baddoo@mit.edu).

##  1 Introduction

Integrating partial knowledge of physical principles into data-driven
techniques is a primary goal of the scientific machine learning (ML) community
[1]. Physical principles – such as conservation laws, symmetries, and
invariances – can be incorporated into ML algorithms in the form of inductive
biases, thereby ensuring that the learned models are constrained to the
correct physics. Recent successful examples of ML algorithms that have been
modified to respect physical principles include neural networks [2, 3, 4, 5,
6, 7, 8, 9, 10], kernel methods [11, 12], deep generative models [13], and
sparse regression [14, 15, 16, 17, 18]. These examples demonstrate that
incorporating partially-known physical principles into machine learning
architectures can increase the accuracy, robustness, and generalizability of
the resulting models, while simultaneously decreasing the required training
data. In this work, we integrate knowledge of physical principles into one of
the most widely-used methods in data-driven dynamical systems research: the
dynamic mode decomposition [19, 20, 21, 22, 23, 24].

The dynamic mode decomposition (DMD) is a data diagnostic technique that
extracts coherent spatial-temporal patterns from high-dimensional time series
data [19, 24]. Although DMD originated in the fluid dynamics community [19],
the algorithm has since been applied to a wealth of dynamical systems
including in epidemiology [25], robotics [26, 27], neuroscience [28], quantum
control [29], power grids [30], and plasma physics [31, 32]. Despite its
widespread successes, DMD is highly sensitive to noise [33, 34, 35], fails to
capture travelling wave physics, and can produce overfit models that do not
generalize. Herein, we demonstrate that integrating physics into the learning
framework can help address these challenges.

Suppose that we are studying a dynamical system defined by
$\dot{\boldsymbol{x}}=\boldsymbol{F}(\boldsymbol{x})$ (continuous time) or
$\boldsymbol{x}_{k+1}=\boldsymbol{F}(\boldsymbol{x}_{k})$ (discrete time)
where $\boldsymbol{F}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is unknown. DMD
identifies the best low-rank linear approximation of $\boldsymbol{F}$ given a
collection of $m$ pairs of measurements
$\\{\boldsymbol{x}_{j},\boldsymbol{y}_{j}\\}$; in other words, DMD seeks a
rank $r$ matrix $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ such that

$\displaystyle\boldsymbol{y}_{j}\approx\boldsymbol{A}\boldsymbol{x}_{j}$ (1)

for $j=1,\dots,m$. Arranging the data measurements into $n\times m$ matrices
$\boldsymbol{X}$ and $\boldsymbol{Y}$ allows us to phrase the above formally
as

$\displaystyle\operatorname*{arg\rm{}min}_{\textrm{rank
}(\boldsymbol{A})=r}\|\boldsymbol{Y}-\boldsymbol{A}\boldsymbol{X}\|_{F}.\qquad\qquad\qquad$
(2)

After approximately solving (2), the DMD process computes the dominant
spectral properties of the learned linear operator [19, 20, 21, 22, 23, 24].
The rank-$r$ constraint in (2) is motivated by the assumed modal structure of
the system but doesn’t account for other important physical properties. For
example, one limitation of DMD is that the solution of (2) lies within the
span of $\boldsymbol{Y}$, so the learned model rarely generalizes outside the
training regime. We address this limitation and others by embedding partial
physics knowledge into the learning process.

In this work, we incorporate physical principles into the optimization (2) by
constraining the solution matrix $\boldsymbol{A}$ to lie on a matrix manifold
$\mathcal{M}\subseteq\mathbb{R}^{n\times n}$:

piDMD regression:whitepiDMD regression:
$\displaystyle\operatorname*{arg\rm{}min}_{\boldsymbol{A}\in\mathcal{M}}\|\boldsymbol{Y}-\boldsymbol{A}\boldsymbol{X}\|_{F}.\qquad\qquad\quad\qquad$
(3)

Figure 1: Comparing standard dynamic mode decomposition (specifically,
optimized DMD [23]) and physics-informed dynamic mode decomposition applied to
the advection equation. The data is contaminated with 2% additive Gaussian
noise. Having trained DMD models, we then perform predictions for different
initial conditions. Standard DMD fails whereas piDMD produces a faithful
prediction. In the figure, $\boldsymbol{\mathcal{S}}_{\phi}$ represents the
shift operator (see §4.1).

The matrix manifold $\mathcal{M}$ is dictated by the known physics of the
system at hand. For example, we may select $\mathcal{M}$ such that its members
satisfy certain symmetries that we known are obeyed by the system at hand. We
call (3) ‘ _physics-informed_ DMD’ (piDMD) as the optimization integrates
underlying knowledge of the system physics into the learning framework. Again,
the spectral features of the solution of (3) can be computed to give insight
into the modes that dominate the system and their dynamics. The low-rank DMD
regression (2) is a special case of piDMD where $\mathcal{M}$ is the manifold
of rank-$r$ matrices. However, piDMD models are not always low rank (see, for
example, §4.4) but nevertheless have few degrees of freedom due to the
physics-informed constraints. Constraining the solution matrix improves the
generalizability of the learned model, reduces sensitivity to noise and
reduces the demand for large training sets. More broadly, constraining a
linear, dimensionality-reduced model to a manifold $\mathcal{M}$ has been
shown to result in a variety of optimization formulations and techniques which
can improve low-rank approximations [36]. In figure 1 we compare the
performance of standard (2) and physics-informed (3) DMD for a simple
travelling wave. We constrain the piDMD model to the shift-invariant symmetry
of the underlying system. As a result, the piDMD model produces an accurate
prediction of the future dynamics, whereas the standard DMD model blows up
immediately. This simple example illustrates that embedding basic physics into
DMD can substantially improve the algorithm’s performance.

The optimization problem (3) is well-known in the statistical literature as a
_Procrustes problem_ [37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49]. Our
recasting of the DMD regression as a Procrustes problem is a new connection,
and is the basis of piDMD. This perspective enables us to leverage the
substantial extant work on Procrustes problems into new application areas. The
literature contains many exact solutions for Procrustes problems, including
notable cases of orthogonal matrices [45], and symmetric matrices [41].

A significant challenge in physics-informed ML architectures is to develop
implementations that incorporate known physics but also scale well to higher
dimensions [1]. A major contribution of this work is several exact solutions,
in terms of standard linear algebra operations, of Procrustes problems of
physical relevance. Where possible, we have made use of matrix factorizations
and low-rank representations to alleviate the computational burden of
implementation; it is rarely necessary to form $\boldsymbol{A}$ explicitly
when computing its spectral properties or other diagnostics. For complicated
matrix manifolds, an exact solution of the Procrustes problem (3) may be
intractable and we must resort to algorithmic approaches. Fortunately,
optimization on matrix manifolds is a mature field [50], and many algorithms
are implemented in the open source software package ‘Manopt’ [51,
www.manopt.org].

The remainder of the paper is arranged as follows. In section 2 we provide
background information on DMD and Procrustes problems. Then, in section 3, we
describe the broad framework of piDMD. We consider a range of applications in
section 4, with a focus on shift-invariant, conservative, self-adjoint, local
and causal systems. Section 5 concludes with a discussion of the limitations
of piDMD and suggests several future research directions. An open source
implementation of piDMD with support for over 30 matrix manifolds is available
in Matlab at www.github.com/baddoo/piDMD.

##  2 Mathematical background

In this section we provide further details on DMD and Procrustes problems.
Throughout the article, we assume that we have access to $m$ snapshots pairs
of $n$ features each:
$\\{(\boldsymbol{x}_{j},\,\boldsymbol{y}_{j}),\,j=1,\dots,m\\}$. For example,
$\boldsymbol{x}_{j}$ may be a discretized fluid flow field at time $t_{j}$ and
$\boldsymbol{y}_{j}$ may be the flow field at the next time step $t_{j+1}$. It
is convenient to arrange the data into $n\times m$ snapshot matrices of the
form

$\displaystyle\boldsymbol{X}$ $\displaystyle=\begin{bmatrix}|&|&|\\\
\boldsymbol{x}_{1}&\cdots&\boldsymbol{x}_{m}\\\
|&|&|\end{bmatrix}=\begin{bmatrix}\rule[1.93748pt]{10.79999pt}{0.55pt}&\tilde{\boldsymbol{x}}_{1}&\rule[1.93748pt]{10.79999pt}{0.55pt}\\\\[-3.87498pt]
\rule[3.09999pt]{10.79999pt}{0.55pt}&{\vdots}&\rule[3.09999pt]{10.79999pt}{0.55pt}\\\
\rule[1.93748pt]{10.79999pt}{0.55pt}&\tilde{\boldsymbol{x}}_{n}&\rule[1.93748pt]{10.79999pt}{0.55pt}\\\
\end{bmatrix},\qquad\boldsymbol{Y}=\begin{bmatrix}|&|&|\\\
\boldsymbol{y}_{1}&\cdots&\boldsymbol{y}_{m}\\\
|&|&|\end{bmatrix}=\begin{bmatrix}\rule[1.93748pt]{10.79999pt}{0.55pt}&\tilde{\boldsymbol{y}}_{1}&\rule[1.93748pt]{10.79999pt}{0.55pt}\\\\[-3.87498pt]
\rule[3.09999pt]{10.79999pt}{0.55pt}&{\vdots}&\rule[3.09999pt]{10.79999pt}{0.55pt}\\\
\rule[1.93748pt]{10.79999pt}{0.55pt}&\tilde{\boldsymbol{y}}_{n}&\rule[1.93748pt]{10.79999pt}{0.55pt}\\\
\end{bmatrix}$ (4)

so that rows $i$ of $\boldsymbol{X}$ and $\boldsymbol{Y}$ are the measurements
of the $i$-th features and columns $j$ are the $j$-th temporal snapshots.
Henceforth, we use $\tilde{\cdot}$ to represent a row vector.

###  2.1 Dynamic mode decomposition

DMD was initially proposed as a dimensionality reduction technique that
extracts dominant spatio-temporal coherent structures from high-dimensional
time-series data [19]. In particular, DMD identifies the leading-order spatial
eigenmodes of the matrix $\boldsymbol{A}$ in (1), along with a linear model
for how the amplitudes of these coherent structures evolve in time. DMD has
been applied to a range of systems, as summarized in the monograph [22] and
the recent review by Schmid [24].

To address the challenges associated with DMD, researchers have derived many
variations on the original algorithm, including sparsity promoting DMD [52],
DMD with control [53], noise-robust variants [33, 34, 35, 23], recursive DMD
[54], online DMD [55, 56], and versions for under-resolved data in space or
time [57, 58, 59]. At present, the most widely used variant is ‘exact DMD’
[21], which phrases the DMD solution in terms of the Moore–Penrose
pseudoinverse. This solution produces a more even distribution of errors
between the terms. Due to its simplicity and widespread use, most of the
comparisons in this paper are made between exact DMD and piDMD.

It is assumed that each $\boldsymbol{x}_{j}$ and $\boldsymbol{y}_{j}$ are
connected by an unknown dynamical system of the form
$\boldsymbol{y}_{j}=\boldsymbol{F}(\boldsymbol{x}_{j})$; for discrete-time
dynamics $\boldsymbol{y}_{j}=\boldsymbol{x}_{j+1}$ and for continuous-time
dynamics $\boldsymbol{y}_{j}=\dot{\boldsymbol{x}}_{j}$. DMD aims to learn the
dominant behaviour of $\boldsymbol{F}$ by finding the best linear
approximation for $\boldsymbol{F}$ given the data and then performing
diagnostics on that approximation. Thus, DMD seeks the linear operator
$\boldsymbol{A}$ that best maps the snapshots in the set
$\\{\boldsymbol{x}_{j}\\}$ to those in the set $\\{\boldsymbol{y}_{j}\\}$:

$\displaystyle\boldsymbol{y}_{j}\approx\boldsymbol{A}\boldsymbol{x}_{j}\qquad\qquad\textrm{for
}j=1,\dots,m.$ (5)

Expressed in terms of the snapshot matrices in (4), the linear system in (5)
becomes

$\displaystyle{\boldsymbol{Y}\approx\boldsymbol{A}\boldsymbol{X}},$ (6)

and the optimization problem for $\boldsymbol{A}$ is given by (2). The
minimum-norm solution for $\boldsymbol{A}$ is given by

$\boldsymbol{A}=\boldsymbol{Y}\boldsymbol{X}^{\dagger}=\boldsymbol{Y}\boldsymbol{V}\boldsymbol{\Sigma}^{\dagger}\boldsymbol{U}^{\ast},$
(7)

where $\dagger$ indicates the Moore–Penrose pseudoinverse [60] and
$\boldsymbol{X}=\boldsymbol{U\Sigma V}^{\ast}$ is the singular value
decomposition. In many applications, the state dimension $n$ is very large,
and forming or storing $\boldsymbol{A}$ explicitly becomes impractical.
Instead, we use a rank-$r$ approximation for $\boldsymbol{A}$, denoted by
$\hat{\boldsymbol{A}}$, where $r\ll n$. To form $\hat{\boldsymbol{A}}$, we
construct the optimal rank-$r$ approximation for $\boldsymbol{X}$ using the
truncated singular value decomposition [61]:
$\boldsymbol{X}\approx\boldsymbol{U}_{r}\boldsymbol{\Sigma}_{r}\boldsymbol{V}_{r}^{\ast}$.
We then project $\boldsymbol{A}$ onto the leading $r$ principal components of
$\boldsymbol{X}$ as

$\hat{\boldsymbol{A}}=\boldsymbol{U}_{r}^{\ast}\boldsymbol{A}\boldsymbol{U}_{r}=\boldsymbol{U}_{r}^{\ast}\boldsymbol{Y}\boldsymbol{V}_{r}\boldsymbol{\Sigma}_{r}^{-1}.$
(8)

It is now computationally viable to compute the eigendecomposition of
$\hat{\boldsymbol{A}}$ as

$\hat{\boldsymbol{A}}\hat{\boldsymbol{\Psi}}=\hat{\boldsymbol{\Psi}}\boldsymbol{\Lambda}.$
(9)

The eigenvectors of $\boldsymbol{A}$ can be approximated from the reduced
eigenvectors $\boldsymbol{\Psi}$ by [21]

$\boldsymbol{\Psi}=\boldsymbol{Y}\boldsymbol{V}\boldsymbol{\Sigma}^{-1}\hat{\boldsymbol{\Psi}}.$
(10)

This eigendecomposition is connected to the Koopman operator of the system,
and allows reconstructions and predictions [20, 62, 22, 63]. For example, for
a discrete-time system (i.e., $\boldsymbol{y}_{k}=\boldsymbol{x}_{k+1}$) with
evenly spaced samples in time, then the eigenvectors form a linearly
independent set and

$\displaystyle\boldsymbol{x}_{j}=\boldsymbol{\Psi}\boldsymbol{\Lambda}^{j-1}\boldsymbol{b}$
(11)

where the vector $\boldsymbol{b}$ contains the weights of the modes in the
initial condition:
$\boldsymbol{b}=\boldsymbol{\Psi}^{\dagger}\boldsymbol{x}_{1}$. From the
above, it is clear that the eigenvalues $\boldsymbol{\Lambda}$ govern the
temporal behaviour of the system and the eigenvectors $\boldsymbol{\Psi}$ are
the spatial modes.

###  2.2 Procrustes problems

Procrustes problems (3) comprise of finding the optimal transformation between
two matrices subject to certain constraints on the class of admissible
transformations [45, 40, 41]. According to Greek mythology, Procrustes was a
bandit who would stretch or amputate the limbs of his victims to force them to
fit onto his bed. Herein, $\boldsymbol{X}$ plays the role of Procrustes’
victim, $\boldsymbol{Y}$ is the bed, and $\boldsymbol{A}$ is the ‘treatment’
(stretching or amputation).111This terminology was first introduced by Hurley
and Cattell [37]. Procrustes problems (3) seek to learn the treatment
$\boldsymbol{A}$ that best represents the data measurements $\boldsymbol{X}$
and $\boldsymbol{Y}$. The minimization is usually phrased in the Frobenius
norm. Procrustes analysis finds relevance in a wealth of fields including
sensory analysis [64], data visualization [65], neural networks [66], climate
science [67], and solid mechanics [41]. A summary is available in the
monograph [40].

The earliest and most common formulation is the ‘orthogonal Procrustes
problem’ [45, 40]. Suppose that we have two sets of measurements
$\boldsymbol{X}$ and $\boldsymbol{Y}$ that we know are related by an unitary
(orthogonal) transformation (i.e., a rotation or reflection). The goal is to
learn the best unitary transformation that relates the data measurements.
Thus, $\boldsymbol{A}$ is constrained to be a unitary matrix and the
minimization problem is

$\displaystyle\operatorname*{arg\rm{}min}_{\boldsymbol{A}^{\ast}\boldsymbol{A}=\boldsymbol{I}}\|\boldsymbol{Y}-\boldsymbol{A}\boldsymbol{X}\|_{F}.$
(12)

The solution to (12) was derived by Schönemann [45] as

$\displaystyle\boldsymbol{A}=\boldsymbol{U}_{YX}\boldsymbol{V}_{YX}^{\ast}$
(13)

where
$\boldsymbol{U}_{YX}\boldsymbol{\Sigma}_{YX}\boldsymbol{V}_{YX}^{\ast}=\boldsymbol{Y}\boldsymbol{X}^{\ast}$
is a full singular value decomposition. Alternatively,
$\boldsymbol{A}=\boldsymbol{U}_{P}$, where
$\boldsymbol{U}_{P}\boldsymbol{H}_{P}=\boldsymbol{Y}\boldsymbol{X}^{\ast}$ is
a polar decomposition. The solution is unique if and only if
$\boldsymbol{Y}\boldsymbol{X}^{\ast}$ is full-rank. As we show in section 4.2,
a unitary matrix corresponds to an energy-preserving operator.

There are many solutions for Procrustes problems available in the literature
for different matrix constraints [41]. When exact solutions are not possible,
algorithmic solutions can be effective [51].

##  3 Physics-informed dynamic mode decomposition (piDMD)

Incorporating physics into ML algorithms involves supplementing an existing
technique with additional biases. Usually, ML practitioners take one of three
approaches [1]. First, observational biases can be embedded through data
augmentation techniques; however, augmented data is not always available, and
incorporating additional training samples can become computationally
expensive. Second, physical constraints can be included by suitably penalising
the loss function [8]. Third, inductive biases can be incorporated directly
into the machine learning architecture in the form of mathematical constraints
[16]. As noted by [1], this third approach is arguably the most principled
method since it produces model that strictly satisfy the physical constraints.
piDMD falls into the final category.

piDMD incorporates physical principles by constraining the matrix manifold of
the DMD regression problem (3). We abandon the low-rank representation of
$\boldsymbol{A}$ typically sought in model order reduction in favour of
alternative or additional matrix structures more relevant to the physical
problem at hand. In figure 2 we illustrate the matrix structures used in this
paper, along with their corresponding physical principle and references for
the optimal solutions.

Figure 2: Visual illustrations of the types of matrices, their corresponding
physical principles and references to the solutions of the corresponding
optimization problem (3).

Analysing a system with the piDMD framework requires four steps: modeling,
interpretation, optimization, and diagnostics:

  1. 1.

ModelingwhiteModeling: When studying a system from a data-driven perspective,
it is common to have some partial knowledge of the underlying physics of the
system. This first step in piDMD asks the user to summarise the known or
suspected physical properties of the system at hand.

  2. 2.

Interpretation:whiteInterpretation: Having determined the physical principle
we would like to enforce, we must translate these laws into the matrix
manifold to which the linear model should be constrained. This step involves
knowledge of the data collection procedure such as the spatial grid or the
energy inner product.

  3. 3.

Optimization:whiteOptimization: Equipped with a target matrix manifold, we may
now solve the relevant Procrustes problem (3). By suitably modifying the
optimization routine, we can guarantee that the resulting model satisfies the
physical principle identified in step 1.

  4. 4.

Diagnostics:whiteDiagnostics: The final step of piDMD involves extracting
physical information from the learned model $\boldsymbol{A}$. For example, one
may wish to analyse the spectrum or modes, compute the resolvent modes [68],
perform predictions, or investigate other diagnostics.

As an example, consider the travelling wave solution explored in the
introduction (figure 1). We began with the physical principle that the system
is shift invariant (‘modeling’), and this lead us to seek a circulant matrix
model (‘interpretation’). We found a solution to the corresponding Procrustes
problem (‘optimization’) and investigated the predictive behaviour of the
model (‘diagnostics’).

While each of the four steps can be problematic, the optimization step is
arguably the most conceptually and computationally difficult. In particular,
finding a solution for the desired Procrustes problem can be challenging.
Moreover, any solutions are quite peculiar to the matrix manifold under
consideration: the symmetric Procrustes problem is largely unconnected to,
say, the tridiagonal Procrustes problem. Thus, in this article, we present
many new solutions for (3) with different physics-informed matrix manifold
constraints. When exact solutions are not possible (e.g., if the manifold
constraint is quite complicated) there are sophisticated algorithmic solutions
available [51].

Standard DMD exploits the low-rank structure of $\boldsymbol{A}$ to
efficiently perform diagnostics on the learned model. Some of the manifolds we
consider (such as circulant, tridiagonal, upper triangular) do not have an
obvious or useful low-rank approximation. Instead, these matrices often have
an alternative structure that can be exploited to perform fast diagnostic
tests. For example, tridiagonal matrices rarely have a meaningful low-rank
approximation, but nevertheless admit a fast eigendecomposition and singular
value decomposition [60].

Some of the matrix manifolds we consider (such as symmetric, triangular,
tridiagonal or circulant) can be phrased as linear equality constraints and
can, in principle, be solved with linear-equality constrained least squares
[60]. However, the number of equality constraints needed is
$\mathcal{O}(n^{2})$ so the resulting least squares matrix will have
$\mathcal{O}(n^{2})$ rows, which is intractably large in most applications.
Thus, we avoid phrasing the piDMD constraints in terms of linear equality
constraints and instead exploit properties of the matrix manifold to find
solutions that can be efficiently implemented (see sections 4.3 and 4.1, for
example).

##  4 Applying piDMD to enforce canonical physical structures

We now present a range of examples of piDMD. Each of the following sections
investigates a detailed application of piDMD for a specific physical
principle. A summary is illustrated in figure 3. The physical principles,
corresponding matrix structures, and optimal solutions and extensions are
summarised below:

  1. §4.1:

Shift-invariant:whiteShift-invariant: circulant (and symmetric, skew-
symmetric, or unitary, §A.1), low rank (§A.2), non-equally spaced samples
(§A.3), total least squares (§A.4), Toeplitz & Hankel (§A.5)

  2. §4.2:

Conservative:whiteConservative: unitary (§2.2)

  3. §4.3:

Self-adjoint:whiteSelf-adjoint: symmetric (§4.3), skew-symmetric (§B.1),
symmetric in a subspace (§4.3.1)

  4. §4.4:

Local:whiteLocal: tridiagonal (§4.4), variable diagonals (§C.1), periodic
(§C.2), symmetric tridiagonal (§C.3), total least squares (§C.4), regularized
locality (§C.5).

  5. §4.5:

Causal:whiteCausal: upper triangular (§D.1 D.2)

Figure 3:  A comparison of the models learned by exact DMD and piDMD for a
range of applications. piDMD identifies the spectrum of the true operator with
higher accuracy than exact DMD. The structure of the model matrices are also
illustrated; piDMD models are generally more coherent than those learned by
exact DMD. Details are given in the corresponding sections. In the spectrum
subplots, we plot the true eigenvalues as , the DMD eigenvalues as  , and the
piDMD eigenvalues as .

###  4.1 Shift-invariant systems

Shift invariance (also called ‘translation invariance’) is ubiquitous in the
physical sciences. Spatially homogeneous systems – such as constant
coefficient PDEs or convolution operators – are shift invariant, and thus
appear identical with respect to any (stationary) observer. In view of
Noether’s theorem [69], the conserved quantity related to shift invariance is
linear momentum; as such, incorporating shift invariance into a DMD model
means that the model preserves linear momentum. While it may be unusual for
practical engineering problems to have a truly homogeneous direction, many
problems of basic physical interest – such as wall-bounded turbulent flow –
frequently possess a spatially homogeneous dimension. In the sequel, we derive
a piDMD formulation for shift-invariant operators.

We define $\mathcal{S}_{\phi}$ as the $\phi$-shift operator i.e.
$\mathcal{S}_{\phi}v(\xi)=v(\xi+\phi)$ for all test functions $v$. We say that
a space-continuous linear operator $\mathcal{A}$ is shift invariant if
$\mathcal{A}$ commutes with the $\phi$-shift operator for all shifts $\phi$:

$\displaystyle\mathcal{S}_{\phi}\mathcal{A}=\mathcal{A}\mathcal{S}_{\phi}.$
(14)

An application of the shift operator shows that if $\mathcal{A}$ is shift
invariant then the quantity
$\textrm{e}^{-\lambda\xi}\mathcal{A}\textrm{e}^{\lambda\xi}$ is constant for
all $\xi$. Thus, neglecting boundary conditions for now,
$\\{\textrm{e}^{\lambda\xi}\\}$ are eigenfunctions of $\mathcal{A}$. If the
domain is normalised to $[-1,1]$ and has periodic boundaries, then
$\lambda=l\pi i$ for integer $l$. In this case the corresponding
eigenfunctions form an orthogonal basis and the operator is diagonalised by
the known eigenfunctions.

In the following, we assume that we are studying a shift-invariant system on a
periodic domain $\xi\in[-1,1]$. We move from the continuous formulation to a
discretized space and assume that we are studying a function $u(\xi,t)$ and
have access to evenly spaced samples of $u$ at
$\boldsymbol{\xi}=[-1,\,-1+\Delta\xi,\dots,1-\Delta\xi]$. We discuss the case
of non-equally spaced samples in appendix A.3. If we define the state variable
as
$\boldsymbol{x}(t)=\left[u(\xi_{1},t),\,u(\xi_{2},t),\cdots\,,u(\xi_{n},t)\right]$
then the discrete-space linear operator $\boldsymbol{A}$ that generates
$\boldsymbol{x}(t)$ is diagonalized by its eigenvectors:

$\displaystyle\boldsymbol{A}=\boldsymbol{\mathcal{F}}\textrm{diag}\left(\hat{\boldsymbol{a}}\right)\boldsymbol{\mathcal{F}}^{-1}$
(15)

where
$\boldsymbol{\mathcal{F}}_{j,k}=\textrm{e}^{2\pi\textrm{i}(j-1)(k-1)/n}/\sqrt{n}$
and $\boldsymbol{\mathcal{F}}^{-1}=\boldsymbol{\mathcal{F}}^{\ast}$ and
$\\{\hat{a}_{j}\\}$ are the unknown eigenvalues. Equation (15) is equivalent
to stating that $\boldsymbol{A}$ is circulant:

$\displaystyle\boldsymbol{A}=\begin{bmatrix}a_{0}&a_{n-1}&\dots&a_{1}\\\
a_{1}&a_{0}&\ddots&\vdots\\\ \vdots&\ddots&\ddots&a_{n-1}\\\
a_{n-1}&\dots&a_{1}&a_{0}\end{bmatrix},$ (16)

i.e. $\boldsymbol{A}_{j,k}=a_{(j-k)\,\textrm{mod}\,n}$. If the boundary
conditions are instead Dirichlet or Neumann then $\boldsymbol{A}$ is not
circulant but Toeplitz; we solve the corresponding Procrustes problem in §A.5.

Substituting (15) into (3) and noting that the Frobenius norm is invariant to
unitary transformations allows (3) to be transformed to

$\displaystyle\operatorname*{arg\rm{}min}_{\hat{\boldsymbol{a}}}\left\|\textrm{diag}(\hat{\boldsymbol{a}})\boldsymbol{\mathcal{X}}-\boldsymbol{\mathcal{Y}}\right\|_{F},$
(17)

where $\boldsymbol{\mathcal{X}}=\boldsymbol{\mathcal{F}}^{\ast}\boldsymbol{X}$
and $\boldsymbol{\mathcal{Y}}=\boldsymbol{\mathcal{F}}^{\ast}\boldsymbol{Y}$
are the spatial discrete Fourier transforms of $\boldsymbol{X}$ and
$\boldsymbol{Y}$. As such, $\boldsymbol{\mathcal{X}}$ and
$\boldsymbol{\mathcal{Y}}$ can be efficiently formed in $O(mn\log(n))$
operations using Fast Fourier Transform (FFT) methods [70].

The rows of the cost function (17) now decouple to produce $n$ minimization
problems:

$\displaystyle\operatorname*{arg\rm{}min}_{\hat{a}_{j}}\|\hat{a}_{j}\tilde{\boldsymbol{\mathcal{X}}}_{j}-\tilde{\boldsymbol{\mathcal{Y}}}_{j}\|_{F}\qquad\textnormal{for
}1\leq j\leq n,$ (18)

where $\tilde{\boldsymbol{\mathcal{X}}}_{j}$ and
$\tilde{\boldsymbol{\mathcal{Y}}}_{j}$ are the $j$th rows of
$\boldsymbol{\mathcal{X}}$ and $\boldsymbol{\mathcal{Y}}$, respectively. An
alternative view of this step is that spatial wavenumbers decouple in shift-
invariant systems, so we may analyse each wavenumber individually. The optimal
solution for each eigenvalue follows as

$\displaystyle\hat{a}_{j}=\tilde{\boldsymbol{\mathcal{Y}}}_{j}\tilde{\boldsymbol{\mathcal{X}}}_{j}^{\dagger}=\left.{\tilde{\boldsymbol{\mathcal{Y}}}_{j}\tilde{\boldsymbol{\mathcal{X}}}_{j}^{\ast}}\middle/{\|\tilde{\boldsymbol{\mathcal{X}}}_{j}\|^{2}_{2}}\right..$
(19)

In principle, standard DMD methods [19, 21] could be encouraged to respect
shift invariance by augmenting the data matrices with their shifted
counterparts. However, one would strictly need to include $n$ shifts, thus
producing $nm$ samples which is usually too large in applications. In
contrast, the above solution is extremely efficient via the use of the FFT and
the decoupling of wavenumbers.

The appendices include alternative solutions for cases where the system is
shift-invariant and symmetric, skew-symmetric or unitary (§A.1), low rank
(§A.2), the samples are not equally spaced (§A.3), the total least squares
case (§A.4), and when the system is Toeplitz or Hankel (§A.5).

####  4.1.1 Example: plane channel flow

Our next example is the incompressible flow inside a plane channel of size
$2\pi\times 2\times 2\pi$ along the spatial $x$, $y$, and $z$ coordinates that
indicate the streamwise, wall-normal, and spanwise directions, respectively.
The configuration considers a Reynolds number of $Re=2000$ based on the
channel half-height and the centerline velocity, and periodic boundary
conditions in the open faces of the channel, hence the flow is homogeneous in
the $x$ and $z$ directions. We use the spectral code _Channelflow_ [71] to
perform direct numerical simulations (DNS) with the same numerical
configurations as that presented in [68].

The dataset investigated is generated from a DNS of the response of the
laminar flow with parabolic velocity profile to a localized perturbation in
the wall-normal velocity component of the form

$v(x,y,z,0)=\left(1-\frac{r^{2}}{c_{r}^{2}}\right)\left(\cos(\pi
y)+1\right)e^{\left(-r^{2}/c_{r}^{2}-y^{2}/c_{y}^{2}\right)},$ (20)

where $r^{2}=(x-\pi)^{2}+(z-\pi)^{2}$, the parameters are set to $c_{r}=0.7$
and $c_{y}=0.6$, and the amplitude of the perturbation was scaled to have an
energy-norm of $10^{-5}$ to ensure that the effect of nonlinearity is
negligible. This initial condition, which was first studied in [72] and also
in [68], is a model of a disturbance that could be generated in experiments
using a spanwise and streamwise periodic array of axisymmetric jets injecting
fluid perpendicular to the wall. The dataset obtained from this simulation is
then comprised of a sequence of $300$ snapshots of the three-dimensional
velocity-perturbation field recorded every $0.5$ time units.

Figure 4: Learning the spectrum of the linearized Navier–Stokes equations from
velocity measurements of the response to a localized disturbance in a channel
flow. piDMD embeds the shift-invariant structure of the Navier–Stokes
equations into the model learning framework, thus learning the spectrum of the
linearized operator with improved accuracy over exact DMD.

Our aim is to learn the spectrum of the underlying linear operator that _best_
describes the evolution of the DNS snapshots. When directly applying exact DMD
to this dataset, the algorithm attempts to find global modes for the three-
dimensional velocity field and their spectrum. This is challenging because
modes associated to different streamwise and spanwise wavenumbers may be
mixed, leading to spurious eigenvalues, as shown in figure 4. However, we know
that the flow is homogeneous in the $x$ and $z$ coordinates, and is therefore
shift-invariant in these directions. Hence, we can leverage piDMD to
incorporate this property, forcing the resulting matrix to respect, by
construction, the structure with three nested levels shown in figure 4. In
practice, this amounts to reshaping every snapshot, taking the FFT in the
dimensions corresponding $x$ and $z$, and performing DMD on the $y$-dependent
Fourier amplitudes for every streamwise and spanwise wavenumber tuple.
Although simple, this has a tremendous impact on the quality of the learned
spectrum. Shift-invariant piDMD produces modes that are associated with a
single wavenumber tuple, which forces the modes to respect the spatial
periodicity of the underlying system.

The spectra learned from both approaches are compared in figure 4, showing
that piDMD results in a far more accurate eigenvalue spectrum. The DMD
calculation considers a truncated SVD of rank $200$. The piDMD results were
obtained with a rank truncation of $50$ at every wavenumber tuple.

###  4.2 Conservative systems

Conservation laws are foundational to all areas of science. Identifying
quantities that remain constant in time – such as mass, momentum, energy,
electric charge, and probability – allow us to derive universal governing
equations, apply inductive reasoning, and understand basic physical
mechanisms. In this section we demonstrate how conservation laws can be
incorporated into the DMD framework.

Suppose that we are studying a system that we know conserves energy. In
applications of DMD, it is implicitly assumed that measurements of the state
have been suitably weighted so that the square of the 2-norm corresponds to
the energy of the state [68]: $E(\boldsymbol{x})=\|\boldsymbol{x}\|_{2}^{2}$.
In these variables, the original optimization problems (2, 3) equate to
finding the model $\boldsymbol{A}$ that minimises the energy of the error
between the true and predicted states ($\boldsymbol{y}_{k}$ and
$\boldsymbol{A}\boldsymbol{x}_{k}$ respectively). Thus, if $\boldsymbol{A}$
represents a discrete-time linear dynamical system
($\boldsymbol{y}_{k}=\boldsymbol{x}_{k+1}=\boldsymbol{A}\boldsymbol{x}_{k}$)
then $\boldsymbol{A}$ is an energy preserving operator if and only if

$\displaystyle
E(\boldsymbol{A}\boldsymbol{x})=\|\boldsymbol{A}\boldsymbol{x}\|_{2}^{2}=\|\boldsymbol{x}\|_{2}^{2}=E(\boldsymbol{x})\qquad\textnormal{for
all }\boldsymbol{x}\in\mathbb{R}^{n}.$ (21)

In words, (21) states that $\boldsymbol{A}$ does not change the energy of the
system but merely redistributes energy between the states. In matrix
terminology, (21) holds if and only if $\boldsymbol{A}$ is _unitary_.
Therefore, for conservative systems the optimization problem (3) is the
orthogonal Procrustes problem described in section 2.2, and the solution is
given by (13). The eigenvalues of $\boldsymbol{A}$ lie on the unit circle, and
the eigenvectors are orthogonal. Thus, the eigenvectors oscillate in time,
with no growth or decay. Since the solution of the orthogonal Procrustes
problem (13) requires the full SVD of an $n\times n$ matrix, it can be more
computationally efficient to first project onto the leading POD modes and
build a model therein. In this case, the model is only energy preserving
within the subspace spanned by the leading POD modes.

Noise is a substantial problem for most DMD methods, and a common remedy is to
phrase the DMD optimization (2) as a total least squares problem [34]. Perhaps
surprisingly, the solution to the orthogonal Procrustes problem (13) is also
the solution to the total least squares problem when the solution is
constrained to be orthogonal [43]. Thus, the solution (13) is optimal even
when there is noise in both $\boldsymbol{X}$ and $\boldsymbol{Y}$, as long as
the noise has the same distribution in $\boldsymbol{X}$ and $\boldsymbol{Y}$.
There are various problems closely related to the orthogonal Procrustes
problem, including the case where $\boldsymbol{A}$ lies on the Steifel
manifold (i.e., when $\boldsymbol{A}$ is rectangular with orthonormal columns
[48]), there is a weighting matrix [47], missing data [42], or high amounts of
noise [46].

####  4.2.1 Example: $Re=100$ flow past a cylinder

Figure 5: Learning the spectrum of the flow past a cylinder with noisy
measurements with DMD and energy-preserving piDMD.

We now apply this energy-preserving piDMD to study the flow past a cylinder,
which is a benchmark problem for modal decompositions [22]. The Reynolds
number is $Re=DU/\nu=100$ where $D$ is the diameter of the cylinder, $U$ is
the free-stream velocity and $\nu$ is the kinematic viscosity. The data
consist of 151 samples of vorticity measurements, corresponding to five
periods of vortex shedding, taken on a grid of 199$\times$449 points. We note
that the 2-norm of the measurements is approximately constant in time. The
data are contaminated with 20% Gaussian noise, as illustrated in the top left
panel of figure 5. We truncate the data to the first 15 POD modes and learn a
standard DMD model, and an energy-preserving piDMD model.

piDMD learns a more accurate representation of the leading-order dynamics than
standard DMD. For the high-order modes, the eigenvalues learned by DMD exhibit
spurious damping whereas the eigenvalues learned by piDMD are purely
oscillatory and remain on the imaginary axis. The trajectories of the POD
coefficients are also more accurate for piDMD, and do not exhibit the spurious
energy loss caused by noise.

###  4.3 Self-adjoint systems

Self-adjoint systems are another important class of linear systems that arise
frequently in solid mechanics, heat and mass transfer, fluid mechanics, and
quantum mechanics. When studying a system known to be self-adjoint we can use
piDMD to constrain our model $\boldsymbol{A}$ to lie in the manifold of
symmetric (Hermitian) matrices, such that
$\boldsymbol{A}=\boldsymbol{A}^{\ast}$. Symmetric matrices have real
eigenvalues and, by the spectral theorem, are diagonalisable. This restriction
places a significant constraint on the learned model and can substantially
reduce the sensitivity of DMD to noise.

The minimum-norm solution of the symmetric Procrustes problem [73] is

$\displaystyle\boldsymbol{A}=\boldsymbol{U}_{X}\,\boldsymbol{L}\,\boldsymbol{U}_{X}^{\ast},$
(22)

where the entries of $\boldsymbol{L}$ are

$\displaystyle\boldsymbol{L}_{i,j}=\overline{\boldsymbol{L}_{j,i}}=\begin{cases}\dfrac{\sigma_{i}\overline{\boldsymbol{C}_{j,i}}+\sigma_{j}\boldsymbol{C}_{i,j}}{\sigma_{i}^{2}+\sigma_{j}^{2}}&\textnormal{if
}\sigma_{i}^{2}+\sigma_{j}^{2}\neq 0,\\\ \qquad\qquad
0&\textnormal{otherwise},\end{cases}$ (23)

where
$\boldsymbol{C}=\boldsymbol{U}_{X}^{\ast}\boldsymbol{Y}\boldsymbol{V}_{X}$ and
$\boldsymbol{X}=\boldsymbol{U}_{X}\boldsymbol{\Sigma}_{X}\boldsymbol{V}_{X}^{\ast}$
is the SVD of $\boldsymbol{X}$, and $\sigma_{i}$ is the $i$th singular value.
The solution of the skew-symmetric case is similar (§B.1), and algorithmic
solutions are available when $\boldsymbol{A}$ is instead constrained to be
positive definite [44]. For large scale systems, a low-rank approximation to
$\boldsymbol{A}$ can be obtained by truncating the SVD of $\boldsymbol{X}$. If
the SVD is truncated to rank $r$, then $\boldsymbol{L}$ is an $r\times r$
matrix. This low-rank approximation preserves the self-adjointness of the
model.

The symmetric Procrustes problem is connected to the orthogonal Procrustes
problem (§4.2) via the fact that every unitary matrix $\boldsymbol{U}$ can be
expressed as the exponential of a Hermitian matrix $\boldsymbol{H}$ as
$\boldsymbol{U}=\exp(i\boldsymbol{H})$.

####  4.3.1 Learning energy states of the Schrödinger equation

Possibly the most famous self-adjoint operator is the quantum Hamiltonian,
$\hat{H}$ [74]. In quantum mechanics, physical quantities of interest, such as
position, momentum, energy, and spin, are represented by self-adjoint linear
operators called ‘observables’, of which one example is the Hamiltonian. The
Hamiltonian describes the evolution of the probabilistic wave function via the
time-dependent Schrödinger equation:

$\displaystyle i\hbar{\frac{d}{dt}}|\Psi(t)\rangle={\hat{H}}|\Psi(t)\rangle.$
(24)

Solutions of (24) take the form of an eigenmode expansion

$\displaystyle\Psi(\boldsymbol{\xi},t)=\sum^{\infty}_{j=1}\alpha_{j}\,\psi_{j}(\boldsymbol{\xi})\,\textrm{e}^{-\textrm{i}E_{j}t/\hbar},$
(25)

where $E_{j}$ is an energy level (eigenvalue) of $\hat{H}$ with corresponding
eigenstate $\psi_{j}$, the coefficients $\alpha_{j}$ are determined by the
initial distribution of the wave function, and $\boldsymbol{\xi}$ is the
spatial coordinate.

Figure 6: Learning the eigenstates of a quantum Hamiltonian of a hexagonal
well with the self-adjoint piDMD.

We consider a wave function evolving according to the two-dimensional
Schrödinger’s equation subject to an unknown Hamiltonian function. The unknown
potential function is taken to be a finite well in a hexagonal shape. We
collect measurements of the wave function and its velocity and train an exact
DMD model and a piDMD model. In practice, empirical measurements of the
wavefunction can be obtained through optical tomographic methods [75, 76] or
sequential measurements of complementary variables [77]. The data matrices
$\boldsymbol{X}$ and $\boldsymbol{Y}$ are formed from measurements of
$|\Psi(t)\rangle$ and $i\hbar\frac{d}{dt}|\Psi(t)\rangle$ respectively, which
are themselves formed by superposing eigenfunctions obtained by a finite
difference method. Thus, the matrices $\boldsymbol{X}$ and $\boldsymbol{Y}$
are connected by an unknown self-adjoint Hamiltonian $\hat{H}$ and we may
apply our self-adjoint DMD optimization. The measurements of $\boldsymbol{Y}$
are contaminated with Gaussian noise but the measurements of $\boldsymbol{X}$
are clean.

Constraining the learned Hamiltonian to be self-adjoint reduces the
sensitivity of the learning process to noise, as evidenced by figure 6. From
the bottom left panel, we see that the energy levels of the piDMD model are
physically consistent insofar as they are all real and positive. In contrast,
the energy levels learned by standard DMD exhibit a spurious imaginary
component, which will produce unrealistic growth/decay of the eigenstates.
Both methods miss one eigenstate ($E_{n}\sim 7.5$) since this particular mode
is not well-represented in the data. Additionally, the eigenstates learned by
piDMD are much less noisy than those learned by standard DMD. This can be
explained by noting that the eigenstates learned by piDMD lie in the span of
$\boldsymbol{X}$, which is clean, whereas the eigenstates learned by standard
DMD lie in the span of $\boldsymbol{Y}$, which is noisy. Lemma B.1 proves that
the symmetric piDMD model (22) is less sensitive to noise than the exact DMD
model (7). In summary, a physics-aware approach enables a more accurate
identification of quantum energy levels and the associated eigenstates.

###  4.4 Spatially local systems

A hallmark of most physical systems is _spatial locality_. Points close to one
another will usually have a stronger coupling than points that are far from
one another. For example, basic physical processes such as convection and
diffusion are spatially local. In most DMD applications, the data are
collected from a spatial grid, but the grid rarely plays a role in the DMD
analysis beyond forming the energy norm [68]. In fact, the output of the exact
DMD algorithm remains invariant to unitary transformations on the rows and
columns of the data matrices, so that randomly shuffling the rows (i.e., the
spatial locations) will result in identical DMD models with the corresponding
shuffling of the rows of the modes [57, 78]. Knowledge of the underlying grid
from which the data were collected enables us to bias the learning process to
prioritise relationships that we expect to have a strong local coupling.

Herein we consider the one-dimensional case; the analysis generalises
straightforwardly to higher dimensions. We assume that the entries of
$\boldsymbol{x}$ are samples of a function $u$ taken at $n$ grid points
$\\{\xi_{i}\\}$. The grid points can be arranged sequentially so that:
$\xi_{i}<\xi_{i+1}$. By the spatial locality principle, we expect states that
are spatially close to one another to have a stronger coupling than states
that are spatially far from one another. For example, we may expect entries
close to the diagonal of $\boldsymbol{A}$ to be larger than entries far from
the diagonal. Then the entries of $\boldsymbol{A}$ would satisfy

$\displaystyle|A_{i,j}|\geq|A_{i,k}|\quad\quad\textnormal{for
}\quad|j-i|<|k-i|.$ (26)

Equation (26) is merely a heuristic, and we would expect it to hold on average
rather than for each $(i,j,k)$. Additionally, (26) is a difficult and
expensive condition to implement in practice as it involves
$\mathcal{O}(n^{2})$ inequality constraints. Instead, we consider an
alternative version of spatial locality where we only permit coupling between
states that are sufficiently close:

$\displaystyle|A_{i,k}|=0\quad\quad\textnormal{for }\quad d<|k-i|.$ (27)

Equation (27) describes a $d$-diagonal matrix. If $d=0$ then states can only
affect themselves and $\boldsymbol{A}$ is diagonal (comparable to (17)). A
more interesting case is if $d=1$ and we only allow coupling between adjacent
states. Then $\boldsymbol{A}$ is tridiagonal: the entries of $\boldsymbol{A}$
are zeros except the leading, upper and lower diagonals:

$\displaystyle\boldsymbol{A}=\begin{bmatrix}\beta_{1}&\gamma_{1}\\\
\alpha_{2}&\beta_{2}&\gamma_{2}\\\ &\alpha_{3}&\ddots&\ddots\\\
&&\ddots&\ddots&\gamma_{n-1}\\\ &&&\alpha_{n}&\beta_{n}\end{bmatrix}.$ (28)

We now solve the optimization problem (3) when $\mathcal{M}$ is the manifold
of tridiagonal matrices. Appendix C includes many solutions for more general
diagonal-like structures including longer-range and variable coupling (§ C.1),
spatially periodic local systems (§ C.2), self-adjoint (§ C.3), total-least
squares (§ C.4), and weaker local structures (§ C.5). Since the rows of
$\boldsymbol{A}$ are decoupled, we can row-wise expand the Frobenius norm in
(3) to obtain $n$ smaller minimization problems:

$\displaystyle\operatorname*{arg\rm{}min}_{\alpha_{i},\,\beta_{i},\,\gamma_{i}}\|\alpha_{i}\tilde{\boldsymbol{x}}_{i-1}+\beta_{i}\tilde{\boldsymbol{x}}_{i}+\gamma_{i}\tilde{\boldsymbol{x}}_{i+1}-\tilde{\boldsymbol{y}}_{i}\|_{2}\qquad\qquad\textnormal{for
}1\leq i\leq n,$ (29)

with the convention $\alpha_{1}=\gamma_{n}=0$ and
$\tilde{\boldsymbol{x}}_{0}=\tilde{\boldsymbol{x}}_{n+1}=\boldsymbol{0}$.
Recall that $\tilde{\boldsymbol{x}}_{i}$ and $\tilde{\boldsymbol{y}}_{i}$ are
the $i$-th rows of $\boldsymbol{X}$ and $\boldsymbol{Y}$ respectively (4).
Each minimization problem (29) has the minimum-norm solution

$\displaystyle\begin{bmatrix}\alpha_{i}&\beta_{i}&\gamma_{i}\end{bmatrix}=\tilde{\boldsymbol{y}}_{i}\begin{bmatrix}\tilde{\boldsymbol{x}}_{i-1}\\\
\tilde{\boldsymbol{x}}_{i}\\\
\tilde{\boldsymbol{x}}_{i+1}\end{bmatrix}^{\dagger}\qquad\qquad\textnormal{for
}1\leq i\leq n.$ (30)

Each of the $n$ minimizations costs $\mathcal{O}(m)$ so the total cost is
$\mathcal{O}(mn)$.

####  4.4.1 Example: data-driven resolvent analysis

We consider a system of the form

$\displaystyle u_{t}=\mathcal{A}u+f(\xi,t)$ (31)

and aim to design a control strategy for the forcing $f$ via input-output
analysis [79, 68, 80, 81]. Fourier transforming in time and rearranging
produces $\hat{u}=R_{\omega}\hat{f}$ where
$R_{\omega}=\left(i\omega-\mathcal{A}\right)^{-1}$ is the resolvent of
$\mathcal{A}$ and the hats $\,\hat{\cdot}\,$ denote the Fourier transform in
time. We wish to understand the form of the forcing $\hat{f}$ that produces
the largest response in $\hat{u}$. Formally, the optimal forcing ($\phi_{1}$),
response ($\psi_{1}$), and gain ($\sigma_{1}$) for each frequency component
satisfy

$\displaystyle\phi_{1}=\operatorname*{arg\rm{}max}_{\|\phi_{1}\|_{2}=1}\|R_{\omega}\phi_{1}\|_{2},\qquad\sigma_{1}=\|R_{\omega}\phi_{1}\|_{2},\qquad\psi_{1}=\sigma_{1}^{-1}R_{\omega}\phi_{1}.$
(32)

We may also investigate the higher-order singular triplets
$(\sigma_{j},\,\phi_{j},\,\psi_{j})$ that satisfy equivalent conditions with
the additional requirement of orthogonality of $\\{\phi_{j}\\}$ and
$\\{\psi_{j}\\}$. Equivalently, we wish to determine the Hilbert–Schmidt
decomposition of $R_{\omega}$ [82]:

$\displaystyle
R_{\omega}=\sum_{j=1}^{\infty}\sigma_{j}\psi_{j}(\xi)\left<\phi_{j},\,\cdot\right>,$
(33)

which is analogous to the SVD of the discretization of $R_{\omega}$.

Recently, the authors have proposed ‘data-driven resolvent analysis’ as an
efficient method for learning the singular triplets of a system purely from
data measurements [68]. The method is analogous to DMD, except instead of
computing the eigenmodes via eigendecomposition, we compute the resolvent
modes via a singular value decomposition. In the following example, we
demonstrate that piDMD can be integrated into data-driven resolvent analysis
to provide physics-informed diagnostic information about the optimal forcings
and responses of the system.

We consider the problem of designing a strategy to control the concentration
of a solute that is governed by unknown dynamics. The solute concentration $u$
is governed by the convection-diffusion operator
$\mathcal{A}u=u_{\xi\xi}+a(\xi)u_{\xi}$ where $a$ represents a non-constant
convection term (here taken to be a random function) and boundary conditions
$u_{\xi}(-1)=u(1)=0$. A typical simulation is illustrated in the left panel of
figure 7. The convection-diffusion operator is a canonical example of a non-
normal system where the typical eigenvalue analysis fails to provide
meaningful information [83, 84] . Our goal is to determine the forcings
$\hat{f}$ that produce the largest increase in solute in the frequency domain.

Figure 7: Data-driven resolvent analysis of the convection-diffusion equation
for DMD and piDMD. The piDMD model is constrained to be tri-diagonal and
therefore respects the spatial locality of the system. Incorporating this
physical structure results a significant improvement in the identification of
the resolvent norm, response modes, and forcing modes.

We solve the convection-diffusion equation in Chebfun [85] and collect 2,000
snapshots of $u$ evaluated at 100 evenly spaced grid points. The snapshots
matrices $\boldsymbol{X}$ and $\boldsymbol{Y}$ are formed from measurements of
$u$ and $u_{t}$ respectively. We apply data-driven resolvent analysis with DMD
and piDMD and compute the learned operators, resolvent norm, and forcing and
response modes. For both algorithms, we approximate $u_{t}$ from the
measurements of $u$. The piDMD model is constrained to be tridiagonal, and the
DMD model has rank 40. From the visualizations of the learned models in the
center of figure 7, we see that the standard DMD identifies spurious
relationships between the states, whereas piDMD, by construction, identifies a
stronger coupling between adjacent states. The model learned by piDMD is
reminiscent of the finite difference matrix; indeed, local piDMD may be viewed
as a data-driven discovery of a finite difference stencil.

The benefit of piDMD is evidenced in its accurate identification of the
resolvent modes (for $\omega=1$), as seen in the right panel of figure 7. The
true resolvent modes are computed in Chebfun (see section 12 of [86]). The
modes learned by standard DMD are highly inaccurate for all $j$, and are
completely incoherent for $j>2$. In contrast, piDMD accurately determines the
mode shapes, and associated gains for a range of frequencies. The physics-
informed approach even determines the correct boundary conditions of the
response and forcing modes, and resolves the complicated structures associated
with the non-normality of the system. In particular, enforcing spatial
locality enables piDMD to identify modes that are poorly represented by the
data, which allows for diagnostic analysis that generalises outside of the
training regime. Equipped with these data-driven modes, one could now design a
low-rank control strategy for the solute concentration.

Analogously to finite difference methods, the accuracy of local piDMD involves
a trade-off between the grid sizes in space and time. Future work should focus
on clarifying these issues.

###  4.5 Causal systems

Causality is the process by which the behaviour of one state (the cause)
influences another state (the effect). Many systems posses a spatially causal
structure where each state depends only on ‘upstream’ states and are
unaffected by ‘downstream’ states. Specifically, we may express the $j$th
state as a function of the $n-j$ upstream states only:
$y_{j}=f(x_{j},\,x_{j+1},\,\dots\,x_{n})$. For example, $y_{n}$ depends only
on $x_{n}$ whereas $y_{1}$ depends only all the elements of $\boldsymbol{x}$.
Accordingly, if we know that the system at hand possesses a spatially causal
structure we may seek a piDMD model that is upper triangular. Such structures
are typical of a system with causal features, for example, time-delay
coordinates [87], feed-forward/strict-feedback systems [88, 89], and other
control systems [90].

As in section 4.4, we may decouple the rows of $\boldsymbol{A}$ in (3) and
expand the cost function row-wise to obtain $n$ smaller minimization problems:

$\displaystyle\operatorname*{arg\rm{}min}_{\boldsymbol{A}_{j,j:n}}\|\boldsymbol{A}_{j,j:n}\boldsymbol{X}_{j:n,\,:}-\tilde{\boldsymbol{y}}_{j}\|_{2}\qquad\qquad\textnormal{for
}1\leq j\leq n.$ (34)

where we have employed the Matlab colon notation ${j\,{:}\,n}$ to indicate
entries from $j$ to $n$ and $:$ to indicate entries from $1$ to $n$. Here,
$\tilde{\boldsymbol{y}}_{j}$ represents the $j$-th row of $\boldsymbol{Y}$.
The (minimum norm) solution of (34) is obtained via the pseudoinverse of each
sub-block of $\boldsymbol{X}$ as

$\displaystyle\boldsymbol{A}_{j,j:n}=\tilde{\boldsymbol{y}}_{j}\boldsymbol{X}_{j:n,\,:}^{\dagger}\qquad\qquad\textnormal{for
}1\leq j\leq n.$ (35)

Naively evaluating (35) for each $j$ would require building $n$ pseudoinverses
for a total cost of $\mathcal{O}(n^{2}m\min(m,n))$. Fortunately, consecutive
blocks of $\boldsymbol{X}$ (e.g. $\boldsymbol{X}_{j-1:n,\,:}$ and
$\boldsymbol{X}_{j:n,\,:}$) are related by a rank-1 update, so we can solve
every problem in (35) for a total of $\mathcal{O}(nm\min(n,m))$ operations.
(see appendix D.1 for further details). However, this approach is very
numerically unstable when applied to realistic data with high condition
number. Thus, in appendix D.2 we derive a more stable alternative.

We illustrate this causal structure on a simple example defined by the
Volterra-type integro-differential equation

$\displaystyle\frac{\partial u}{\partial
t}(\xi,t)=\int_{-1}^{\xi}K(\xi,\nu)u(\nu,t),\textrm{d}\nu\qquad\qquad-1\leq\xi\leq
1.$ (36)

Similar systems have been used to model transmission lines in neural networks
during bursting activity [91] and the spread of disease in epidemics [92].
Taking $K(\xi,\nu)=\sqrt{1-\xi^{2}}\sqrt{1-\nu^{2}}$ and
$u(\xi,0)=\exp(-\xi^{2})$, we numerically simulate (36) and construct
$\boldsymbol{X}$ and $\boldsymbol{Y}$ from measurements of $u$ only (velocity
measurements are not used). The results of piDMD are illustrated in the
penultimate row of figure 3, and indicate that piDMD is able to learn the
leading eigenvalues while standard DMD fails. Note that the true spectrum here
is continuous.

##  5 Conclusions and outlook

This work presents physics-informed dynamic mode decomposition (piDMD), a
physics-aware modal decomposition technique that extracts coherent structures
from high-dimensional time-series data. Rephrasing the DMD regression (2) as a
Procrustes problem (3) shows that partially-known physics can be incorporated
into the DMD framework by enforcing a matrix manifold constraint determined by
known physics. We have applied piDMD to five of the most fundamental physical
principles: conservation laws, self-adjointness, shift-invariance, locality,
and causality. In several of these cases, we have presented new ‘exact’
solutions of the corresponding optimization problems. The examples presented
demonstrate that piDMD can exhibit superior performance compared to classical
DMD methods that do not account for the physics of the system.

The framework developed herein offers exciting new opportunities and
challenges for data-driven dynamical systems, numerical linear algebra, convex
optimization, and statistics. Fortunately, many research advances developed
for classical DMD can also be leveraged to enhance piDMD. We conclude the
paper with a brief discussion of the limitations of piDMD, and suggest
promising future research directions.

###  5.1 Limitations and extensions

Each of the fours steps of the piDMD framework outlined in section 3 pose
distinct challenges. In some scenarios where the physics is poorly understood,
determining suitable physical laws to impose on the model (step 1, ‘modeling’)
can be challenging. For problems with intricate geometries and multiple
dimensions, interpreting the physical principle as a matrix manifold (step 2,
‘interpretation’) can be the roadblock, as the manifold can become exceedingly
complicated. We have already commented on the challenges and opportunities
posed by step 3 (‘optimization’), and suggested algorithmic solutions [51].
Finally, the typically large state dimension and number of samples can
obfuscate step 4 when the matrix manifold does not easily admit the desired
diagnostics. For example, if $\boldsymbol{A}$ has a complicated banded
structure, as in local piDMD, it may not be amenable to a fast
eigendecomposition or SVD.

The solutions presented herein exhibit various levels of sensitivity to noise.
As is typical of DMD methods, most solutions are unbiased with respect to
$\boldsymbol{Y}$, but can be quite sensitive to noise in $\boldsymbol{X}$.
Some of the solutions we have presented are relatively insensitive (§4.2) or
can be reformulated in the total least squares sense (§A.4, §C.4). A full
characterization of the sensitivity of piDMD should be performed on a case-by-
case basis for different manifold constraints, and will be the subject of
future investigations.

The Procrustes problem (3) is not the only means of incorporating physical
principles into the DMD process. Closely related to Procrustes problems are
‘nearest matrix’ problems [93] where, given a matrix $\hat{\boldsymbol{A}}$,
we seek the closest matrix on some matrix manifold $\mathcal{M}$:

$\displaystyle\operatorname*{arg\rm{}min}_{\boldsymbol{A}\in\mathcal{M}}\|\boldsymbol{A}-\hat{\boldsymbol{A}}\|.$
(37)

Equation (37) can be viewed as a special case of the Procrustes problem (3)
where $\boldsymbol{Y}=\hat{\boldsymbol{A}}$ and
$\boldsymbol{X}=\boldsymbol{I}$. As such, nearest matrix problems are
generally easier to solve than Procrustes problems. The nearest matrix problem
(37) can be used to project any DMD model onto the nearest physically
consistent model, or within optimization routines, such as proximal gradient
methods or constrained gradient descent, to constrain the solution to the
feasible set.

As presented, piDMD strictly requires that the model $\boldsymbol{A}$ lies on
the chosen matrix manifold $\mathcal{M}$. In many applications, such as when
the data are very noisy or the physical laws and/or constraints are only
approximately understood, it may be more appropriate to merely ‘encourage’
$\boldsymbol{A}$ toward $\mathcal{M}$. In such cases, the piDMD regression (3)
becomes

$\displaystyle\operatorname*{arg\rm{}min}_{\boldsymbol{A}}\|\boldsymbol{Y}-\boldsymbol{A}\boldsymbol{X}\|_{F}+\lambda\,R(\boldsymbol{A}),$
(38)

where the first term is the reconstruction loss, $R$ is a physics-informed
regularizer, and $\lambda$ is a user-defined constant that tunes the relative
importance of the reconstruction and regularization. For example,
$R(\boldsymbol{A})$ could represent the distance between $\boldsymbol{A}$ and
$\mathcal{M}$, which can be computed by solving the nearest matrix problem
(37). An example of such a problem is solved in appendix C.5.

DMD models may be trained in discrete or continuous time. In discrete time
with a constant sampling rate, the snapshots in $\boldsymbol{Y}$ are
$\boldsymbol{y}_{k}=\boldsymbol{x}_{k+1}$, whereas in the continuous time we
have $\boldsymbol{y}_{k}=\dot{\boldsymbol{x}}_{k}$. If the continuous time
operator is $\boldsymbol{A}_{C}$, then the discrete time operator is
$\boldsymbol{A}_{D}=\exp(\Delta t\boldsymbol{A}_{C})$. In special cases, the
matrices $\boldsymbol{A}_{C}$ and $\boldsymbol{A}_{D}$ may lie on the same
manifold, such as the upper-triangular and circulant manifolds; however, in
general, they do not. For example, if $\boldsymbol{A}_{C}$ is tridiagonal,
then $\boldsymbol{A}_{D}$ is not generally tridiagonal. Accordingly, if only
discrete-time measurements are available but the matrix manifold is imposed on
the continuous-time operator, then the optimization problem to solve is

$\displaystyle\operatorname*{arg\rm{}min}_{\boldsymbol{A}\in\mathcal{M}}\|\boldsymbol{Y}-\exp(\Delta
t\boldsymbol{A})\boldsymbol{X}\|_{F},$ (39)

which is a much more difficult optimization problem.

## Acknowledgements

P.J.B. acknowledges insightful conversations with Suvrit Sra, Tasuku Soma and
Andrew Horning. S.L.B. acknowledges valuable discussions with Jean-Christophe
Loiseau. The authors acknowledge support from the Army Research Office (ARO
W911NF-17-1-0306) and the National Science Foundation AI Institute in Dynamic
Systems (Grant No. 2112085).

## Appendices

The following appendices provide further details of the optimization problems
considered in this paper.

##  Appendix A Further details of shift-invariant systems

Here we detail several extensions to the shift-invariant DMD problem
introduced in section 4.1. Figure 8 compares these alternative solutions for
training data from the 2-D advection equation.

###  A.1 Symmetric, skew-symmetric and unitary circulant matrices

Now consider the case where $\boldsymbol{A}$ is symmetric and circulant. An
application of (15) reveals that a circulant matrix is symmetric if and only
if its eigenvalues $\\{\hat{\boldsymbol{a}}_{i}\\}$ are real. Thus, the
minimization problem is

$\displaystyle\operatorname*{arg\rm{}min}_{\Im[\hat{\boldsymbol{a}}]=\boldsymbol{0}}\|\textrm{diag}(\hat{\boldsymbol{a}}){\boldsymbol{\mathcal{X}}}-{\boldsymbol{\mathcal{Y}}}\|_{F}^{2}.$
(40)

Expanding the Frobenius norm row-wise and enforcing stationarity yields the
optimal value as

$\displaystyle\hat{a}_{j}=\frac{\Re\left[\tilde{\boldsymbol{\mathcal{Y}}}_{j}\tilde{\boldsymbol{\mathcal{X}}}_{j}^{\ast}\right]}{\left\|\tilde{\boldsymbol{\mathcal{X}}}_{j}\right\|^{2}_{2}}.$
(41)

The case of a skew-symmetric circulant matrix is similar except the
eigenvalues are now imaginary; the solution is

$\displaystyle\hat{a}_{j}=\textrm{i}\frac{\Im\left[\tilde{\boldsymbol{\mathcal{Y}}}_{j}\tilde{\boldsymbol{\mathcal{X}}}_{j}^{\ast}\right]}{\left\|\tilde{\boldsymbol{\mathcal{X}}}_{j}\right\|^{2}_{2}}.$
(42)

Finally, for a unitary circulant matrix (when the eigenvalues lie on the unit
circle), we have

$\displaystyle\hat{a}_{j}=\frac{\tilde{\boldsymbol{\mathcal{Y}}}_{j}\tilde{\boldsymbol{\mathcal{X}}}_{j}^{\ast}}{\left\|\tilde{\boldsymbol{\mathcal{Y}}}_{j}\tilde{\boldsymbol{\mathcal{X}}}_{j}^{\ast}\right\|_{2}}.$
(43)

Figure 8: Comparison of exact DMD to shift-invariant piDMD and some variants.
The model are trained on snapshots of a solution to the 2-D advection
equation. The snapshots are contaminated by 2% Gaussian noise. The low-rank
model has rank 45.

###  A.2 Best rank-$r$ circulant approximant

Now consider the problem where $\boldsymbol{A}$ is both circulant and rank
$r\leq n$. Since the rank of a matrix is the number of non-zero eigenvalues,
the rank of a circulant matrix is the number of non-zero entries in
$\hat{\boldsymbol{a}}$. Thus, only $r$ entries in $\hat{\boldsymbol{a}}$ are
permitted to be non-zero and the minimization problem (17) becomes

$\displaystyle\operatorname*{arg\rm{}min}_{\|\hat{\boldsymbol{a}}\|_{0}=r}\|\textrm{diag}(\hat{\boldsymbol{a}})\boldsymbol{\mathcal{X}}-\boldsymbol{\mathcal{Y}}\|_{F}.$
(44)

Equivalently, given the solution (15), we must select $r$ eigenvalues to
maintain and set all other eigenvalues to zero. To minimise the cost function,
the $r$ chosen eigenvalues should contribute to the largest reduction in the
cost function. The extent to which each eigenvector minimises the cost
function is characterised by their residual:

$\displaystyle{\rho}_{j}=\sqrt{\|\tilde{\boldsymbol{\mathcal{Y}}}_{j}\|_{2}^{2}-\|\tilde{\boldsymbol{\mathcal{X}}}_{j}-\hat{a}_{j}\tilde{\boldsymbol{\mathcal{X}}}_{j}\|_{2}^{2}}=\frac{\Re\left[\tilde{\boldsymbol{\mathcal{Y}}}_{j}\tilde{\boldsymbol{\mathcal{X}}}_{j}^{\ast}\right]}{\|\tilde{\boldsymbol{\mathcal{X}}}_{j}\|_{2}}.$
(45)

A large residual indicates a large reduction in the cost function and vice
versa. Therefore, the $r$ non-zero eigenvalues to maintain are the $r$ with
the largest residuals.

###  A.3 Learning shift-invariant operators with non-equispaced samples

Here we extend the method of section 4.1 to allow spatial samples that are not
equally spaced. We begin by formulating the Procrustes problem (3) in
continuous space. Suppose that we have access to measurements of functions
$u(\xi,t)$ and $v(\xi,t)$ for all $\xi\in[-1,1]$ but at discrete points
$t\in\\{t_{j}\ |\ j=1,\dots,m\\}$. Then, the continuous-space piDMD regression
(3) is

$\displaystyle\operatorname*{arg\rm{}min}_{\mathcal{A}\in\mathcal{M}}\sum_{j=1}^{m}\left\|v(\xi,t_{j})-\mathcal{A}u(\xi,t_{j})\right\|_{2}^{2},$
(46)

where the continuous $2$-norm is

$\displaystyle\|f(\xi)\|_{2}^{2}=\int_{1}^{1}|f(\xi)|^{2}\textrm{d}\xi.$ (47)

We now constrain $\mathcal{A}$ to be a shift-invariant operator on a periodic
domain. Such operators are diagonalized by exponential functions (§4.1) so we
may write

$\displaystyle\mathcal{A}f(\xi)=\sum_{k=-\infty}^{\infty}\hat{a}_{k}\textrm{e}^{\textrm{i}\pi
k\xi}\hat{f}(k),$ (48)

where $\hat{f}(k)$ is the finite Fourier transform of $f$ at wavenumber $k$:

$\displaystyle\hat{f}(k)=\int_{-1}^{1}f(\eta)\textrm{e}^{-\textrm{i}\pi
k\eta}\textrm{d}\eta.$ (49)

Furthermore, if $u$ and $v$ are sufficiently smooth in $\xi$ (e.g. they are
Hölder continuous [94]) then we may write

$\displaystyle
u(\xi,t)=\sum_{k=-\infty}^{\infty}\hat{u}(k,t)\textrm{e}^{\textrm{i}\pi
k\xi},\qquad
v(\xi,t)=\sum_{k=-\infty}^{\infty}\hat{v}(k,t)\textrm{e}^{\textrm{i}\pi
k\xi}.$ (50)

Applying Parseval’s theorem in concert with (48) and (50) shows that the
minimization problem (46) is equivalent to

$\displaystyle\operatorname*{arg\rm{}min}_{\hat{a}_{k}}\sum_{j=1}^{m}\left\|\sum_{k=-\infty}^{\infty}\left(\hat{v}(k,t_{j})-\hat{a}_{k}\hat{u}(k,t_{j})\right)\textrm{e}^{\textrm{i}\pi
k\xi}\right\|_{2}^{2}$
$\displaystyle=\operatorname*{arg\rm{}min}_{\hat{a}_{k}}\sum_{k=-\infty}^{\infty}\sum_{j=1}^{m}\left|\hat{v}(k,t_{j})-\hat{a}_{k}\hat{u}(k,t_{j})\right|^{2}.$
(51)

Furthermore, writing
$\hat{\boldsymbol{u}}=\begin{bmatrix}\hat{u}(k,t_{1})&\cdots&\hat{u}(k,t_{m})\end{bmatrix}$
and
$\hat{\boldsymbol{v}}=\begin{bmatrix}\hat{v}(k,t_{1})&\cdots&\hat{v}(k,t_{m})\end{bmatrix}$
transforms the right side of (51) to

$\displaystyle\operatorname*{arg\rm{}min}_{\hat{a}_{k}}\sum_{k=-\infty}^{\infty}\left\|\hat{\boldsymbol{v}}(k)-\hat{a}_{k}\hat{\boldsymbol{u}}(k)\right\|_{2}^{2},$

whose solution is

$\displaystyle\hat{a}_{k}=\hat{\boldsymbol{v}}(k)\hat{\boldsymbol{u}}^{\dagger}(k)=\left.\hat{\boldsymbol{v}}(k)\hat{\boldsymbol{u}}^{\ast}(k)\middle/\|\hat{\boldsymbol{u}}(k)\|_{2}^{2}\right..$
(52)

Thus, (52) solves the continuous-space shift-invariant piDMD regression (46).
Note that the solution (52) depends only on the Fourier coefficients of $u$
and $v$ (50).

Now, suppose that we only have measurements of $u$ and $v$ at arbitrary
discrete grid points $\xi\in\\{\xi_{l}\\}$. We can then approximate the
solution (52) by approximating the coefficients of the Fourier coefficients
$\hat{\boldsymbol{v}}(k)$ and $\hat{\boldsymbol{u}}(k)$. The samples at
$\\{\xi_{l}\\}$ correspond to the quadrature rule for the Fourier coefficients
(49)

$\displaystyle\hat{f}(k)\approx\sum_{l=1}^{n}w_{l}\
f(\xi_{l})\textrm{e}^{-\textrm{i}\pi k\xi_{l}}$ (53)

where $w_{k}=(\xi_{k+1}-\xi_{k-1})/2\ (\mathrm{mod}\ 1)$ are the quadrature
nodes for the grid. We can now use (53) to approximate $\hat{\boldsymbol{v}}$
and $\hat{\boldsymbol{v}}$ and thereby approximate the solution (52). The sums
(53) can be computed efficiently with the non-uniform discrete Fourier
transform (NUDFT–II, [70, 95]).

###  A.4 Total least-squares circulant problem

A equivalent statement of the Procrustes problem (3) is

$\displaystyle\operatorname*{arg\rm{}min}_{\boldsymbol{A}\in\mathcal{M}}\|\boldsymbol{R}\|_{F}\quad\textnormal{subject
to }\boldsymbol{A}\boldsymbol{X}=\boldsymbol{Y}+\boldsymbol{R}.$ (54)

As such, the formulation of (3) assumes assumes that any errors are in
$\boldsymbol{Y}$ only. Alternatively, we can consider errors in both
$\boldsymbol{X}$ and $\boldsymbol{Y}$ and reformulate (3) as

$\displaystyle\operatorname*{arg\rm{}min}_{\boldsymbol{A}\in\mathcal{M}}\left\|\left[\boldsymbol{R}\;\;\;\boldsymbol{S}\right]\right\|_{F}\quad\textnormal{subject
to
}\boldsymbol{A}\left(\boldsymbol{X}+\boldsymbol{S}\right)=\boldsymbol{Y}+\boldsymbol{R}$
(55)

This approach is usually more principled as there will be inevitable errors in
measurements of $\boldsymbol{X}$ and $\boldsymbol{Y}$. Here we show that we
can adapt the circulant solution of §4.1 to account for noise in both sets of
measurements. Again, by the invariance of the Frobenius norm with respect to
unitary transformations, (55) is equivalent to

$\displaystyle\operatorname*{arg\rm{}min}_{\hat{\boldsymbol{a}}}\left\|\left[\boldsymbol{\mathcal{R}}\;\;\;\boldsymbol{\mathcal{S}}\right]\right\|_{F}\quad\textnormal{subject
to
}\qquad\textrm{diag}\left(\hat{\boldsymbol{a}}\right)\left(\boldsymbol{\mathcal{X}}+\boldsymbol{\mathcal{S}}\right)=\boldsymbol{\mathcal{Y}}+\boldsymbol{\mathcal{R}}$
(56)

where $\boldsymbol{\mathcal{R}}=\boldsymbol{\mathcal{F}}^{\ast}\boldsymbol{R}$
and $\boldsymbol{\mathcal{S}}=\boldsymbol{\mathcal{F}}^{\ast}\boldsymbol{S}$.
Again, the rows decouple to produce $n$ smaller problems:

$\displaystyle\operatorname*{arg\rm{}min}_{\hat{a}_{i}}\left\|\left[\tilde{\boldsymbol{\mathcal{R}}}_{i}\;\;\;\tilde{\boldsymbol{\mathcal{S}}}_{i}\right]\right\|_{F}\quad\textnormal{subject
to
}\quad\hat{a}_{i}\left(\tilde{\boldsymbol{\mathcal{X}}}_{i}+\tilde{\boldsymbol{\mathcal{S}}}_{i}\right)=\tilde{\boldsymbol{\mathcal{Y}}}_{i}+\tilde{\boldsymbol{\mathcal{R}}}_{i}.$
(57)

for $1\leq i\leq n$. These problems may now be solved in the total least
squares sense [35, 96].

We will briefly characterize the scenario where the total least-squares
circulant model is statistically optimal. By optimal here, we mean that the
model converges to the true model in probability as the number of samples
increases [97]. In the unconstrained case, it is known that this is the case
when the entries of $\boldsymbol{R}$ and $\boldsymbol{S}$ are random variables
with zero mean and the covariance matrix of
$\textnormal{vec}\left([\boldsymbol{R}\;\;\;\boldsymbol{S}]\right)$ is a
multiple of the identity matrix [97]. Therefore, the solution of the total
least squares problem (57) will be statistically optimal if

$\displaystyle\textnormal{Var}\left({\boldsymbol{\mathcal{R}}}_{i,j}\right)=\textnormal{Var}\left({\boldsymbol{\mathcal{S}}}_{i,j}\right)\qquad\textnormal{and}\qquad\textnormal{Cov}\left({\boldsymbol{\mathcal{R}}}_{i,j},{\boldsymbol{\mathcal{S}}}_{i,j}\right)=0.$
(58)

In words the circulant-TLS problem is statistically optimal if errors in
$\boldsymbol{R}$ and $\boldsymbol{S}$ have equal variance for a given
wavenumber and the errors between
$\textnormal{Cov}\left({\boldsymbol{\mathcal{R}}}_{i,j},{\boldsymbol{\mathcal{S}}}_{i,j}\right)$
are independent. These conditions are met if $\boldsymbol{R}$ and
$\boldsymbol{S}$ are independent random variables of equal variance $\sigma$,
in which case we have
$\textnormal{Var}\left({\boldsymbol{\mathcal{R}}}_{i,j}\right)=\textnormal{Var}\left({\boldsymbol{\mathcal{S}}}_{i,j}\right)=\sigma^{2}$.
Therefore, if the unconstrained total least squares model is statistically
optimal then so is the circulant piDMD model.

###  A.5 Toeplitz and Hankel matrices

Here we adapt the optimization of section §4.1 to Toeplitz matrices, which
take the form

$\displaystyle\boldsymbol{A}=\begin{bmatrix}a_{0}&a_{-1}&\cdots&a_{-(n-1)}\\\
a_{1}&a_{0}&\ddots&\vdots\\\ \vdots&\ddots&\ddots&a_{-1}\\\
a_{n-1}&\cdots&a_{1}&a_{0}\end{bmatrix}.$ (59)

Similarly to the circulant case, $\boldsymbol{A}_{i,j}=a_{i-j}$ but now the
subscripts not are considered modulo $n$. The Procrustes problem for a
Toeplitz matrix was considered by [98]. We produce an alternative solution
based on the fast Fourier transform, analogous to the solution for the
circulant Procrustes problem (§4.1). A Toeplitz matrix can be embedded in a
circulant matrix as

$\displaystyle{\boldsymbol{C}=\left[\begin{array}[]{c|c}\boldsymbol{A}&\tilde{\boldsymbol{A}}\\\
\hline\cr\tilde{\boldsymbol{A}}&\boldsymbol{A}\end{array}\right]}\qquad\textrm{
where }\qquad\tilde{\boldsymbol{A}}=\begin{bmatrix}{0}&a_{n-1}&\cdots&a_{1}\\\
a_{-(n-1)}&0&\ddots&\vdots\\\ \vdots&\ddots&\ddots&a_{n-1}\\\
a_{-1}&\cdots&a_{-(n-1)}&0\end{bmatrix}$ (62)

Similarly to §4.1, the circulant matrix $\boldsymbol{C}$ is diagonalized by
the DFT matrix $\boldsymbol{\mathcal{F}}$ (which is now $2n\times 2n$). Thus,
we may write
$\boldsymbol{C}=\boldsymbol{\mathcal{F}}^{\ast}\textrm{diag}(\hat{\boldsymbol{a}})\boldsymbol{\mathcal{F}}$
and the minimization (3) becomes

$\displaystyle\operatorname*{arg\rm{}min}_{\hat{\boldsymbol{a}}}\left\|\hat{\boldsymbol{I}}^{\ast}\textrm{diag}(\hat{\boldsymbol{a}}){\boldsymbol{\mathcal{X}}}-\boldsymbol{Y}\right\|_{F}^{2}$
(63)

where
$\hat{\boldsymbol{I}}=\boldsymbol{\mathcal{F}}\begin{bmatrix}\boldsymbol{I}\\\
\boldsymbol{0}\end{bmatrix}$ and
${\boldsymbol{\mathcal{X}}}=\boldsymbol{\mathcal{F}}\begin{bmatrix}\boldsymbol{X}\\\
\boldsymbol{0}\end{bmatrix}$. It can be shown that the optimal coefficients
$\hat{\boldsymbol{a}}$ solve the linear system

$\displaystyle\boldsymbol{H}\hat{\boldsymbol{a}}=\boldsymbol{d}\qquad\textrm{
where
}\qquad\boldsymbol{H}=\left(\hat{\boldsymbol{I}}\hat{\boldsymbol{I}}^{\ast}\odot\overline{{\boldsymbol{\mathcal{X}}}{\boldsymbol{\mathcal{X}}}^{\ast}}\right),\quad\boldsymbol{d}=\textrm{diag}\left(\hat{\boldsymbol{I}}\boldsymbol{Y}{\boldsymbol{\mathcal{X}}}^{\ast}\right)$
(64)

and $\odot$ represents element-wise multiplication (i.e. the Hadamard
product). Note that a similar optimization arises in the computation of the
amplitudes of DMD modes [99]. While mathematically correct, numerical
implementations of (64) can be difficult. Specifically, issues arise because
the Gramian matrix
${\boldsymbol{\mathcal{X}}}{\boldsymbol{\mathcal{X}}}^{\ast}$ will typically
be large and ill-conditioned. Thus, we expect that one could find more
efficient algorithms to calculate $\boldsymbol{A}$, possibly by adapting
methods for fast solutions of Toeplitz systems [100, 101, 102].

This approach also applies to Hankel matrices, which take the form

$\displaystyle\boldsymbol{A}=\begin{bmatrix}a_{0}&a_{1}&\ldots&a_{n-1}\\\
a_{1}&a_{2}&\iddots&\vdots\\\ \vdots&\iddots&\iddots&a_{2n-3}\\\
a_{n-1}&\ldots&a_{2n-3}&a_{2n-2}\end{bmatrix}$ (65)

i.e. each ascending skew-diagonal is constant:
$\boldsymbol{A}_{i,j}=\boldsymbol{a}_{i+j-2}.$ Hankel matrices are upside-down
Toeplitz matrices: if $\boldsymbol{A}$ is Hankel, then
$\boldsymbol{J}\boldsymbol{A}$ is Toeplitz, where $\boldsymbol{J}$ is ths
upside-down identity matrix. Thus, the Procrustes problem for a Hankel matrix
is equivalent to the for a Toeplitz matrix with a suitable flipping of the
data.

##  Appendix B Further details of self-adjoint systems

###  B.1 Skew-symmetric systems

If $\boldsymbol{A}$ is constrained to be skew-symmetric then (22) and still
holds and
$\boldsymbol{C}=\boldsymbol{U}_{X}^{\ast}\boldsymbol{Y}\boldsymbol{V}_{X}$ but
now $\boldsymbol{L}$ is now defined by

$\displaystyle\boldsymbol{L}_{i,j}=-\overline{\boldsymbol{L}}_{j,i}=\begin{cases}\dfrac{-\sigma_{i}\overline{\boldsymbol{C}}_{j,i}+\sigma_{j}\boldsymbol{C}_{i,j}}{\sigma_{i}^{2}+\sigma_{j}^{2}}&\textnormal{if
}\sigma_{i}^{2}+\sigma_{j}^{2}\neq 0,\\\ \qquad\qquad
0&\textnormal{otherwise}.\end{cases}$ (66)

###  B.2 Sensitivity of symmetric piDMD

In light of the results of section 4.3, we now compare the sensitivities of
the symmetric piDMD model (22) and exact DMD model (7). In particular, we will
show that the elements of the piDMD model have lower variance than those of
the corresponding exact DMD model. Higham [41] also investigated the
sensitivity of the symmetric Procrustes problem, but from the perspective of
numerical stability; here we take a statistical perspective. In this section
only, we denote the exact DMD solution (7) as $\boldsymbol{A}^{e}$ and the
symmetric piDMD solution (22) as $\boldsymbol{A}^{s}$.

We assume that measurements of $\boldsymbol{Y}$ are contaminated with I. I. D.
Gaussian noise. For now, we neglect noise in $\boldsymbol{X}$, which could be
incorporated via a total least square approach [35, 34], By linearity, it is
sufficient to consider $\boldsymbol{Y}_{i,j}\sim\mathcal{N}(0,1)$. Again,
since exact DMD and symmetric piDMD are linear models, they are unbiased
estimators:

$\displaystyle\mathbb{E}\left(\boldsymbol{A}^{s}_{i,j}\right)=\mathbb{E}\left(\boldsymbol{A}^{e}_{i,j}\right)=0.$

We now compare the variance of the models. Extensive algebra produces explicit
expressions for the variances of $\boldsymbol{A}^{e}_{i,j}$ and
$\boldsymbol{A}^{s}_{i,j}$ as

$\displaystyle\mathrm{Var}\left({\boldsymbol{A}^{e}}_{i,j}\right)=\sum_{l=1}^{r}\frac{\boldsymbol{U}_{j,l}^{2}}{\sigma_{l}^{2}},$
(67)
$\displaystyle\mathrm{Var}(\boldsymbol{A}^{s}_{i,j})=\frac{1}{2}\sum^{r}_{k=1}\sum^{r}_{l=1}\frac{\left(\boldsymbol{U}_{i,k}\boldsymbol{U}_{j,l}+\boldsymbol{U}_{i,l}\boldsymbol{U}_{j,k}\right)^{2}}{\sigma_{k}^{2}+\sigma_{l}^{2}},$
(68)

where $\boldsymbol{U}$ is the matrix of left singular vectors of
$\boldsymbol{X}$, $\\{\sigma_{l}\\}$ are the associated singular values, and
$\textrm{rank}(\boldsymbol{X})=r$. Figure 9 compares these theoretical values
to empirical values for random data. These results are in agreement with
figure 6 and indicate that the variance of the piDMD model is smaller than
that of the exact DMD model.

Figure 9: Distribution of errors in the entries of the learned models
($A^{e,s}_{2,2}$, $A^{e,s}_{13,3}$ and $A^{e,s}_{17,6}$) for exact DMD and
self-adjoint piDMD. The block colors represent the frequency and the dashed
lines are the theoretical curves (67, 68).

We formalize this result in the following lemma and proof.

######  Lemma B.1.

If the measurements of $\boldsymbol{Y}$ are contaminated with I. I. D.
Gaussian noise, the variance of each element of the symmetric piDMD model (22)
is bounded above by the variance of the corresponding element of the exact DMD
model (7):

$\displaystyle\mathrm{Var}(\boldsymbol{A}_{i,j}^{s})\leq\mathrm{Var}({\boldsymbol{A}}_{i,j}^{e}).$
(69)

###### Proof.

Sedrakyan’s inequality [103] produces

$\displaystyle\frac{\left(\boldsymbol{U}_{i,k}\boldsymbol{U}_{j,l}+\boldsymbol{U}_{i,l}\boldsymbol{U}_{j,k}\right)^{2}}{\sigma_{k}^{2}+\sigma_{l}^{2}}$
$\displaystyle\leq\frac{\boldsymbol{U}_{i,k}^{2}\boldsymbol{U}_{j,l}^{2}}{\sigma_{l}^{2}}+\frac{\boldsymbol{U}_{i,l}^{2}\boldsymbol{U}_{j,k}^{2}}{\sigma_{k}^{2}}.$

Therefore,

$\displaystyle\mathrm{Var}(\boldsymbol{A}^{s}_{i,j})\leq\frac{1}{2}\sum^{r}_{k=1}\sum^{r}_{l=1}\left(\frac{\boldsymbol{U}_{i,k}^{2}\boldsymbol{U}_{j,l}^{2}}{\sigma_{l}^{2}}+\frac{\boldsymbol{U}_{i,l}^{2}\boldsymbol{U}_{j,k}^{2}}{\sigma_{k}^{2}}\right)=\sum^{r}_{k=1}\sum^{r}_{l=1}\frac{\boldsymbol{U}_{i,k}^{2}\boldsymbol{U}_{j,l}^{2}}{\sigma_{l}^{2}}=\left(\sum^{r}_{k=1}\boldsymbol{U}_{i,k}^{2}\right)\mathrm{Var}(\boldsymbol{A}^{e}_{i,j}).$
(70)

The result now follows since the rows of $\boldsymbol{U}$ have unit norm so
$\sum_{k=1}^{r}\boldsymbol{U}_{i,k}^{2}\leq 1$. ∎

This result indicates that the symmetric piDMD model is less sensitive to
noise than the exact DMD model when the noise is only in $\boldsymbol{Y}$.
Future work should investigate the case when both $\boldsymbol{X}$ and
$\boldsymbol{Y}$ are contaminated with noise, and clarify the effect of noise
on the learned spectrum.

##  Appendix C Further details of local systems

###  C.1 More general local structures

The local solution of section 4.4 naturally extends to more general local
structures. For example, it may be necessary to consider longer range
interactions between states and allow, say, a state to depend on other states
three grid points away. Additionally, there may be a variable local dependence
where different states can depend on different numbers of elements. For
example, a state near boundary (e.g. $x_{1}$ or $x_{n}$) may need to depend on
more than one other state in order to properly resolve the boundary behaviour.
To make this precise, suppose that state $i$ depends on $l(i)$ states to the
right and $u(i)$ states to the left above the leading diagonal (assuming that
the data is organised from left to right). Then, the solution for the optimal
matrix entries is

$\displaystyle\begin{bmatrix}A_{i-l(i),i}&\cdots&A_{i,i}&\cdots&A_{i+u(i),i}\end{bmatrix}=\tilde{\boldsymbol{y}}_{i}\begin{bmatrix}\tilde{\boldsymbol{x}}_{i-l(i)}\\\
\vdots\\\ \tilde{\boldsymbol{x}}_{i}\\\ \vdots\\\
\tilde{\boldsymbol{x}}_{i+u(i)}\end{bmatrix}^{\dagger}\qquad\qquad\textnormal{for
}1\leq i\leq n$ (71)

with the convention $A_{i,j}=0$ if the entry does not exist and
$\tilde{\boldsymbol{x}}_{i}=\boldsymbol{0}$ for $i<1$ or $i>n$. We recover the
tridiagonal solution (30) when $l(i)=u(i)=1$.

###  C.2 Local and periodic

The tridiagonal solution in section 4.4 effectively assumed that the system
had two endpoints that were disconnected from one another. However, the
spatial domain could be periodic, in which case the first state can affect the
$n$th state and vice versa. If we allow this structure then the matrix
$\boldsymbol{A}$ takes the form

$\displaystyle\boldsymbol{A}=\begin{bmatrix}\beta_{1}&\gamma_{1}&&&\alpha_{1}\\\
\alpha_{2}&\beta_{2}&\gamma_{2}\\\ &\alpha_{3}&\ddots&\ddots\\\
&&\ddots&\ddots&\gamma_{n-1}\\\
\gamma_{n}&&&\alpha_{n}&\beta_{n}\end{bmatrix}.$ (72)

The solution for the optimal constants is identical to those stated in (30)
except now the subscripts are evaluated modulo $n$.

###  C.3 Symmetric tridiagonal

Many spatially local systems are also self-adjoint. In such cases, we may
enforce that $\boldsymbol{A}$ is both tri-diagonal and symmetric:

$\displaystyle\boldsymbol{A}=\begin{bmatrix}\alpha_{1}&\beta_{1}\\\
\beta_{1}&\alpha_{2}&\beta_{2}\\\ &\beta_{2}&\ddots&\ddots\\\
&&\ddots&\ddots&b_{n-1}\\\ &&&b_{n-1}&\alpha_{n}\end{bmatrix}$ (73)

. The cost function (3) is then expressible as

$\displaystyle\|\boldsymbol{Y}-\boldsymbol{A}\boldsymbol{X}\|_{F}^{2}$
$\displaystyle=\|\tilde{\boldsymbol{y}}_{1}-\alpha_{1}\tilde{\boldsymbol{x}}_{1}-\beta_{1}\tilde{\boldsymbol{x}}_{2}\|_{2}^{2}+\|\tilde{\boldsymbol{y}}_{n}-b_{n-1}\tilde{\boldsymbol{x}}_{n-1}-\alpha_{n}\tilde{\boldsymbol{x}}_{n}\|_{2}^{2}$
(74)
$\displaystyle+\sum_{i=2}^{n-1}\|\tilde{\boldsymbol{y}}_{i}-b_{i-1}\tilde{\boldsymbol{x}}_{i-1}-\alpha_{i}\tilde{\boldsymbol{x}}_{i}-b_{i}\tilde{\boldsymbol{x}}_{i+1}\|_{2}^{2}.$
(75)

Enforcing stationarity in
$\boldsymbol{c}=\begin{bmatrix}\alpha_{1}&\cdots&\alpha_{n}&\beta_{1}&\cdots&b_{n-1}\end{bmatrix}^{T}$
produces the linear system

$\displaystyle\boldsymbol{T}\boldsymbol{c}=\boldsymbol{d}$ (76)

where $\boldsymbol{T}$ is the symmetric block-tridiagonal matrix

$\displaystyle{\boldsymbol{T}=\left[\begin{array}[]{c|c}\boldsymbol{T}_{1}&\boldsymbol{T}_{2}\\\
\hline\cr\boldsymbol{T}_{2}^{\ast}&\boldsymbol{T}_{3}\end{array}\right]}.$
(79)

where the blocks are

$\displaystyle\boldsymbol{T}_{1}=\mathrm{diag}\left(\|\tilde{\boldsymbol{x}}_{1}\|_{2}^{2}\,,\,\cdots\,,\,\|\tilde{\boldsymbol{x}}_{n}\|_{2}^{2}\right),$
(80)
$\displaystyle\boldsymbol{T}_{2}=\begin{bmatrix}\tilde{\boldsymbol{x}}_{2}^{\ast}\tilde{\boldsymbol{x}}_{1}&\\\
\tilde{\boldsymbol{x}}^{\ast}_{2}\tilde{\boldsymbol{x}}_{1}&\tilde{\boldsymbol{x}}_{3}^{\ast}\tilde{\boldsymbol{x}}_{2}&\\\
&\tilde{\boldsymbol{x}}_{4}^{\ast}\tilde{\boldsymbol{x}}_{3}&\ddots&\\\
&&\ddots&\ddots&\\\
&&&\tilde{\boldsymbol{x}}_{n-1}^{\ast}\tilde{\boldsymbol{x}}_{n-2}&\tilde{\boldsymbol{x}}_{n}^{\ast}\tilde{\boldsymbol{x}}_{n-1}\\\
&&&&\tilde{\boldsymbol{x}}_{n}^{\ast}\tilde{\boldsymbol{x}}_{n-1}\end{bmatrix},$
(81)
$\displaystyle\boldsymbol{T}_{3}=\begin{bmatrix}\tilde{\boldsymbol{x}}_{2}^{\ast}\tilde{\boldsymbol{x}}_{2}+\tilde{\boldsymbol{x}}_{1}^{\ast}\tilde{\boldsymbol{x}}_{1}&\tilde{\boldsymbol{x}}_{3}^{\ast}\tilde{\boldsymbol{x}}_{1}\\\
\tilde{\boldsymbol{x}}_{3}^{\ast}\tilde{\boldsymbol{x}}_{1}&\tilde{\boldsymbol{x}}_{3}^{\ast}\tilde{\boldsymbol{x}}_{3}+\tilde{\boldsymbol{x}}_{2}^{\ast}\tilde{\boldsymbol{x}}_{2}&\tilde{\boldsymbol{x}}_{4}^{\ast}\tilde{\boldsymbol{x}}_{2}\\\
&\tilde{\boldsymbol{x}}_{4}^{\ast}\tilde{\boldsymbol{x}}_{2}&\ddots&\ddots&&\\\
&&\ddots&\ddots&\tilde{\boldsymbol{x}}_{n}^{\ast}\tilde{\boldsymbol{x}}_{n-1}\\\
&&&\tilde{\boldsymbol{x}}_{n}^{\ast}\tilde{\boldsymbol{x}}_{n-1}&\tilde{\boldsymbol{x}}_{n}^{\ast}\tilde{\boldsymbol{x}}_{n-1},\end{bmatrix}$
(82)

and

$\displaystyle\boldsymbol{d}=\begin{bmatrix}\tilde{\boldsymbol{x}}_{1}^{\ast}\tilde{\boldsymbol{y}}_{1}&\cdots&\tilde{\boldsymbol{x}}_{n}^{\ast}\tilde{\boldsymbol{y}}_{n}&&\tilde{\boldsymbol{x}}_{1}^{\ast}\tilde{\boldsymbol{y}}_{2}+\tilde{\boldsymbol{x}}_{2}^{\ast}\tilde{\boldsymbol{y}}_{1}&\cdots&\tilde{\boldsymbol{x}}_{n-1}^{\ast}\tilde{\boldsymbol{y}}_{n}+\tilde{\boldsymbol{x}}_{n}^{\ast}\tilde{\boldsymbol{y}}_{n-1}\end{bmatrix}^{\ast}.$

Note that $\boldsymbol{T}_{1}\in\mathbb{R}^{n\times n}$,
$\boldsymbol{T}_{2}\in\mathbb{R}^{n\times(n-1)}$, and
$\boldsymbol{T}_{3}\in\mathbb{R}^{(n-1)\times(n-1)}$. The system (76) can be
solved very efficiently by exploiting the sparsity and symmetry of
$\boldsymbol{T}$.

###  C.4 Local total least squares

This estimator is unbiased if there is only noise in $\boldsymbol{Y}$. If
$\boldsymbol{X}$ and $\boldsymbol{Y}$ are both contaminated by noise then (29)
can be suitably adapted and solved via total least squares [97], provided that
the noise is normally distributed and has the same (diagonal) covariance
matrix for $\boldsymbol{X}$ and $\boldsymbol{Y}$.

###  C.5 Weaker spatial locality

A weaker version of spatial locality can also be enforced by penalising
entries via

$\displaystyle\operatorname*{arg\rm{}min}_{\boldsymbol{A}}\|\boldsymbol{Y}-\boldsymbol{A}\boldsymbol{X}\|_{F}+\lambda\|\boldsymbol{H}\odot\boldsymbol{A}\|_{F}$
(83)

where $\boldsymbol{H}$ is a matrix that penalises the solution for grid points
that are spatially far from one another and $\lambda$ is a regularization
parameter. For example, a suitable choice for $\boldsymbol{H}$ could be the
Gaussian kernel matrix

$\displaystyle\boldsymbol{H}_{i,j}=\exp\left(\|\xi_{i}-\xi_{j}\|_{2}/(2\sigma^{2})\right)$
(84)

for some parameter $\sigma$, where $\\{\xi_{i}\\}$ are the known grid points.
For general $\boldsymbol{H}$, the solution to (83) is

$\displaystyle\tilde{\boldsymbol{A}}_{i}=\tilde{\boldsymbol{y}}_{i}\boldsymbol{X}^{\ast}\left(\boldsymbol{X}\boldsymbol{X}^{\ast}+\textnormal{diag}(\tilde{\boldsymbol{h}}_{i})\right)^{-1}$

for $1\leq i\leq n$ where $\tilde{\boldsymbol{h}}_{i}$ is the $i$-th row of
$\boldsymbol{H}$. Unfortunately, since the regularizer penalises every entry
of $\boldsymbol{A}$, the implementation of (C.5) costs $\mathcal{O}(n^{4})$
operations to compute the full model $\boldsymbol{A}$. Thus, we recommend
pursuing inexact but efficient solutions, possibly using descent methods
[104].

##  Appendix D Further details for causal systems

###  D.1 Updating equations for upper triangular system

Here we present an efficient implementation of the solution (35). Consecutive
blocks of $\boldsymbol{X}$ (e.g. $\boldsymbol{X}_{j-1:n,:}$ and
$\boldsymbol{X}_{j:n,:}$) are related by the rank-1 update

$\displaystyle\boldsymbol{X}_{j-1:n,:}=\begin{bmatrix}\boldsymbol{0}\\\
\boldsymbol{X}_{j:n,:}\end{bmatrix}+\begin{bmatrix}1\\\
\boldsymbol{0}\end{bmatrix}\tilde{\boldsymbol{x}}_{j-1}.$ (85)

This observation enables us to efficiently calculate the pseudoinverse
$\tilde{\boldsymbol{X}}_{j-1:n,:}^{\dagger}$ from
$\tilde{\boldsymbol{X}}_{j:n,:}^{\dagger}$. The following formulas come from
[105] (after some simplification). Let
$\hat{\boldsymbol{x}}_{j-1}=\tilde{\boldsymbol{x}}_{j-1}\left(\boldsymbol{I}-\boldsymbol{X}_{j:n,:}^{\dagger}\boldsymbol{X}_{j:n,:}\right)$
be the orthogonal projection of $\tilde{\boldsymbol{x}}_{j}$ onto the
complement of the row space of $\boldsymbol{X}_{j:n,:}$. If
$\tilde{\boldsymbol{x}}_{j-1}$ lies in the row space of
$\boldsymbol{X}_{j:n,:}$ (i.e. $\|\hat{\boldsymbol{x}}_{j-1}\|_{2}\neq 0$)
then

$\displaystyle\boldsymbol{X}_{j-1:n,:}^{\dagger}=\begin{bmatrix}\hat{\boldsymbol{x}}_{j-1}^{\dagger}&&\left(\boldsymbol{I}-\hat{\boldsymbol{x}}_{j-1}^{\dagger}\tilde{\boldsymbol{x}}_{j-1}\right)\boldsymbol{X}_{j:n,:}^{\dagger}\end{bmatrix}$
(86)

otherwise,

$\displaystyle\boldsymbol{X}_{j-1:n,:}^{\dagger}=\begin{bmatrix}\boldsymbol{0}&\boldsymbol{X}_{j:n,:}^{\dagger}\end{bmatrix}\left(\boldsymbol{I}-\boldsymbol{v}_{j}^{\dagger}\boldsymbol{v}_{j}\right)$
(87)

where
$\boldsymbol{v}_{j}=\begin{bmatrix}1&-\tilde{\boldsymbol{x}}_{j-1}\boldsymbol{X}_{j:n,:}^{\dagger}\end{bmatrix}$.
Note that (86) and (87) can be formed using only matrix-vector products. It
turns out that the above formulation is unstable for realistic data
measurements (with high condition number). Below, we present an alternative
algorithm that is as fast but more stable.

###  D.2 Alternative (more) stable solution

We now derive an alternative solution to the upper-triangular piDMD problem of
section 4.5. In our experience, this version is more stable than that of the
previous section. We use the economy RQ decomposition of $\boldsymbol{X}$ and
write $\boldsymbol{X}=\boldsymbol{R}\boldsymbol{Q}$ where
$\boldsymbol{R}\in\mathbb{C}^{n\times n}$ is upper triangular, and
$\boldsymbol{Q}\in\mathbb{C}^{n\times m}$ satisfies
$\boldsymbol{Q}\boldsymbol{Q}^{\ast}=\boldsymbol{I}$. Note that

$\displaystyle\|\boldsymbol{Y}-\boldsymbol{A}\boldsymbol{X}\|_{F}^{2}$
$\displaystyle=\|\boldsymbol{Y}\|_{F}^{2}-\|\boldsymbol{Y}\boldsymbol{Q}^{\ast}\|_{F}^{2}+\|\boldsymbol{Y}\boldsymbol{Q}^{\ast}-\boldsymbol{A}\boldsymbol{R}\|_{F}^{2}.$
(88)

The first two terms in (88) are independent of $\boldsymbol{A}$ and, by the
submultiplicativity of the Frobenius norm222In particular,
$\|\boldsymbol{C}\boldsymbol{D}\|_{F}\leq\|\boldsymbol{C}\|_{F}\|\boldsymbol{D}\|_{2}$
for any matrices $\boldsymbol{C}$ and $\boldsymbol{D}$. Note the $2$-norm in
the second term on the right., have a non-negative sum. Thus, the upper-
triangular Procrustes problem may be phrased as

$\displaystyle\operatorname*{arg\rm{}min}_{\boldsymbol{A}\in\mathcal{M}}\|\boldsymbol{Y}-\boldsymbol{A}\boldsymbol{X}\|_{F}=\operatorname*{arg\rm{}min}_{\boldsymbol{A}\in\mathcal{M}}\|\boldsymbol{Y}\boldsymbol{Q}^{\ast}-\boldsymbol{A}\boldsymbol{R}\|_{F}.$
(89)

Now, since the product of upper triangular matrices is upper triangular,
$\boldsymbol{A}\boldsymbol{R}$ is also upper triangular. Thus, the lower
triangular component of $\boldsymbol{Y}\boldsymbol{Q}^{\ast}$ is irrelevant to
the minimization problem and we may write the problem as

$\displaystyle\operatorname*{arg\rm{}min}_{\boldsymbol{A}\in\mathcal{M}}\|\texttt{triu}\left(\boldsymbol{Y}\boldsymbol{Q}^{\ast}\right)-\boldsymbol{A}\boldsymbol{R}\|_{F},$
(90)

where, in Matlab notation, $\texttt{triu}(\cdot)$ extracts the upper
triangular portion of its argument. Since the rows of $\boldsymbol{A}$ are
independent, the Frobenius norm may be expanded row-wise and minimized
separately. Thus, the rows of $\boldsymbol{A}$ solve the $n$ smaller
minimization problems

$\displaystyle\operatorname*{arg\rm{}min}_{{\boldsymbol{A}}_{i,:}}\|\left(\boldsymbol{Y}\boldsymbol{Q}^{\ast}\right)_{i,i:n}-\boldsymbol{A}_{i,i:n}{\boldsymbol{R}}_{i}\|_{F}\qquad\qquad\textnormal{for
}1\leq i\leq n,$ (91)

where $\boldsymbol{R}_{i}=\boldsymbol{R}_{i:n,i:n}\in\mathbb{R}^{(n-i+1)\times
n-i+1}$ is the $i$-th bottom right block of $\boldsymbol{R}$ and
$\boldsymbol{A}_{i,i:n}\in\mathbb{R}^{n-i+1}$ is the $i$-th row of
$\boldsymbol{A}$ with the $i-1$ leading zeros removed. Each minimization
admits a unique solution if each $\boldsymbol{R}_{i}$ is full rank, which
occurs if the pivots of $\boldsymbol{R}$ are non zero, i.e. $\boldsymbol{X}$
is full rank. If the data is rank deficient, as is often the case in
applications, there are infinitely many solutions. In any case, we seek the
smallest solution (in the Frobenius norm sense), which is given by

$\displaystyle\boldsymbol{A}_{i,i:n}=\left(\boldsymbol{Y}\boldsymbol{Q}^{\ast}\right)_{i,i:n}{\boldsymbol{R}}_{i}^{\dagger}.$
(92)

Calculating $\boldsymbol{A}$ via (92) naively requires computing each
pseudoinverse $n$ times for a total of $\mathcal{O}(n^{4})$ operations. As in
§D.1, this complexity can be lowered with an updating procedure: following
[105], we can efficiently compute ${\boldsymbol{R}}_{i}^{\dagger}$ using
${\boldsymbol{R}}_{i+1}^{\dagger}$. Thus, we iterate through (92) backwards,
starting at $i=n$. when the (scalar) pseudoinverse is given by
${\boldsymbol{R}}_{n,n}^{\dagger}=1/{\boldsymbol{R}}_{n,n}$. Proceeding
inductively, we note that the next matrix block can be expressed as a rank-one
update to the previous block:

$\displaystyle{\boldsymbol{R}}_{i}=\begin{bmatrix}0&\boldsymbol{0}\\\
\boldsymbol{0}&{\boldsymbol{R}}_{i+1}\end{bmatrix}+\begin{bmatrix}1\\\
\boldsymbol{0}\end{bmatrix}{\boldsymbol{R}}_{i,i:n}.$ (93)

The updating formula takes one of two forms depending on whether the pivot
${\boldsymbol{R}}_{i,i}$ is smaller than some tolerance
$\epsilon_{\textrm{tol}}$ close to machine precision:

##### Case 1: $|{\boldsymbol{R}}_{i,i}|>\epsilon_{\textrm{tol}}$

$\displaystyle{\boldsymbol{R}}_{i}^{\dagger}=\begin{bmatrix}1/{\boldsymbol{R}}_{i,i}&{\boldsymbol{R}}_{i,i:n}^{\ast}{\boldsymbol{R}}_{i+1}^{\dagger}/{\boldsymbol{R}}_{i,i}\\\
\boldsymbol{0}&{\boldsymbol{R}}_{i+1}^{\dagger}\end{bmatrix}.$ (94)

##### Case 2: $|{\boldsymbol{R}}_{i,i}|<\epsilon_{\textrm{tol}}$

$\displaystyle{\boldsymbol{R}}_{i}^{\dagger}=\begin{bmatrix}0&\boldsymbol{0}\\\
\boldsymbol{0}&{\boldsymbol{R}}_{i+1}^{\dagger}\end{bmatrix}\left(\boldsymbol{I}-\frac{\boldsymbol{\alpha}_{i}\boldsymbol{\alpha}_{i}^{\ast}}{\|\boldsymbol{\alpha}_{i}\|_{2}^{2}}\right)$
(95)

where
$\boldsymbol{\alpha}_{i}=({\boldsymbol{R}}_{i,i:n}{\boldsymbol{R}}_{i+1}^{\dagger})^{\ast}-\begin{bmatrix}1&\boldsymbol{0}\end{bmatrix}^{T}$.
While mathematically correct, the accuracy of (94) and (95) can deteriorate
over many iterations. Thus, in applications it may be sensible to reform the
pseudoinverse ${\boldsymbol{R}}_{i}^{\dagger}$ from scratch every so often.

The eigenvalues of $\boldsymbol{A}$ can also be computed without forming
$\boldsymbol{A}$ explicitly. Note that if ${\boldsymbol{R}_{i,i}}\neq 0$ then
the upper-triangular structure of ${\boldsymbol{R}}_{i}$ implies that
$\boldsymbol{A}_{i,i}=\left(\boldsymbol{Y}\boldsymbol{Q}^{\ast}\right)_{i,i}/\boldsymbol{R}_{i,i}$.
Further, the eigenvalues of a triangular matrix are its diagonal entries, so
$\\{\boldsymbol{A}_{i,i}\\}$ are the eigenvalues of $\boldsymbol{A}$ (for
$\boldsymbol{R}_{i,i}\neq 0$) and can be computed without forming
$\boldsymbol{A}$ fully.

## References

  * [1] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, “Physics-informed machine learning,” Nat. Rev. Phys., vol. 3, no. 6, pp. 422–440, 2021.
  * [2] J. Behler and M. Parrinello, “Generalized neural-network representation of high-dimensional potential-energy surfaces,” Phys. Rev. Lett., vol. 98, no. 14, p. 146401, 2007.
  * [3] S. Greydanus, M. Dzamba, and J. Yosinski, “Hamiltonian neural networks,” in Adv. Neural Inf. Process. Syst., vol. 32, 2019.
  * [4] M. Reichstein, G. Camps-Valls, B. Stevens, M. Jung, J. Denzler, N. Carvalhais, and Prabhat, “Deep learning and process understanding for data-driven Earth system science,” Nature, vol. 566, no. 7743, pp. 195–204, 2019.
  * [5] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” J. Comput. Phys., vol. 378, pp. 686–707, 2019.
  * [6] M. Cranmer, S. Greydanus, S. Hoyer, P. Battaglia, D. Spergel, and S. Ho, “Lagrangian neural networks,” arXiv preprint arXiv:2003.04630, 2020\.
  * [7] M. Raissi, A. Yazdani, and G. E. Karniadakis, “Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations,” Science, vol. 367, no. 6481, pp. 1026–1030, 2020.
  * [8] R. Wang, D. Maddix, C. Faloutsos, Y. Wang, and R. Yu, “Bridging physics-based and data-driven modeling for learning dynamical systems,” Proc. Mach. Learn. Res., vol. 144, pp. 385–398, 2020.
  * [9] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis, “Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators,” Nat. Mach. Intell., vol. 3, no. 3, pp. 218–229, 2021.
  * [10] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, “Fourier Neural Operator for Parametric Partial Differential Equations,” in 9th Int. Conf. Learn. Represent., 2021.
  * [11] P. J. Baddoo, B. Herrmann, B. J. McKeon, and S. L. Brunton, “Kernel Learning for Robust Dynamic Mode Decomposition: Linear and Nonlinear Disambiguation Optimization (LANDO),” arXiv preprinv arXiv:2106.01510, 2021.
  * [12] S. Klus, P. Gelß, F. Nüske, and F. Noé, “Symmetric and antisymmetric kernels for machine learning problems in quantum physics and chemistry,” Mach. Learn.: Sci. Technol., vol. 2, no. 4, 2021.
  * [13] V. Shah, A. Joshi, S. Ghosal, B. Pokuri, S. Sarkar, B. Ganapathysubramanian, and C. Hegde, “Encoding Invariances in Deep Generative Models,” arXiv preprint arXiv:1906.01626, 2019.
  * [14] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” Proc. Natl. Acad. Sci., vol. 113, no. 15, pp. 3932–3937, 2016\.
  * [15] S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Data-driven discovery of partial differential equations,” Sci. Adv., vol. 3, no. 4, p. e1602614, 2017.
  * [16] J. C. Loiseau and S. L. Brunton, “Constrained sparse Galerkin regression,” J. Fluid Mech., vol. 838, pp. 42–67, 2018.
  * [17] Y. Guan, S. L. Brunton, and I. Novosselov, “Sparse nonlinear models of chaotic electroconvection,” Royal Society Open Science, vol. 8, no. 8, p. 202367, 2021.
  * [18] L. Zanna and T. Bolton, “Data-driven equation discovery of ocean mesoscale closures,” Geophysical Research Letters, vol. 47, no. 17, p. e2020GL088376, 2020.
  * [19] P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” J. Fluid Mech., vol. 656, pp. 5–28, 2010.
  * [20] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson, “Spectral analysis of nonlinear flows,” J. Fluid Mech., vol. 641, pp. 115–127, 2009.
  * [21] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz, “On dynamic mode decomposition: Theory and applications,” J. Comput. Dyn., vol. 1, no. 2, pp. 391–421, 2014.
  * [22] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems.  SIAM, 2016.
  * [23] T. Askham and J. N. Kutz, “Variable projection methods for an optimized dynamic mode decomposition,” SIAM J. Appl. Dyn. Syst., vol. 17, no. 1, pp. 380–416, 2018.
  * [24] P. J. Schmid, “Dynamic Mode Decomposition and Its Variants,” Annu. Rev. of Fluid Mech., vol. 54, pp. 225–254, 2022.
  * [25] J. L. Proctor and P. A. Eckhoff, “Discovering dynamic patterns from infectious disease data using dynamic mode decomposition,” International health, vol. 7, no. 2, pp. 139–145, 2015.
  * [26] E. Berger, M. Sastuba, D. Vogt, B. Jung, and H. B. Amor, “Estimation of perturbations in robotic behavior using dynamic mode decomposition,” Journal of Advanced Robotics, vol. 29, no. 5, pp. 331–343, 2015.
  * [27] I. Abraham and T. D. Murphey, “Active learning of dynamics for data-driven control using Koopman operators,” IEEE Transactions on Robotics, vol. 35, no. 5, pp. 1071–1083, 2019.
  * [28] B. W. Brunton, L. A. Johnson, J. G. Ojemann, and J. N. Kutz, “Extracting spatial–temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition,” Journal of Neuroscience Methods, vol. 258, pp. 1–15, 2016.
  * [29] A. Goldschmidt, E. Kaiser, J. L. DuBois, S. L. Brunton, and J. N. Kutz, “Bilinear dynamic mode decomposition for quantum control,” New J. Phys., vol. 23, no. 3, p. 033035, 2021.
  * [30] Y. Susuki, I. Mezić, and T. Hikihara, “Coherent Swing Instability of Power Grids,” Journal of Nonlinear Science, vol. 21, no. 3, pp. 403–439, 2011.
  * [31] R. Taylor, J. N. Kutz, K. Morgan, and B. A. Nelson, “Dynamic mode decomposition for plasma diagnostics and validation,” Rev. Sci. Instrum., vol. 89, no. 5, p. 053501, 2018.
  * [32] A. A. Kaptanoglu, K. D. Morgan, C. J. Hansen, and S. L. Brunton, “Characterizing magnetized plasmas with dynamic mode decomposition,” Physics of Plasmas, vol. 27, p. 032108, 2020.
  * [33] S. Bagheri, “Effects of weak noise on oscillating flows: Linking quality factor, Floquet modes, and Koopman spectrum,” Physics of Fluids, vol. 26, no. 9, p. 094104, 2014.
  * [34] S. T. Dawson, M. S. Hemati, M. O. Williams, and C. W. Rowley, “Characterizing and correcting for the effect of sensor noise in the dynamic mode decomposition,” Exp. Fluids, vol. 57, no. 3, p. 42, 2016.
  * [35] M. S. Hemati, C. W. Rowley, E. A. Deem, and L. N. Cattafesta, “De-biasing the dynamic mode decomposition for applied Koopman spectral analysis,” TCFD, vol. 31, no. 4, pp. 349–368, 2017.
  * [36] J. P. Cunningham and Z. Ghahramani, “Linear dimensionality reduction: Survey, insights, and generalizations,” The Journal of Machine Learning Research, vol. 16, no. 1, pp. 2859–2900, 2015.
  * [37] J. R. Hurley and R. B. Cattell, “The Procrustes Program: Producing direct rotation to test a hypothesized factor structure,” Behav. Sci., vol. 7, no. 2, pp. 258–262, 1962.
  * [38] I. Söderkvist, “Perturbation analysis of the orthogonal procrustes problem,” BIT, vol. 33, no. 4, pp. 687–694, 1993.
  * [39] N. T. Trendafilov and R. A. Lippert, “The multimode Procrustes problem,” Linear Algebra Appl., vol. 349, no. 1-3, pp. 245–264, 2002.
  * [40] J. C. Gower and G. B. Dijksterhuis, Procrustes Problems.  Oxford University Press, 2007.
  * [41] N. J. Higham, “The symmetric Procrustes problem,” BIT, vol. 28, no. 1, pp. 133–143, 1988.
  * [42] J. M. ten Berge, H. A. Kiers, and J. J. Commandeur, “Orthogonal Procrustes rotation for matrices with missing values,” Br. J. Math. Stat. Psychol., vol. 46, no. 1, pp. 119–134, 1993.
  * [43] K. S. Arun, “A Unitarily Constrained Total Least Squares Problem in Signal Processing,” SIAM J. Matrix Anal. Appl., vol. 13, no. 3, pp. 729–745, 1992.
  * [44] N. Gillis and P. Sharma, “A semi-analytical approach for the positive semidefinite Procrustes problem,” Linear Algebra Appl., vol. 540, pp. 112–137, 2018.
  * [45] P. H. Schönemann, “A generalized solution of the orthogonal procrustes problem,” Psychometrika, vol. 31, no. 1, pp. 1–10, 1966.
  * [46] T. Pumir, A. Singer, and N. Boumal, “The generalized orthogonal Procrustes problem in the high noise regime,” Inf. Inference A J. IMA, vol. 10, no. 3, pp. 921–954, 2021.
  * [47] T. Viklands, Algorithms for the Weighted Orthogonal Procrustes Problem and other Least Squares Problems.  PhD thesis, Umeå University, 2006.
  * [48] L. Eldén and H. Park, “A Procrustes problem on the Stiefel manifold,” Numer. Math., vol. 82, no. 4, pp. 599–619, 1999.
  * [49] L. E. Andersson and T. Elfving, “A constrained Procrustes problem,” SIAM J. Matrix Anal. Appl., vol. 18, no. 1, pp. 124–139, 1997.
  * [50] P. A. Absil, R. Mahony, and R. Sepulchre, Optimization algorithms on matrix manifolds.  Princeton University Press, 2009.
  * [51] N. Boumal, B. Mishra, P. A. Absil, and R. Sepulchre, “Manopt, a Matlab toolbox for optimization on manifolds,” J. Mach. Learn. Res., vol. 15, no. 42, pp. 1455–1459, 2014.
  * [52] M. R. Jovanović, P. J. Schmid, and J. W. Nichols, “Sparsity-promoting dynamic mode decomposition,” Phys. Fluids, vol. 26, no. 2, p. 024103, 2014\.
  * [53] J. L. Proctor, S. L. Brunton, and J. N. Kutz, “Dynamic mode decomposition with control,” SIAM J. Appl. Dyn. Syst., vol. 15, no. 1, pp. 142–161, 2016.
  * [54] B. R. Noack, W. Stankiewicz, M. Morzynski, and P. J. Schmid, “Recursive dynamic mode decomposition of a transient cylinder wake,” J. Fluid Mech., vol. 809, pp. 843–872, 2016.
  * [55] M. S. Hemati, M. O. Williams, and C. W. Rowley, “Dynamic mode decomposition for large and streaming datasets,” Physics of Fluids (1994-present), vol. 26, no. 11, p. 111701, 2014.
  * [56] H. Zhang, C. W. Rowley, E. A. Deem, and L. N. Cattafesta, “Online dynamic mode decomposition for time-varying systems,” SIAM J. Appl. Dyn. Syst., vol. 18, no. 3, pp. 1586–1609, 2019.
  * [57] S. L. Brunton, J. L. Proctor, J. H. Tu, and J. N. Kutz, “Compressed sensing and dynamic mode decomposition,” Journal of Computational Dynamics, vol. 2, no. 2, pp. 165–191, 2015.
  * [58] F. Gueniat, L. Mathelin, and L. Pastur, “A dynamic mode decomposition approach for large and arbitrarily sampled systems,” Physics of Fluids, vol. 27, no. 2, p. 025113, 2015.
  * [59] J. H. Tu, C. W. Rowley, J. N. Kutz, and J. K. Shang, “Spectral analysis of fluid flows using sub-Nyquist-rate PIV data,” Experiments in Fluids, vol. 55, no. 9, pp. 1–13, 2014.
  * [60] G. H. Golub and C. F. Van Loan, Matrix Computations, vol. 3.  JHU Press, 2013.
  * [61] C. Eckart and G. Young, “The approximation of one matrix by another of lower rank,” Psychometrika, vol. 1, no. 3, pp. 211–218, 1936.
  * [62] I. Mezić, “Analysis of fluid flows via spectral properties of the Koopman operator,” Ann. Rev. Fluid Mech., vol. 45, no. 1, pp. 357–378, 2013.
  * [63] S. L. Brunton, M. Budišić, E. Kaiser, and J. N. Kutz, “Modern Koopman Theory for Dynamical Systems,” arXiv preprint arXiv:2102.12086, 2021.
  * [64] G. Dijksterhuis, “Procrustes analysis in studying sensory-instrumental relations,” Food Qual. Prefer., vol. 5, no. 1-2, pp. 115–120, 1994.
  * [65] M. A. A. Cox and T. F. Cox, “Multidimensional Scaling,” in Handb. Data Vis., pp. 315–347, Springer, Berlin, Heidelberg, 2008.
  * [66] X. Zhang, J. Zou, K. He, and J. Sun, “Accelerating Very Deep Convolutional Networks for Classification and Detection,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 38, no. 10, pp. 1943–1955, 2016.
  * [67] M. B. Richman, “Rotation of principal components,” J. Climatol., vol. 6, no. 3, pp. 293–335, 1986.
  * [68] B. Herrmann, P. J. Baddoo, R. Semaan, S. L. Brunton, and B. J. McKeon, “Data-driven resolvent analysis,” J. Fluid Mech., vol. 918, p. A10, 2021\.
  * [69] E. Noether, “Invariante variationsprobleme,” Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse, vol. 1918, pp. 235–257, 1918.
  * [70] A. Dutt and V. Rokhlin, “Fast Fourier Transforms for Nonequispaced Data,” SIAM J. Sci. Comput., vol. 14, no. 6, pp. 1368–1393, 1993.
  * [71] J. F. Gibson, “Channelflow: A spectral Navier-Stokes simulator in C++,” tech. rep., U. New Hampshire, 2014.  Channelflow.org.
  * [72] M. Ilak and C. W. Rowley, “Modeling of transitional channel flow using balanced proper orthogonal decomposition,” Phys. Fluids, vol. 20, p. 034103, 2008.
  * [73] N. J. Higham, Functions of matrices: theory and computation.  Society for Industrial and Applied Mathematics, 2008.
  * [74] D. J. Griffiths, Introduction to quantum mechanics.  Pearson Prentice Hall, 1995.
  * [75] D. T. Smithey, M. Beck, M. G. Raymer, and A. Faridani, “Measurement of the Wigner distribution and the density matrix of a light mode using optical homodyne tomography: Application to squeezed states and the vacuum,” Phys. Rev. Lett., vol. 70, no. 9, pp. 1244–1247, 1993.
  * [76] K. Vogel and H. Risken, “Determination of quasiprobability distributions in terms of probability distributions for the rotated quadrature phase,” Phys. Rev. A, vol. 40, no. 5, pp. 2847–2849, 1989.
  * [77] J. S. Lundeen, B. Sutherland, A. Patel, C. Stewart, and C. Bamber, “Direct measurement of the quantum wavefunction,” Nature, vol. 474, no. 7350, pp. 188–191, 2011.
  * [78] S. L. Brunton and J. N. Kutz, Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control.  Cambridge University Press, 2019.
  * [79] B. J. McKeon and A. S. Sharma, “A critical-layer framework for turbulent pipe flow,” J. Fluid Mech., vol. 658, p. 336–382, 2010.
  * [80] M. R. Jovanović, “From bypass transition to flow control and data-driven turbulence modeling: An input–output viewpoint,” Ann. Rev. Fluid Mech., vol. 53, pp. 311–345, 2021.
  * [81] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll, “Hydrodynamic stability without eigenvalues,” Science, vol. 261, no. 5121, pp. 578–584, 1993.
  * [82] T. Kato, Perturbation Theory for Linear Operators.  Classics in Mathematics, Berlin, Heidelberg: Springer Berlin Heidelberg, 1966.
  * [83] S. C. Reddy and L. N. Trefethen, “Pseudospectra of the convection-diffusion operator,” SIAM J. Appl. Math., vol. 54, no. 6, pp. 1634–1649, 1994.
  * [84] L. N. Trefethen and M. Embree, Spectra and pseudospectra: the behavior of nonnormal matrices and operators.  Princeton University Press, 2005.
  * [85] T. A. Driscoll, N. Hale, and L. N. Trefethen, Chebfun Guide.  Oxford: Pafnuty Publications, 2014.
  * [86] J. L. Aurentz and L. N. Trefethen, “Block operators and spectral discretizations,” SIAM Rev., vol. 59, no. 2, pp. 423–446, 2017.
  * [87] S. L. Brunton, B. W. Brunton, J. L. Proctor, E. Kaiser, and J. Nathan Kutz, “Chaos as an intermittently forced linear system,” Nat. Commun., vol. 8, no. 1, p. 19, 2017.
  * [88] P. Krishnamurthy and F. Khorrami, “Generalized state scaling and applications to feedback, feedforward, and nontriangular nonlinear systems,” IEEE Trans. Automat. Contr., vol. 52, no. 1, pp. 102–108, 2007.
  * [89] A. R. Teel, “A nonlinear small gain theorem for the analysis of control systems with saturation,” IEEE Trans. Automat. Contr., vol. 41, no. 9, pp. 1256–1270, 1996.
  * [90] A. M. Annaswamy and J. Baillieul, “Adaptive Control of Nonlinear Systems with a Triangular Structure,” IEEE Trans. Automat. Contr., vol. 39, no. 7, pp. 1411–1428, 1994.
  * [91] Z. Jackiewicz, M. Rahman, and B. D. Welfert, “Numerical solution of a Fredholm integro-differential equation modelling $\dot{\theta}$-neural networks,” Appl. Math. Comput., vol. 195, no. 2, pp. 523–536, 2008.
  * [92] J. Medlock and M. Kot, “Spreading disease: Integro-differential equations old and new,” Math. Biosci., vol. 184, no. 2, pp. 201–222, 2003.
  * [93] N. J. Higham, “Matrix nearness problems and applications,” Appl. Matrix Theory, pp. 1–27, 1989.
  * [94] A. Zygmund, Trigonometric series, vol. I, II.  Cambridge University Press, 1988.
  * [95] D. Ruiz-Antolin and A. Townsend, “A nonuniform fast fourier transform based on low rank approximation,” SIAM J. Sci. Comput., vol. 40, no. 1, pp. A529–A547, 2018.
  * [96] G. H. Golub and C. F. Van Loan, Matrix computations, vol. 3.  JHU Press, 2012.
  * [97] S. V. Huffel and J. Vandewalle, The total least squares problem: computational aspects and analysis.  SIAM, 1991.
  * [98] J. Yang and Y. B. Deng, “Procrustes problems for general, triangular, and symmetric toeplitz matrices,” J. Appl. Math., 2013.
  * [99] M. R. Jovanović, P. J. Schmid, and J. W. Nichols, “Sparsity-promoting dynamic mode decomposition,” Physics of Fluids, vol. 26, no. 2, p. 024103, 2014.
  * [100] N. Levinson, “The Wiener (Root Mean Square) Error Criterion in Filter Design and Prediction,” J. Math. Phys., vol. 25, no. 1-4, pp. 261–278, 1946\.
  * [101] R. P. Brent, F. G. Gustavson, and D. Y. Yun, “Fast solution of toeplitz systems of equations and computation of Padé approximants,” J. Algorithms, vol. 1, no. 3, pp. 259–295, 1980.
  * [102] M. Stewart, “A superfast toeplitz solver with improved numerical stability,” SIAM J. Matrix Anal. Appl., vol. 25, no. 3, pp. 669–693, 2004\.
  * [103] N. Sedrakyan, “About the applications of one useful inequality,” Kvant J., vol. 97, no. 2, pp. 42–44, 1997.
  * [104] S. Boyd and L. Vandenberghe, Convex optimization.  Cambridge University Press, 2004.
  * [105] C. D. Meyer, Jr., “Generalized Inversion of Modified Matrices,” SIAM J. Appl. Math., vol. 24, no. 3, pp. 315–323, 1973.