Physics-informed dynamic mode decomposition

royalsocietypublishing.org/journal/rspa

Research

Cite this article:

Baddoo PJ, Herrmann B,

McKeonBJ,NathanKutzJ,BruntonSL.2023

Physics-informed dynamic mode

decomposition.

Proc.R.Soc.A

479

: 20220576.

https://doi.org/10.1098/rspa.2022.0576

Received: 1 September 2022

Accepted: 23 January 2023

Subject Areas:

applied mathematics, computational

mathematics, fluid mechanics

Keywords:

machine learning, dynamic mode

decomposition, data-driven dynamical

systems

Author for correspondence:

Steven L. Brunton

e-mail:

sbrunton@uw.edu

Electronic supplementary material is available

online at

https://doi.org/10.6084/m9.figshare.

c.6423942

Physics-informed dynamic

mode decomposition

Peter J. Baddoo

, Benjamin Herrmann

Beverley J. McKeon

,J.NathanKutz

and

Steven L. Brunton

Department of Mathematics, Massachusetts Institute of

Technology, Cambridge, MA 02139, USA

Department of Mechanical Engineering, University of Chile,

Beauchef 851, Santiago, Chile

Graduate Aerospace Laboratories, California Institute of

Technology, Pasadena, CA 91125, USA

Department of Applied Mathematics, and

Department of

Mechanical Engineering, University of Washington, Seattle,

WA 98195, USA

PJB,

0000-0002-8671-6952

;BJM,

0000-0003-4220-1583

;

JNK ,

0000-0002-6004-2275

;SLB,

0000-0002-6565-5118

this

work,

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 can produce 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-equivariance—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

2023 The Authors. Published by the Royal Society under the terms of the

Creative Commons Attribution License

http://creativecommons.org/licenses/

by/4.0/

, which permits unrestricted use, provided the original author and

source are credited.

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

DMD models. We demonstrate piDMD on a range of problems, including energy-preserving

fluid flow, the Schrödinger equation, solute advection-diffusion and three-dimensional

transitional channel flow. In each case, piDMD outperforms standard DMD algorithms in

metrics such as spectral identification, state prediction and estimation of optimal forcings and

responses.

1. Introduction

Integrating partial knowledge of physical principles into data-driven techniques is a primary goal

of the scientific machine learning (ML) community [

]. 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 [

–

], kernel methods [

], deep generative models

[

], data assimilation [

] and sparse regression [

–

]. These examples demonstrate that

incorporating partially known physical principles into ML 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 (DMD) [

–

DMD is a data diagnostic technique that extracts coherent spatial–temporal patterns from

high-dimensional time-series data [

]. Although DMD originated in the fluid dynamics

community [

], the algorithm has since been applied to a wealth of dynamical systems including

in epidemiology [

], robotics [

], neuroscience [

], quantum control [

], power grids [

]

and plasma physics [

]. Despite its widespread successes, DMD is highly sensitive to noise

[

–

], fails to capture travelling wave physics, and can produce overfit models that do not

generalize. We demonstrate that integrating physics into the learning framework can help address

these challenges.

Suppose that we are studying a dynamical system defined by

(

) for an unknown

function

. Given a collection of measurements

{

...

}

, DMD identifies the best low-rank

linear approximation of

. In other words, DMD seeks a low-rank matrix

such that

≈

for

...

. Arranging the data measurements into matrices

and

(see (2.1) for details)

allows us to phrase the above formally as

arg min

rank (

)

≤

−

(1.1)

After approximately

solving (1.1), the DMD process computes the dominant spectral

properties of the learned linear operator [

–

The rank-

constraint in (1.1) is motivated by the assumed modal structure of the system but

does not account for other important physical properties. For example, one limitation of DMD

is that the solution of (1.1) lies within the span of

, so the learned model can fail to generalize

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 in (1.1) by constraining

the solution matrix

to lie on a matrix manifold:

arg min

∈

M

−

(1.2)

Although the problem (1.1) admits an exact solution [

], the vast majority of DMD implementations solve (1.1)

approximately

by projecting the problem into the

-dimensional subspace spanned by the first

principal components of

.For

that reason, we focus on comparing our new method to the commonly used exact DMD [

] and optimized DMD algorithms

[

] as opposed to the method of [

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

The matrix manifold

is dictated by the known physics of the system at hand. For example, we

may select

such that its members satisfy certain symmetries that are known to be obeyed by

the system. We call (1.2)

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 (1.2) can be computed to give insight into the modes that dominate the system

and their dynamics. Similarly to DMD, these modes can usually be computed without forming

explicitly. The low-rank DMD regression (1.1) is a special case of piDMD where

is the manifold

of matrices with rank bounded by

. However, piDMD models are not always low rank (see,

for example, §4d) 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, as evidenced by a

series of applications in §4. Similarly to standard DMD, our goal is to decompose the data into

its most important modes; as such, our results focus on comparing the spectrum of the learned

operator

with that of the true underlying operator.

The optimization problem (1.2) is well known in the statistical literature as a

Procrustes problem

[

–

], which comprises of finding the optimal transformation between two matrices subject to

certain constraints on the class of admissible transformations. 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,

plays the role of Procrustes’ victim,

is the bed and

is the ‘treatment’

(stretching or amputation).

Procrustes problems (1.2) seek to learn the treatment

that best

relates the data measurements

and

. 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 [

], and symmetric matrices [

In a similar vein to the present work, [

] deployed techniques from group theory and

representation theory to derive a data-driven method to study the Koopman operator for

equivariant dynamical systems. By contrast, piDMD does not work with the Koopman operator

directly but permits a broader class of physical principles, including those that cannot be

expressed as equivariances.

A significant challenge in physics-informed ML is to develop implementations that

incorporate known physics but also scale to higher dimensions [

]. One contribution of this

work is to develop several exact solutions, in terms of standard linear algebra objects, of

physically relevant Procrustes problems. Where possible, we use matrix factorizations and

low-rank representations to alleviate the computational burden of implementation; it is rarely

necessary to form

explicitly when computing its spectral properties or other diagnostics. For

complicated matrix manifolds, an exact solution of the Procrustes problem (1.2) may be intractable

and we must resort to algorithmic approaches. Fortunately, optimization on matrix manifolds is

amaturefield[

], and many algorithms are implemented in the open source software package

‘Manopt’ [

The remainder of the paper is arranged as follows. In §2, we provide background information

on DMD and Procrustes problems. Then, in §3, we describe the broad framework of piDMD. We

consider a range of applications in §4, with a focus on shift-equivariant, conservative, self-adjoint,

local and causal systems. Section 5 concludes with a discussion of the limitations of piDMD

and suggests several future research directions. Further information is available in electronic

supplementary material.

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

snapshots pairs of

features each:

{

(

An open source implementation of piDMD is available in M

ATL AB

www.github.com/baddoo/piDMD

This terminology was first introduced by Hurley & Cattell [

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

...

}

. For example,

may be a discretized fluid flow field at time

and

may be the flow

field at the next time step

. It is convenient to arrange the data into

snapshot matrices of

the form

⎡

⎢

⎣

|| |

···

|| |

⎤

⎥

⎦

⎡

⎢

⎣

⎤

⎥

⎦

and

⎡

⎢

⎣

|| |

···

|| |

⎤

⎥

⎦

⎡

⎢

⎣

⎤

⎥

⎦

(2.1)

so that the

th rows of

and

are the measurements of the

th features and the

th columns are

the

th temporal snapshots. Henceforth, we use

to represent a row vector.

(a) 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 [

]. In particular,

DMD identifies the leading-order spatial eigenmodes of the matrix

in (1.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 [

] and the recent review by

Schmid [

To address the challenges associated with DMD, researchers have derived many variations on

the original algorithm, including sparsity promoting DMD [

], DMD with control [

], noise-

robust variants [

–

], recursive DMD [

], online DMD [

], low-rank DMD [

]and

versions for under-resolved data in space or time [

–

]. At present, the most widely used

variant is ‘exact DMD’ [

], which phrases the DMD solution in terms of the Moore–Penrose

pseudoinverse. Owing 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

and

are connected by an unknown dynamical system of the form

(

); for discrete-time dynamics

and for continuous-time dynamics

.DMD

aims to learn the dominant behaviour of

by finding the best linear approximation for

given

the data and then performing diagnostics on that approximation. Thus, DMD seeks the linear

operator

that best maps the snapshots in the set

{

}

to those in the set

{

}

≈

for

...

(2.2)

Expressed in terms of the snapshot matrices in (2.1), the linear system in (2.2) becomes

≈

(2.3)

and the optimization problem for

is given by (1.1). The minimum-norm solution for

is given

†

∗

(2.4)

where † indicates the Moore–Penrose pseudoinverse [

]and

∗

is the singular value

decomposition.

In many applications, the state dimension

is very large, and forming or storing

explicitly

becomes impractical. Instead, we use a rank-

approximation for

, denoted by

,where

.To

form

, we construct the optimal rank-

approximation for

using the truncated singular value

decomposition [

≈

∗

. We then project

onto the leading

principal components of

∗

−

(2.5)

It is now computationally viable to compute the eigendecomposition of

ΨΛ

(2.6)

The eigenvectors of

can be approximated from the reduced eigenvectors by [

]

−

(2.7)

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

This eigendecomposition is connected to the Koopman operator of the system, and allows

reconstructions and predictions [

]. For example, for a discrete-time system (i.e.

) with evenly spaced samples in time, if the eigenvectors form a linearly independent

set then

ΨΛ

−

(2.8)

where the vector

contains the weights of the modes in the initial condition:

†

. From the

above, it is clear that the eigenvalues

govern the temporal behaviour of the system and the

eigenvectors

are the spatial modes.

(b) Procrustes problems

Procrustes problems (1.2) comprise finding the optimal transformation between two matrices

subject to certain constraints on the class of admissible transformations [

]. The

minimization is usually phrased in the Frobenius norm. Procrustes analysis finds relevance in

a wealth of fields including sensory analysis [

], data visualization [

], neural networks [

climate science [

] and solid mechanics [

]. A summary is available in the monograph [

The earliest and most common formulation is the ‘orthogonal Procrustes problem’ [

Suppose that we have two sets of measurements

and

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,

is constrained to be a unitary matrix

and the minimization problem is

arg min

∗

−

(2.9)

The solution to (2.9) was derived by Schönemann [

]as

∗

(2.10)

where

∗

is a full singular value decomposition. Alternatively,

,where

∗

is a polar decomposition. The solution is unique if and only if

∗

is full-rank. As

we show in §4b, a unitary matrix corresponds to an energy-preserving operator.

There are many solutions for Procrustes problems available in the literature for different matrix

constraints [

]. When exact solutions are not possible, algorithmic solutions can be effective [

3. Physics-informed dynamic mode decomposition

Incorporating physics into ML algorithms involves supplementing an existing technique with

additional biases. Usually, ML practitioners take one of three approaches [

]. First, observational

biases can be embedded through data augmentation techniques; however, augmented data are

not always available, and incorporating additional training samples can become computationally

expensive. Second, physical constraints can be included by suitably penalizing the loss function

[

]. Third, inductive biases can be incorporated directly into the ML architecture in the form of

mathematical constraints [

]. As noted by Karniadakis

et al

], this third approach is arguably

the most principled method since it produces models 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 (1.2).

We exchange the low-rank representation

typically sought in model order reduction in favour of alternative or additional matrix

structures more relevant to the physical problem at hand. In

figure 1

, we illustrate the matrix

structures used in this paper, along with their corresponding physical principles and references

for the optimal solutions.

Analysing a system with the piDMD framework requires four steps: modelling, interpretation,

optimization and diagnostics:

Although the focus of this paper is on the underlying regression problem, we are still primarily concerned with decomposing

the data into the form (2.8); hence, piDMD is a modal decomposition algorithm.

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

physics:

modal

matrix: low-rank

solution:

Schmid [21]

conservative

unitary

Schönemann

[49]

self-adjoint

symmetric

Higham [45]

shift-equivariant

circulant/Toeplitz/Hankel

new solutions

causal

upper triangular

new solutions

local

tri-diagonal

new solutions

multi-dimensional

block matrices

new solutions

dynamic mode decomposition

physics-informed dynamic mode decomposition (piDMD)

Figure 1.

Visual illustrations of the types of matrices, their corresponding physical principles and references to the solutions of

the corresponding optimization problem (1.2).

(i) Modelling: 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 summarize the known or suspected physical properties of the system at

hand.

(ii) Interpretation: 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.

(iii) Optimization: Equipped with a target matrix manifold, we may now solve the relevant

Procrustes problem (1.2). By suitably modifying the optimization routine, we can

guarantee that the resulting model satisfies the physical principle identified in step (i).

(iv) Diagnostics: The final step of piDMD involves extracting physical information from the

learned model

. For example, one may wish to analyse the spectrum or modes, compute

the resolvent modes [

], perform predictions, or investigate other diagnostics. Similarly

to the original DMD formulation, the learned modes can be used to decompose the data

into the form (2.8).

These steps are expanded on in electronic supplementary material, §A. As an example,

consider the travelling wave solution explored in the

figure 2

. We begin with the physical

principle that the system is shift equivariant (modelling), and this lead us to seek a circulant

matrix model (interpretation). We find a solution to the corresponding Procrustes problem

(optimization) and investigated the predictive behaviour of the model (diagnostics).

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

training data

standard DMD:

shift-equivariance:

test data

physics-informed DMD:

true system

physical principle

modal structure:

rank (

)



optimization

predictions

truth

rank (

)



arg min

circulant

arg min

Figure2.

Comparingstandarddynamicmodedecomposition(specifically,optimizedDMD[

])andphysics-informeddynamic

mode decomposition applied to the advection equation with incomplete data. The data are contaminated with 2% additive

Gaussian noise. Having trained DMD models, we then perform predictions for different initial conditions. Standard DMD fails

whereaspiDMDproducesafaithfulprediction.Inthefigure,

representstheshiftoperator.TheperformanceofDMDimproves

ifthedatacoversafullcycleofthepulse.

While each of the four steps can be challenging, the optimization step is often the most difficult.

Moreover, any solutions are 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 (1.2) 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 [

Standard DMD exploits the low-rank structure of

to efficiently perform diagnostics on the

learned model. Some of the manifolds we consider (such as tridiagonal and 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 [

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 [

]. However, the number of equality constraints needed is

(

) so the resulting least-squares matrix will have

(

) 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 §4a and §4c, for example).

DMD uses the available degrees of freedom in the model to optimize the reconstruction of the

data in a low-dimensional subspace. This can result in overfitting since the algorithm does not

respect global properties of the system; as such, the modes usually fail to characterize unseen

test data. In piDMD, we exchange some of the available degrees of freedom in the model to

enforce a global physical structure. This regularization produces models that are suboptimal in

the reconstruction norm but are more generalizable to unseen datasets. The generalizability of

piDMD is observed in each of our test cases as the algorithm is able to more accurately learn the

true eigenvalues of the underlying system. We emphasize that we are not claiming that the low-

rank constraint is always suboptimal. In fact, the low-rank constraint is often best and should be

used as a first attempt when analysing data. However, if prior physical knowledge of the system

is available then a physics-informed manifold may be more appropriate.

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

system

data

standard DMD

physics-informed DMD

physics

(

matrix structure

)

cylinder flow

causal

(

upper triangular

)

conservative

(

unitary

)

Schrödinger equation

two-dimensional

advection

self-adjoint

(

Hermitian

)

spectrum

shift-equivariant

(

block circulant

)

Volterra integro-

differential equation

local

(

tri-diagonal

)

advection-diffusion

spectrum

shift-equivariant

(

block circulant

)

channel flow

= 0

= 12

= 0

–

(

ξ, ν

) u (

ν,

) d

ψψ

(

)(

)

Figure 3.

A comparison of the models learned by exact DMD and piDMD for a range of applications. piDMD identifies the

spectrumofthetrueoperatorwithhigheraccuracythanexactDMD.Thestructureofthemodelmatricesisalsoillustrated;piDMD

models are generally more coherent than those learned by exact DMD. Details of problem sizes are given in the corresponding

sections in the electronic supplementary material. In the spectrum subplots, we plot the true eigenvalues as circles, the DMD

eigenvalues as triangles and the piDMD eigenvalues as crosses.

4. Applyingphysics-informeddynamicmodedecompositiontoenforcecanonical

physical structures

We now present a range of applications of piDMD. Each of the following sections investigates

a detailed application of piDMD for a specific physical principle. A summary is illustrated

figure 3

, and the degrees of freedom in each model, and associated computational cost, is

discussed in electronic supplementary material, F. The physical principles, corresponding matrix

structures, and optimal solutions and extensions are summarized below

§4a: Shift equivariant: circulant (and symmetric, skew-symmetric or unitary, electronic

supplementary material, B.1), low rank (electronic supplementary material, B.2), non-

equally spaced samples (electronic supplementary material, B.3), total least squares

(electronic supplementary material, B.4), Toeplitz & Hankel (electronic supplementary

material, B.5)

§4b: Conservative: unitary (§2b)

§4c: Self-adjoint: symmetric (§4c), skew-symmetric (electronic supplementary material, C.1),

symmetric in a subspace (§4c(i))

§4d: Local: tridiagonal (§4d), variable diagonals (electronic supplementary material, D.1),

periodic (electronic supplementary material, D.2), symmetric tridiagonal (electronic

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

supplementary material, D.3), total least squares (electronic supplementary material,

D.4), regularized locality (electronic supplementary material, D.5).

§4e: Causal: upper triangular (electronic supplementary material, E)

(a) Shift-equivariant systems

Shift equivariance (also called ‘translation equivariance’) is ubiquitous in the physical sciences.

Spatially homogeneous systems—such as constant coefficient PDEs or convolution operators—

are shift equivariant, and thus appear identical with respect to any (stationary) observer. In view

of Noether’s theorem [

], the conserved quantity related to shift equivariance is (generalized)

linear momentum; as such, incorporating shift equivariance 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-equivariant operators.

We define

as the

-shift operator i.e.

(

)

(

) for all test functions

. We say that

a space-continuous linear operator

is shift equivariant if

commutes with the

-shift operator

for all shifts:

(4.1)

An application of the shift operator shows that if

is shift equivariant then the quantity e

−

λξ

is constant for all

. Thus, neglecting boundary conditions for now,

{

λξ

}

are eigenfunctions of

. If the domain is normalized to [

−

1, 1] and has periodic boundaries, then

for integer

. In this case, the corresponding eigenfunctions form an orthogonal basis and the operator is

diagonalized by the known eigenfunctions.

In the following, we assume that we are studying a shift-equivariant system on a periodic

domain

∈

[

−

1, 1]. We move from the continuous formulation to a discretized space and assume

that we are studying a function

(

) and have access to evenly spaced samples of

[

−

...

−

]. We discuss the case of non-equally spaced samples in electronic

supplementary material, B.3. If we define the state variable as

(

)

[

(

...

(

)]

then the discrete-space linear operator

that generates

(

) is diagonalized by its eigenvectors

diag(

)

−

(4.2)

where

−

1)(

−

√

and

−

∗

and

{

}

are the unknown eigenvalues. Equation

(4.2) is equivalent to stating that

is circulant

⎡

⎢

⎣

−

...

−

...

⎤

⎥

⎦

(4.3)

i.e.

(

−

)mod

. If the boundary conditions are instead Dirichlet or Neumann then

is not circulant but Toeplitz; we solve the corresponding Procrustes problem in electronic

supplementary material, B.5.

Substituting (4.2) into (1.2) and noting that the Frobenius norm is invariant to unitary

transformations allows (1.2) to be transformed to

arg min

diag(

)

−

(4.4)

where

∗

and

∗

are the spatial discrete Fourier transforms of

and

.Assuch,

and

can be efficiently formed in

(

log(

)) operations using fast Fourier transform (FFT)

methods [

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

The rows of the cost function (4.4) now decouple to produce

minimization problems:

arg min

−

for 1

≤

(4.5)

where

and

are the

th rows of

and

, respectively. An alternative view of this step is that

spatial wavenumbers decouple in periodic shift-equivariant systems, so we may analyse each

wavenumber individually. The coefficient

represents the temporal growth of the

th Fourier

component. The optimal solution for each eigenvalue follows as:

†

(4.6)

In principle, standard DMD methods [

] could be encouraged to respect shift equivariance

by augmenting the data matrices with their shifted counterparts. However, one would strictly

need to include

shifts, thus producing

samples which is usually too large in applications.

By contrast, the above solution is extremely efficient via the use of the FFT and the decoupling of

wavenumbers.

The electronic supplementary material includes alternative solutions for cases where

the system is shift equivariant and symmetric, skew-symmetric or unitary (electronic

supplementary material, B.1), low rank (electronic supplementary material, B.2), the samples

are not equally spaced (electronic supplementary material, B.3), the total least-squares case

(electronic supplementary material, B.4), and when the system is Toeplitz or Hankel (electronic

supplementary material, B.5).

(i) Example: plane channel flow

We now use this shift-equivariant piDMD to compute the modes of the linearized Navier–Stokes

equations from a single simulation. We consider the incompressible flow inside a plane channel

of size 2

along the spatial

and

coordinates that indicate the streamwise, wall-

normal and spanwise directions, respectively. The configuration considers a Reynolds number of

2000 based on the channel half-height and the centreline velocity, and periodic boundary

conditions in the open faces of the channel; hence the flow is homogeneous in the

-and

directions. We use the spectral code

Channelflow

[

] to perform direct numerical simulations

(DNS) with the same numerical configurations as that presented in [

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

(

,0)

−

(cos(

)

1)e

(

−

)

(4.7)

where

(

−

)

(

−

)

, the parameters are set to

0.7 and

0.6, and the amplitude

of the perturbation was scaled to have an energy-norm of 10

−

to ensure that the effect of

nonlinearity is negligible. This initial condition, which was first studied in [

]andalsoin[

], 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. The data is corrupted with 2% Gaussian

noise.

Our aim is to learn the spectrum of the underlying linear operator that

best

describes the

evolution of the DNS snapshots, and to thus obtain a reliable predictive model for the dynamics.

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 with 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

and

coordinates, and is therefore shift-equivariant in these directions. Hence, we

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

shift-equivariant

directions

data : 200 000 states, 300 snapshots, 2% noise

initial condition

transient growth

decay to laminar

= 0

= 15

= 25

= 50





= 12

= 90

A =

truth

DMD

piDMD

piDMD model

DMD model

learned spectrum

0.02

–0.02

–0.04

–0.06

–0.08

0.02

–0.02

–0.04

–0.06

–0.08

–5

test data

piDMD

DMD

predictions:

0.02

–0.02

–0.04

–0.06

–0.08

–0.1

024



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-equivariant 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.

The training data are corrupted with 2% noise. The predictions of the learned DMD and piDMD modes are also compared. The

isosurfaces are surfaces of constant vertical velocity.

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 to

and

and performing exact DMD on the

-dependent Fourier amplitudes for every streamwise and

spanwise wavenumber tuple.

The spectra learned from both approaches are compared in

figure 4

, showing that piDMD

results in a far more accurate and complete eigenvalue spectrum. Importantly, piDMD accurately

identifies the most dynamically significant (least damped) eigenvalues, as illustrated in the

close-up plot. piDMD identifies three unstable eigenvalues, but the corresponding modes

are of low amplitude. By contrast, exact DMD identifies very few relevant eigenvalues

despite the model’s comparatively high rank of 200. As a result, the predictive ability of

exact DMD is poor while piDMD produces very accurate predictions of unseen test data

(bottom row).

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

(b) 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 [

(

)

=||

. In these variables,

the original optimization problems (1.1), (1.2) equate to finding the model

that minimizes the

energy of the errors between the true and predicted states (

and

, respectively). Thus,

represents a discrete-time linear dynamical system (

) then

is an energy

preserving operator if and only if

(

)

=||

(

)forall

∈

R

(4.8)

In words, (4.8) states that

does not change the energy of the system but merely redistributes

energy between the states. In matrix terminology, (4.8) holds if and only if

unitary

. Therefore,

for conservative systems, the optimization problem (1.2) is the orthogonal Procrustes problem

described in §2b, and the solution is given by (2.10). The eigenvalues of

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 (2.10) requires the full SVD of an

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 (1.1) as a total least-squares problem [

]. Perhaps surprisingly, the solution

to the orthogonal Procrustes problem (2.10) is also the solution to the total least-squares problem

when the solution is constrained to be orthogonal [

]. Thus, the solution (2.10) is statistically

optimal even when there is noise in both

and

, as long as the noise has the same Gaussian

distribution in

and

. The solutions obtained by energy-preserving piDMD and total least-

squares DMD [

] are often similar. Further, there are other optimizations closely related to the

orthogonal Procrustes problem including the case where

lies on the Steifel manifold (i.e. when

is rectangular with orthonormal columns [

]), there is a weighting matrix [

], missing data

[

] or high amounts of noise [

(i) Example:

100 flow past a cylinder

We now apply this energy-preserving piDMD to study the flow past a cylinder, which is a

benchmark problem for modal decompositions [

]. The Reynolds number is

/ν

100,

where

is the diameter of the cylinder,

is the free-stream velocity and

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

449 points. We note that the 2-norm of the

measurements is approximately constant in time. Note that the cylinder flow problem is only

approximately conservative due to the balance of the applied pressure and viscous loss; in reality,

the (numerical) energy of the state vector exhibits small oscillations. 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 exhibits (well-known) 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,

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

noisy training data

POD mode trajectories

time

mode 5

100

–100

–20

–10

–40

–0.05

–0.15

–0.10

–1.0

–0.5

0.5

1.0

mode 7

mode 1

–10

mode 9

010203040506070

truth

DMD

piDMD

learned spectra



Figure5.

LearningthespectrumoftheflowpastacylinderwithnoisymeasurementswithDMDandenergy-preservingpiDMD.

and do not exhibit the spurious energy loss caused by noise. Note that total least-squares DMD

exhibits comparable performance to piDMD in this problem.

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

to lie in the

manifold of symmetric (Hermitian) matrices, such that

∗

. Symmetric matrices have real

eigenvalues and, by the spectral theorem, are diagonalizable. 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 [

]is

∗

(4.9)

where the entries of

are

⎧

⎪

⎨

⎪

⎩

otherwise,

(4.10)

where

∗

and

∗

is the SVD of

,and

is the

th singular value. The

solution of the skew-symmetric case is similar (electronic supplementary material, C.1), and

algorithmic solutions are available when

is instead constrained to be positive definite [

]. For

large-scale systems, a low-rank approximation to

can be obtained by projecting the system onto

the first

principal components of

. This low-rank projection preserves the self-adjointness of

the model.

Note that the symmetric Procrustes problem is connected to the orthogonal Procrustes problem

(§4b) via the fact that every unitary matrix

can be expressed as the exponential of a Hermitian

matrix

exp(

(i) Learning energy states of the Schrödinger equation

Possibly the most famous self-adjoint operator is the quantum Hamiltonian,

[

]. In quantum

mechanics, physical quantities of interest, such as position, momentum, energy and spin, are

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

DMD

(rank 15)

piDMD

(rank 15)

truth

learned energy levels

truth

DMD

piDMD

energy eigenstates

–1

6400 states, 400 snapshots,



aussian noise in

data:

system:

(

)(

)



ΨΨ

Figure 6.

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

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

(

)

(

)

(4.11)

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

(

)

∞

(

−

(4.12)

where

is an energy level (eigenvalue) of

with corresponding eigenstate

, the coefficients

are determined by the initial distribution of the wave function, and

is the spatial coordinate.

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 wave function can be obtained through optical tomographic methods

[

] or sequential measurements of complementary variables [

]. The data matrices

and

are formed from measurements of

(

)

and

)

(

)

, respectively, which are themselves

formed by superposing eigenfunctions obtained by a finite difference method. Thus, the matrices

and

are connected by an unknown self-adjoint Hamiltonian

and we may apply our self-

adjoint piDMD optimization. The measurements of

are contaminated with a large amount

(200%) of Gaussian noise but the measurements of

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. By 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 (

∼

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

,whichis

clean, whereas the eigenstates learned by standard DMD lie in the span of

, which is noisy.

In the electronic supplementary material, lemma C.1 proves that the symmetric piDMD model

(4.9) is less sensitive to noise than the exact DMD model (2.4). In summary, a physics-aware

approach enables a more accurate identification of quantum energy levels and the associated

eigenstates.

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023

royalsocietypublishing.org/journal/rspa

Proc.R.Soc.A

479

:20220576

..........................................................

(d) 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 advection 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 [

]. In fact, the output of the exact DMD algorithm remains invariant to

unitary transformations of 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 [

]. Knowledge of the underlying grid from which the data

were collected enables us to bias the learning process to prioritize relationships that we expect to

have a strong local coupling.

Herein we consider the one-dimensional case; the analysis generalizes straightforwardly to

higher dimensions. We assume that the entries of

are samples of a function

taken at

grid

points

{

}

. The grid points can be arranged sequentially so that

<ξ

. 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

to be larger than entries far from the diagonal. Then the entries of

would satisfy

|≥|

for

−

(4.13)

Equation (4.13) is merely a heuristic, and we would expect it to hold on average rather than for

each (

). Additionally, (4.13) is a difficult and expensive condition to implement in practice.

Instead, we consider an alternative version of spatial locality where we only permit coupling

between states that are sufficiently close

0for

−

(4.14)

Equation (4.14) describes a

-diagonal matrix. If

0 then states can only affect themselves and

is diagonal (comparable to (4.4)). A more interesting case is if

1 and we only allow coupling

between adjacent states. Then

is tridiagonal: the entries of

are zero except the leading, upper

and lower diagonals

⎡

⎢

⎣

−

⎤

⎥

⎦

(4.15)

We now solve the optimization problem (1.2) when

is the manifold of tridiagonal matrices.

Electronic supplementary material, D includes many solutions for more general diagonal-like

structures including longer-range and variable coupling (electronic supplementary material, D.1),

spatially periodic local systems (electronic supplementary material, D.2), self-adjoint (electronic

supplementary material, D.3), total-least squares (electronic supplementary material, D.4) and

weaker local structures (electronic supplementary material, D.5). Since the rows of

are

decoupled, we can row-wise expand the Frobenius norm in (1.2) to obtain

smaller minimization

problems

arg min

−

for 1

≤

(4.16)

with the convention

0and

0

. Recall that

and

are the

th rows of

and

, respectively (2.1). Each minimization problem (4.16) has the minimum-norm solution

⎡

⎢

⎣

−

⎤

⎥

⎦

†

for 1

≤

(4.17)

Each of the

minimizations costs

(

) so the total cost is

(

Downloaded from https://royalsocietypublishing.org/ on 05 May 2023