Dynamic mode decomposition
for compressive system identification
Zhe Bai1∗ , Eurika Kaiser1 , Joshua L. Proctor3 , J. Nathan Kutz2 , Steven L. Brunton1

arXiv:1710.07737v2 [cs.SY] 29 Nov 2017

1

Department of Mechanical Engineering, University of Washington, Seattle, WA 98195, United States
2

Department of Applied Mathematics, University of Washington, Seattle, WA 98195, United States
3

Institute of Disease Modeling, Bellevue, WA 98004, United States

December 1, 2017
Abstract
Dynamic mode decomposition has emerged as a leading technique to identify spatiotemporal coherent structures
from high-dimensional data, benefiting from a strong connection to nonlinear dynamical systems via the Koopman
operator. In this work, we integrate and unify two recent innovations that extend DMD to systems with actuation [56]
and systems with heavily subsampled measurements [17]. When combined, these methods yield a novel framework
for compressive system identification 1 . It is possible to identify a low-order model from limited input–output data and
reconstruct the associated full-state dynamic modes with compressed sensing, adding interpretability to the state of the
reduced-order model. Moreover, when full-state data is available, it is possible to dramatically accelerate downstream
computations by first compressing the data. We demonstrate this unified framework on two model systems, investigating
the effects of sensor noise, different types of measurements (e.g., point sensors, Gaussian random projections, etc.),
compression ratios, and different choices of actuation (e.g., localized, broadband, etc.). In the first example, we explore
this architecture on a test system with known low-rank dynamics and an artificially inflated state dimension. The second
example consists of a real-world engineering application given by the fluid flow past a pitching airfoil at low Reynolds
number. This example provides a challenging and realistic test-case for the proposed method, and results demonstrate
that the dominant coherent structures are well characterized despite actuation and heavily subsampled data.
Keywords– Dynamic mode decomposition, compressed sensing, control theory, nonlinear dynamics, Koopman theory

1

Introduction

cillate at a fixed frequency, and grow or decay at a fixed
rate. Thus, DMD yields a set of modes along with a linear
evolution model. Since being introduced in fluid dynamics, DMD has been widely applied in fields as diverse as
epidemiology [58], neuroscience [12], robotics [9], video
processing [27], and financial trading [49].

Dynamic mode decomposition (DMD) is a dimensionality reduction technique introduced by Schmid [61]
in the fluid dynamics community to decompose highdimensional fluid data into dominant spatiotemporal coherent structures [61, 60, 20, 70, 43, 13, 1]. Shortly after
its introduction, DMD was reframed by Rowley et al. [60]
as a numerical technique to approximate the Koopman
operator [39, 50, 51, 3], establishing a strong connection
to the analysis of nonlinear dynamical systems. Unlike other dimensionality reduction techniques, such as
proper orthogonal decomposition (POD) [32], DMD is
designed to extract modes that are spatially coherent, os1 Code

DMD may be thought of as a form of system identification, resulting in reduced-order models that are more
tractable than the original high-dimensional dynamics.
Low-order models are especially important for the control of high-dimensional dynamical systems, such as fluid
flow control [37, 15], where fast control decisions must
be enacted to reduce latency for robust performance.
There are a number of excellent overviews of model reduction [8] and system identification [45, 44] for such

is publicly available at: https://github.com/zhbai/cDMDc

1

Nomenclature
A, Ã
AY , ÃY
B, B̃, B̂
B̂
b
C
F
G
K
P
q
r, r̃
s
t
∆t
tk
U, Û, Ũ
UY , ÛY , ŨY
u
V, V̂, Ṽ
VY , V̂Y , ṼY
W, WY

State transition matrix for x and x̃
State transition matrix for y
Actuation matrix
Estimate of B
Vector of DMD amplitudes
Output measurement matrix
Discrete cosine transform (DCT)
Augmented matrix of A and B
Number of nonzero DCT coeffs.
Projection matrix acting on x
Spatial coordinate
Rank of truncated SVD of X, Ω
DCT coefficients of x
Time
Time step
kth discrete time step
Left singular vectors (POD modes)
of X, X0 , Ω
Left singular vectors (POD modes)
of Y, Y0 , ΩY
Control input, u ∈ Rq
Right singular vectors of X, X0 , Ω
Right singular vectors of Y, Y0 , ΩY
Eigenvectors of Ã, ÃY

X
X0
X̃
ξ
x
x̃
Y
Y0
y
η
λ
ν
Λ, ΛY
Ω
ΩY
ω
Φ, ΦY
ΦS
Φ̂
Ψ
Σ, Σ̂, Σ̃
ΣY , Σ̂Y , Σ̃Y
Θ
Υ

high-dimensional systems. Balanced truncation provides
a principled approach to reducing the system dimension by identifying a reduced-order model with the most
jointly controllable and observable states [52], and extensions include the balanced proper orthogonal decomposition (BPOD) for systems with very large state dimension [72, 59]. It was recently shown [48] that BPOD
is equivalent to the eigensystem realization algorithm
(ERA) [34]. In Tu et al. [70], it was further shown that
DMD may be equivalent to ERA under certain conditions. Proctor et al. [56], further extended DMD to include inputs and control, disambiguating internal state
dynamics from the effect of actuation.

Data matrix of states, X ∈ Rn×m
Shifted data matrix, X0 ∈ Rn×m
Data matrix of reduced states, X̃ ∈ Rr×m
Spatial variable
State vector, x ∈ Rn
Low-rank state vector, x̃ ∈ Rr
Data matrix of measurements, Y ∈ Rp×m
Shifted matrix of measurements, Y0 ∈ Rp×m
Output vector, y ∈ Rp
Noise magnitude
DMD eigenvalue of A
Wavenumber
Matrix of DMD eigenvalues of A, AY
State and input snapshot matrices [XT ΥT ]T
Output and input snapshot matrices [YT ΥT ]T
Continuous-time DMD eigenvalue of A
Matrix of DMD (DMDc) modes of A, AY
Matrix of sparse representation of Φ
Compressed sensing estimate of Φ
Orthonormal basis (e.g. Fourier or POD)
Matrix of singular values of X, X0 , Ω
Matrix of singular values of Y, Y0 , ΩY
Product of measurement matrix and
sparsifying basis, Θ = CΨ
Data matrix of control inputs

surprisingly few spatial measurements, and if full-state
data is available, performing DMD on compressed data
dramatically reduces computation time. In addition to
compressed sensing, randomized linear algebra has been
leveraged to accelerate DMD computations [26, 10].
In this work, we combine the compressed sensing
DMD [17] and DMD with control [56] approaches, resulting in a powerful mathematical framework for compressive system identification. There has been other work on
compressive system identification [29, 40, 21, 73, 38, 53, 6,
16, 36, 42], but this is the first effort in the context of DMD.
We begin by presenting the DMD algorithm, and related
compressed sensing and control extensions, in a common framework and notation, making it possible to unify
these approaches. Combining these algorithms results
in the compressive DMD with control (cDMDc) architecture. cDMDc then makes it possible to identify reducedorder models on downsampled spatial measurements of
a high-dimensional system, and then reconstruct the fullstate dynamic modes associated with the model using
compressed sensing. This approach adds interpretability
to otherwise black-box system identification models. The
resulting cDMDc architecture is demonstrated on two example dynamical systems, including a high-dimensional
fluid flow simulation.

DMD relies on the fact that even high-dimensional dynamics typically evolve on a low-dimensional attractor.
This low-dimensional behavior suggests sparsity in an
appropriate basis, so that sparsity-promoting and randomized techniques may be exploited to reduce measurement resolution, bandwidth requirements, and computational overhead. Sparsity was first used in DMD by Jovanović et al. [33] to select the dominant DMD modes.
Compressed sensing [19, 24, 7] was subsequently used to
compute DMD using snapshots that were sampled below
the Shannon-Nyquist sampling limit in time [69] and in
space [17, 28]. It was shown in Brunton et al. [17] that
it is possible to reconstruct accurate DMD modes with
2

Algorithm 1 Exact DMD [70]
Input: Data matrix X, shifted data matrix X0 , and target rank r.
Output: DMD spectrum Λ and modes Φ.
1: procedure DMD(X, X0 , r)
2:
[U, Σ, V] ← SVD(X, r)
. Truncated r-rank SVD of X.
3:
Ã ← U∗ X0 VΣ−1
. Low-rank approximation of A.
4:
[W, Λ] ← EIG(Ã)
. Eigen-decomposition of Ã.
5:
Φ ← X0 VΣ−1 W
. DMD modes of A.
6: end procedure
Note: If λi = 0, then φi = Uwi for step 5. In the original DMD algorithm [62] all modes are computed as φi = Uwi .

2

Background: The Dynamic Mode
Decomposition

tk , k = 0, . . . , m are arranged into the data matrix


X = x0

Matrix decomposition techniques are ubiquitous in the
data sciences. Their fundamental objective is often to extract low-rank and interpretable patterns from data. Foremost among matrix decomposition methods is the singular value decomposition (SVD), which is the computational engine for principal component analysis (PCA) and
produces a set of ranked orthonormal modes that capture
the dominant correlation structures in the data. However, the SVD fails to correlate both spatial and temporal features of the data together. The DMD provides a
least-square regression architecture whereby both space
and time are jointly correlated by merging a spatial SVD
with a temporal Fourier transform [20]. Specifically, the
DMD algorithm decomposes the data matrix X into the
rank r approximation
X ≈ ΦΛb

x1

and the time-shifted matrix

X0 = x1

x2

...

...



xm−1 


xm  .

(2)

(3)

Typically n  m, i.e. there are many more spatial measurements available than temporal. The vector xk is the
state of a high-dimensional system such as a fluid flow.
Here, we assume evenly sampled snapshots, although
this is generally not required [70, 2].
The DMD algorithm constructs the leading eigendecomposition of the best-fit operator A, chosen to minimize kxk+1 − Axk k2 over the k = 0, 2, 3, · · · , m − 1 snapshots, so that X0 ≈ AX. Computationally, the matrix A
is obtained as

(1)

where the columns of Φ are the DMD modes (spatial
structures), the elements of the diagonal matrix Λ are
the corresponding DMD eigenvalues (with angular frequency λi = exp(ωi ∆t)), and the vector b determines
the weighting of each of the r modes. Thus, each spatial mode in Φ is associated with a single temporal eigenvalue of Λ.

A = X0 X†

(4)

where A ∈ Rn×n and X† is the Moore-Penrose pseudoinverse of X. The dominant eigenvectors of A are the
dynamic modes Φ, and the associated eigenvalues determine how these modes behave in time.

The DMD method was originally used as a low-rank
diagnostic tool for decomposing fluid flow data into
dominant spatiotemporal modes [61, 60, 70, 43], providing a valuable interpretation of coherent structures in
complex systems. Recent innovations have also made
significant progress to employ DMD for robust futurestate prediction [33, 2] and control for input-output systems [56], even when using only a small number of measurements [17].

In practice, the high-dimensional A is not computed
directly. Instead, the SVD can be used to first project
to a low-rank subspace and then compute the matrix
Ã = U∗ AU which has many of the same eigenvalues as
A. This provides an efficient algorithm whose computational expense is bounded by the rank r of the data. The
exact DMD algorithm of Tu et al. [70] is shown in Algorithm 1. Using a variable projection optimization scheme,
the exact DMD method can be modified to handle arbitrary temporal spacing between snapshots. It further produces a more robust, or optimal DMD approximation, for
noisy data [2].

The exact DMD algorithm formulates the decomposition as a least-squares regression. Specifically, a series of
snapshots xk ∈ Rn sampled at discrete instances in time
3

Algorithm 2 DMD with control [56]
Input: Data matrices X, X0 , input snapshot matrix Υ, target rank r of X or X0 and r̃ of Ω.
Optional: Actuation matrix B.
Output: Spectrum Λ, modes Φ, [ and actuation matrix B̂].
. Optional outputs in brackets [·].
1: procedure DMD C (X, X0 , Υ, r, r̃, [B])
2:
if B is known then
3:
[Λ, Φ] ← DMD(X, X0 − BΥ, r)
. Perform DMD (Algorithm 1) adjusted
4:
. for known actuation.
5:
else
 
X
6:
Ω←
. Matrix of the state and input snapshots.
Υ
7:
[Ũ, Σ̃, Ṽ] ← SVD(Ω, r̃)
. Truncated r̃-rank SVD of Ω.
8:
Ũ1 , Ũ2 ← Ũ
. Split Ũ into two components.
9:
[Û, Σ̂, V̂] ← SVD(X0 , r)
. Truncated r-rank SVD of X0 .
−1
10:
Ã ← Û∗ X0 ṼΣ̃ Ũ∗1 Û.
. Low-rank approximation of A.
11:

12:
13:

−1

. Estimate reduced actuation matrix B̃.
. Estimate actuation matrix B̂.
. Eigendecomposition of Ã.

B̃ ← Û∗ X0 ṼΣ̃ Ũ∗2
−1
B̂ ← X0 ṼΣ̃ Ũ∗2
[W, Λ] ← EIG(Ã)
−1

. DMD modes of A.

Φ ← X0 ṼΣ̃ Ũ∗1 ÛW
end if
end procedure
Note: If λi = 0, then φi = Ũ1 Ũ∗1 Ûwi for step 13.
14:
15:
16:

2.1Dynamic mode decomposition with control

of states and inputs. The DMD modes can then be obtained following the procedure of the exact DMD algorithm 1 with X0 replaced with X0 − BΥ. For an unknown
B it is possible to compute the DMD modes and an approximation to the matrix B via regression [56]. For this
case, an augmented matrix containing both the state and
input snapshots is constructed
 
X
Ω=
(8)
Υ

The dynamic mode decomposition with control (DMDc)
method is a critically enabling extension of DMD [56].
DMDc disambiguates between the underlying dynamics
and the effects of actuation, modifying the basic assumption of DMD to include the effect of inputs uk ∈ Rq
xk+1 = Axk + Buk
where B ∈ Rn×q . The matrix form of the actuation is


Υ = u0

u1

...

um−1  .

(5)

along with an augmented matrix containing the two unknown system matrices


G= A B .
(9)

(6)

The regression problem is then formulated as
 

 X
0
X = A B
= GΩ
Υ

The system can be written in matrix form as:
X0 = AX + BΥ

(7)

(10)

where Ω ∈ R(n+q)×m is the combination of the state and
control snapshots. The DMDc algorithm with an unknown B is shown in Algorithm 2. Importantly, the SVD
now constructs a low-rank representation of the state and
input variables, both of which are used to produce approximations of the matrices A and B through low-rank
structures. Much like the DMD algorithm, DMDc relies on least-squares regression of the data to build a linear model A of the state dynamics that is disambiguated
from the discovered actuation matrix B.

where each column uk of the input snapshot matrix Υ is
the input at each snapshot in Eq. (6) and the data matrices
X, X0 are formulated in the same way as in Eq. (2).
A least-squares regression algorithm can once again be
used to determine the matrix A and its associated DMD
modes and eigenvalues. Two distinguishing cases can
be considered, when B is known and when B is unknown. When the input matrix B is known, or can be
well-estimated, the output is a simple linear combination
4

Algorithm 3 Compressive DMD [17]
Input: Measurements Y, Y0 , measurement matrix C, sparsifying basis Ψ, and target rank r.
Optional: X, X0 .
Output: cDMD spectrum Λ and modes Φ̂.
1: procedure C DMD(Y, Y 0 , C, r, [X, X0 ]).
2:
[UY , ΣY , VY ] ← SVD(Y, r)
. Truncated r-rank SVD of Y.
.
Low-rank approximation of AY .
3:
ÃY ← U∗Y Y0 VY Σ−1
Y
4:
[WY , ΛY ] ← EIG(ÃY )
. Eigendecomposition of ÃY .
5:
if X is known then
. Perform compressed DMD.
W
. Estimate DMD modes of A.
6:
Φ̂ ← X0 VY Σ−1
Y
Y
7:
else
. Perform compressed sensing DMD.
8:
ΦY ← Y0 VY Σ−1
. DMD modes of AY .
Y WY
9:
ΦS ← Compressed Sensing(ΦY , C, Ψ)
. Perform l1 minimization on φY i to solve for φS i .
10:
Φ̂ ← ΨΦS
. Estimate DMD modes of A.
11:
end if
12: end procedure
Note: If λi = 0, then φi = UwY,i , φY,i = UY wY,i for steps 6 and 8.

2.2

Compressed sensing and dynamic mode
decomposition

X, X0

Another innovation of the DMD algorithm addresses limitations on measurement and acquisition of a dynamical
system. Such limited data acquisition is often imposed by
physical constraints, such as data-transfer bandwidth in
particle image velocimetry (PIV), or the costs of sensors.
Given the low-rank nature of the spatiotemporal structures exhibited in many complex systems, we can utilize
ideas from compressed sensing [24, 19, 7] to reconstruct
the full high-dimensional state x from a small number
of measurements. Compressive DMD (cDMD) develops
a strategy for computing the dynamic mode decomposition from compressed or subsampled data [17].

C

Y, Y0

X0 = AX

Λ, Φ

Path 1:
Compressed DMD

ΛY , ΦY
Path 2: Compressed sensing DMD

Figure 1: Schematic of DMD and compressive DMD. Path
1 shows compressed DMD and path 2 shows compressed
sensing DMD ( [17]).

Consider compressed or subsampled data Y given by
Y = CX,

DMD

(11)

data. Importantly, when data and modes are sparse in
some transform basis, then there is an invariance of DMD
to measurement matrices that satisfy the restricted isometry property (RIP) from compressed sensing.

where C is a measurement matrix. There are two key
strategies to cDMD as illustrated in Fig. 1. First, it is possible to reconstruct full-state DMD modes from heavily
subsampled or compressed data using compressed sensing. This is called compressed sensing DMD, and it is appropriate to use when access to the full state space is not
possible due to constraints on physical measurements.
Second, if full-state snapshots are available, it is possible
to first compress the data, perform DMD, and then reconstruct by taking a linear combination of the snapshot
data, determined by the DMD on compressed data. This
is called compressed DMD. In this case, it is assumed that
one has access to the full high-dimensional state space
data. The theory for either of these methods relies on relationships between DMD on full-state and compressed

In addition to the the data matrices X and X0 , it is
also now required to specify a measurements matrix C.
These three matrices together are required to execute Algorithm 3. In practice, we often consider point measurements so that the rows of C are rows of the identity matrix. For DMD modes that are global in nature, point measurements naturally satisfy the RIP property as they are
incoherent with respect to the global modes. The success
of the method, both the compressive-sampling DMD and
the cDMD, has been demonstrated by Brunton et al. [17]
to be an effective strategy for subsampling of data and
the reconstruction of DMD modes and eigenvalues.
5

3

Compressive system identification

DMDc

X, X0 , Υ, [B]

A major benefit of dynamic mode decomposition is that
it provides a physically interpretable and highly extensible linear model framework, which enables the incorporation of actuation inputs and sparse output measurements. When combined, these innovations result in the
compressive DMD with control (cDMDc) architecture for
compressive system identification, where low-order models are identified from input–output measurements. In
contrast to traditional system identification, the reduced
states of the cDMDc models may be used to reconstruct
the high-dimensional state space via compressed sensing,
adding physical interpretability to the models. Thus, cDMDc relies on the existence of a few dominant coherent
patterns, which in turn facilitates sparse measurements.

=

Axk + Buk

(12a)

yk

=

Cxk .

(12b)

X = AX + BΥ
Path 1:
Compressed DMDc

C

Y, Y0 , Υ, [B]

ΛY , ΦY , [B̂Y ]

Path 2: Compressed sensing DMDc

Figure 2: Schematic of compressive DMD control. Path 1
shows compressed DMDc and path 2 shows compressed
sensing DMDc.

We now consider a general input–output system with
high-dimensional state x:
xk+1

Λ, Φ, [B̂]

0

Assumption 2. The columns of the full-state output matrix
X0 are in the subspace of the upper left singular vectors Ũ1 of
the full-state augmented matrix Ω so that Ũ1 Ũ∗1 X0 ≈ X0 .
Lemma 1. If Assumption 1 holds, we have an identity

As in DMD, the goal is to obtain the leading eigendecomposition of A, resulting in a low-order model in terms of
a few DMD modes. However, now we must account for
limited measurements in the output y and disambiguate
internal state dynamics from the effect of actuation u.
Writing Eq. (12) in terms of data matrices yields:

−1

ṼΣ̃

−1

= ṼY Σ̃Y Ũ∗Y,1 CŨ1 , similar to that derived in [17]:

Proof.

Y = CX
∗
ŨY,1 Σ̃Y ṼY
= CŨ1 Σ̃Ṽ∗
∗
ṼY
ṼΣ̃

−1

= Σ̃Y Ũ∗Y,1 CŨ1

(13b)

∗
ṼY ṼY
ṼΣ̃

−1

= ṼY Σ̃Y Ũ∗Y,1 CŨ1

Under certain conditions it is possible to apply DMDc
to the compressed data Y and then recover full-state
DMD modes via compressed sensing. As in the compressive DMD algorithm [17], if full-state data X is available, significant computational savings may be attained
by compressing the data and working in the compressed
subspace. If full-state data is unavailable, it may still be
possible to reconstruct full-state modes via convex optimization, exploiting sparsity of the modes in a generic
basis, such as Fourier or wavelets. The cDMDc framework is shown in Fig. 2, and is described in algorithm 4.

ṼΣ̃

−1

= ṼY Σ̃Y Ũ∗Y,1 CŨ1 .

0

=

AX + BΥ

(13a)

Y

=

CX.

X

−1

−1

−1

The next theorem uses the following definitions of A
and AY from the DMDc and cDMDc algorithms:
A = X0 ṼΣ̃

−1

Ũ∗1

−1

AY = Y0 ṼY Σ̃Y Ũ∗Y,1 .
Theorem 1. If Assumption 1 and 2 hold, then
CAX0 = AY CX0 .

(14)

We now prove that when compressed, DMD eigenvectors of the full data X become DMD eigenvectors of the
compressed data Y. This section relies on notation developed above, which is consolidated in the nomenclature.
Matrices with a subscript Y are computed on compressed
data, and matrices with a tilde are computed from the
SVD of the augmented matrix Ω.

Proof. Using Lemma 1, we have

Assumption 1. The measurement matrix C preserves the
∗
temporal information in Ω so that ṼY ṼY
Ṽ = Ṽ. This requires the columns of Ṽ to be in the column space of ṼY and
will only be approximately satisfied with measurement noise.

Thus, the compression matrix C commutes with the
dynamics in A, when applied to data X0 . We use this to
obtain the central result: compressed dynamic modes of
the full data are dynamic modes of the compressed data.

CAX0 = C(X0 ṼΣ̃

−1

Ũ∗1 )X0

−1

= Y0 (ṼY Σ̃Y Ũ∗Y,1 CŨ1 )Ũ∗1 X0
= AY CŨ1 Ũ∗1 X0
= AY CX0 .

6

Algorithm 4 Compressive DMD with control
Input: Measurement snapshot matrices Y, Y0 , sensing matrix C, input snapshot matrix Υ, and target rank r, r̃.
Optional: Full-state data matrices X, X0 and actuation matrix B.
Output: cDMDc spectrum Λ and modes Φ̂, [B̂].
1: procedure C DMD C (Y, Y 0 , C, r, r̃, [X, X0 , B])
2:
if X is known then
3:
if B is known then
4:
Λ, Φ̂ ← cDMD(Y, Y0 − CBΥ, C, X, X0 , r)
. Perform compressed DMD.
5:
else
6:
ΛY , ÛY , ŨY,1 , ŨY,2 , Σ̃Y , ṼY , WY
7:
← DMDc(Y, Y0 , Υ, r, r̃, X, X0 − BΥ)
. Perform DMDc.
−1 ∗
0
8:
Φ̂ ← X ṼY Σ̃Y ŨY,1 ÛY WY
. Estimate DMD modes of A.
−1

. Estimate actuation matrix B̂.

B̂ ← X0 ṼY Σ̃Y Ũ∗Y,2
end if
else
12:
if B is known then
13:
Λ, Φ̂ ← cDMD(Y, Y0 − CBΥ, C, r)
14:
else
15:
ΛY , ΦY , B̂Y ← DMDc(Y, Y0 , Υ, r, r̃)
16:
Φ̂ ← Compressed Sensing(ΦY , Ψ)
17:
B̂ ← Compressed Sensing(B̂Y , Ψ)
18:
end if
19:
end if
20: end procedure
Note: If λi = 0, then φi = ŨŨ∗Y,1 ÛY wY,i for step 7.
9:
10:
11:

. Perform compressed sensing DMD.
. Perform DMDc.
. Estimate DMD modes of A.
. Estimate actuation matrix B̂.

Theorem 2. If Assumptions 1 and 2 hold, then compressing a
full-state DMDc eigenvector φ yields a DMDc eigenvector of
the compressed data with the same eigenvalue:

Compressed recovery of the actuation matrix. We now
establish a similar relationship between the full-state actuation matrix B and the compressed matrix BY .

(15)

Assumption 3. The columns of ŨY,2 are spanned by those of
Ũ2 , so Ũ2 Ũ∗2 ŨY,2 = ŨY,2 and Ũ∗Y,2 Ũ2 Ũ∗2 = Ũ∗Y,2 .

AY Cφ = λCφ.
Proof.

−1

AY Cφ = AY CX0 ṼΣ̃

−1

= CAX0 ṼΣ̃

−1

Ũ∗1 Ûw

Lemma 2. We have an identity ṼΣ̃

Ũ∗1 Ûw

Proof. Expanding Υ in the SVD bases of ΩY and Ω, we
find:

= CAφ
= λCφ.

∗
ŨY,2 Σ̃Y ṼY
= Ũ2 Σ̃Ṽ∗

If C is chosen poorly so that φ is in its nullspace,
then Theorem 2 applies trivially. This theorem does not
guarantee that every eigenvector of AY is a compressed
eigenvector of A, although under reasonable assumptions the dominant eigenvalues of AY will approximate
those of A, as shown in [17]. We then have
φY = CφX = CΨφS ,

−1

= ṼY Σ̃Y Ũ∗Y,2 Ũ2 .

∗
ṼY
ṼΣ̃

−1

= Σ̃Y Ũ∗Y,2 Ũ2

−1

∗
ṼY ṼY
ṼΣ̃

−1

= ṼY Σ̃Y Ũ∗Y,2 Ũ2

−1

−1

−1

ṼΣ̃ = ṼY Σ̃Y Ũ∗Y,2 Ũ2 .
The next theorem uses the following definitions of B
and BY from the DMDc and cDMDc algorithms:

(16)

B = X0 ṼΣ̃

where φS is the sparse representation of φX in a universal basis Ψ, such as a Fourier or wavelet basis. Thus,
it is possible to recover φS , and hence φX , from compressed DMD modes φY , given enough incoherent measurements [17].

−1

Ũ∗2

−1

BY = Y0 ṼY Σ̃Y Ũ∗Y,2 .
Theorem 3. If Assumptions 1 and 3 hold, then
CB = BY .
7

(17)

Proof. Using Lemma 2, we have
−1

CB = CX0 ṼΣ̃

Therefore, if the full-state systems is controllable, i.e.
the controllability matrix C is full rank, and the compressed measurement matrix C is not in the null space of
Eq. (20), then the compressed system is also controllable.
The controllability is preserved after compression under
certain conditions. The observability property regarding
the compressed sensing framework was discussed in [71].

Ũ∗2

−1

= Y0 ṼY Σ̃Y Ũ∗Y,2 Ũ2 Ũ∗2
−1

= Y0 ṼY Σ̃Y Ũ∗Y,2
= BY .
Therefore, we can first compute the DMD modes ΦY
and the compressed actuation matrix BY and then reconstruct the full-state Φ and B through compressed sensing.

3.1

3.2

Similar to cDMD, there are two main paths presented in
compressive DMD with control, depending on the availability of the full-state data matrix X. Specifically, we refer to the first path as compressed DMDc (Path 1 in Figure 2) if X is known and the second path as compressed
sensing DMDc (Path 2 in Figure 2) if X is only partially
known due to the lack of full state measurements:

Relationship to system identification

The result from Theorem 1 above also carries over to general impulse response parameters in the following theorem. These impulse response parameters, also known as
Markov parameters, are used extensively in system identification, for example to construct Hankel matrices in the
eigensystem realization algorithm (ERA) [34].

1. If the full-state measurements X are available, compressed DMDc is advantageous as the expensive calculations are performed on the compressed data Y
and the full-state modes are obtained by linearly
combining the snapshots of X.

Theorem 4. If Assumptions 1 and 2 hold, then
(18)

CAB = AY CB.

2. Without access to full-state measurements X, compressed sensing DMDc extracts the inherent dynamics from sub-sampled/compressed data in the matrix Y, disambiguating the effect of actuation. Fullstate modes may then be recovered, under the standard conditions of compressed sensing.

Proof. Using Lemma 1, we have
−1

CAB = C(X0 ṼΣ̃

Ũ∗1 )B

−1

= Y0 (ṼY Σ̃Y Ũ∗Y,1 CŨ1 )Ũ∗1 B
−1

= AY CŨ1 Ũ∗1 X0 ṼΣ̃
= AY CX0 ṼΣ̃

−1

Ũ∗2

In addition, two approaches are considered in each
path considering the prior knowledge of B. Generally,
the cDMDc algorithm can be simplified as a corrected
cDMD algorithm if B is known. Otherwise, DMDc is utilized to extract the underlying dynamics from the compressed states Y, Y0 and then the modes Φ and actuation
matrix B are reconstructed by either projecting onto the
full states X0 or through compressed sensing.

Ũ∗2

= AY CB.
Corollary 1. Theorem 4 can be expanded to further steps in
dynamics, for k ∈ Z+ , such as:
CAk B = AY k CB = AY k BY .

(19)

This theorem and corollary establish a surprising result: under certain conditions the iterative dynamics in
the impulse response parameters commute with the measurement matrix. Impulse response parameters are a cornerstone of subspace system identification methods, such
as ERA, and the above results simplify the dynamics.

In Algorithm 4, these four scenarios are discussed.
When both X and B are known, the spectrum ΛY and
modes Φ are obtained by performing compressed DMD
on the pre-compressed data Y and the shifted matrix
correcting for the effect of control Y0 − CBΥ. When
X is known and B is unknown, DMDc is computed on
Y, Y0 and the dynamic modes Φ and actuation matrix
B are reconstructed as a linear combination of the fullstate data X0 . When X is only partially known and B
is known, compressed sensing DMD is performed on Y
and Y0 − CBΥ and the full-state modes Φ are reconstructed from the compressed ΦY . When both X and B
are unavailable, DMDc is computed on Y, Y0 and the dynamic modes Φ and actuation matrix B are reconstructed
using compressed sensing.

Corollary 1 also has implications for the controllability
of the projected and full-state systems. Given the controllability matrix


C = B AB · · · An−1 B
(20)

we have the following relationship:


C C = C B AB · · · An−1 B


= BY AY BY · · · An−1
Y BY = CY .

Computational considerations

(21a)
(21b)
8

(a)

(b)

0.5

xi
0.25

x̃2

0

t
-0.25

-0.5
-0.5

0

x̃

(c)

0.1

Pi,1

0

0.25

1

(d)

i
1

P̂i,1 0

-0.1
0.1 0

Pi,2

-0.25

256

512

768

-1
10

1024

256

512

768

1024

256

512

768

1024

P̂i,2 0

0
-0.1

-1
0

256

512

768

1024

0

i
i
Figure 3: Two-dimensional system with known dynamics: (a) phase portrait (color denotes the progression of time), (b)
Figure
3: Two-dimensional
system with known
dynamics:
(a) phase
portrait
(color denotes
of time),
(b)
x-t
diagram
of inflated high-dimensional
system,
X 2 R1024⇥301
(first
31 snapshots
shown),the
(c)progression
two orthogonal
modes
1024×301
x-t
diagram
of
inflated
high-dimensional
system,
X
∈
R
(first
31
snapshots
shown),
(c)
two
orthogonal
modes
of P and (d) DCT coefficients of the two modes of P.
of P and (d) DCT coefficients of the two modes of P.

4 Results
4 Results

The state matrix Ã and input matrix B̃ are designed to
Theastate
matrix
Ã and input
matrix B̃ are designed to
yield
stable,
controllable
system:
yield a stable, controllable system:
 1



x̃ 1 
 0.9 0.2  x̃11 
 0.1 
= 0.9 0.2 x̃2 + 0.1 uk .
(22)
x̃
x̃22 k+1 =
0.1 0.9 x̃ 2 k + 0.01 uk .
(22)
x̃ k+1 | −0.1{z 0.9} x̃ k | 0.01
{z }
{z
}
{z
}
|
|
Ã
B̃

In this section, we present two numerical experiments
In this
section,
we presentDMD
two with
numerical
experiments
that
illustrate
compressive
control.
In the first
that
illustrate
compressive
DMD
with
control.
In the first
example, we investigate a spatially high-dimensional
example,
we
investigate
a
spatially
high-dimensional
system that is lifted from a two-dimensional system with
system dynamics
that is lifted
from
a two-dimensional
known
and
random
control input.system
In the with
secknown
dynamics
and
random
control
input.
theflow
second example, we model the vorticity field of a In
fluid
ond
example,
we
model
the
vorticity
field
of
a
fluid
flow
downstream of a pitching plate at low Reynolds numdownstream
of a pitching
at lowthe
Reynolds
number.
In both examples,
we plate
investigate
effectiveness
ber.
In
both
examples,
we
investigate
the
effectiveness
of different random measurement matrices, compression
of different
random measurement
ratios,
and actuation
input vectors.matrices, compression
ratios, and actuation input vectors.

Ã

B̃

The
The dynamics
dynamics are
are excited
excited via
via Gaussian
Gaussian random
random input
input exexcitation
in
u
,
starting
from
an
initial
condition
k
in ⇤uTk , starting from an initial condition x
x00 =
=
⇥citation
0.25 0.25 T . The system is integrated from t = 0 to
0.25 0.25 . The system is integrated from t = 0 to
tt =
= 30
30 with
with aa time-step
time-step of
of ∆tt =
= 0.1,
0.1, resulting
resulting in
in 301
301 snapsnapshots.
A
sample
trajectory
of
the
low-dimensional
shots. A sample trajectory of the low-dimensional system
system
is
is shown
shown in
in Fig.
Fig. 33 (a).
(a).

To
we
To inflate
inflate the
the1 state
state dimension,
dimension,
we associate
associate each
each of
of
2
the
two
states
x̃
and
x̃
with
a
high-dimensional
the two states x̃1 and x̃2 with a high-dimensional mode
mode
that
that is
is sparse
sparse in
in the
the spatial
spatial wavenumber
wavenumber domain.
domain. These
These
modes,
shown
in
Fig.
3
(c),
are
given
by
the
modes,1024⇥2
shown in Fig. 3 (c), are given by the columns
columns of
of
P
P2
∈R
R1024×2,, where
where each
each column
column is
is constructed
constructed to
to have
have
K
K =
= 44 non-zero
non-zero elements
elements in
in the
the discrete
discrete cosine
cosine transtransform
(DCT)
basis.
The
nonzero
DCT
coefficients
form (DCT) basis. The nonzero DCT coefficients of
of each
each
mode
mode have
have the
the same
same wave
wave number
number with
with different
different magnimagnitudes,
tudes, as
as shown
shown in
in Fig.
Fig. 33 (d).
(d). Thus,
Thus, it
it is
is possible
possible to
to use
use

4.1
4.1 Stochastically
Stochastically forced
forced linear
linear system
system
This
This experiment
experiment is
is designed
designed to
to test
test the
the compressive
compressive DMD
DMD
with
control
framework
on
an
example
with control framework on an example where
where the
the lowlowrank
rank dynamics
dynamics are
are known.
known. Thus,
Thus, it
it is
is possible
possible to
to directly
directly
compare
compare the
the true
true eigenvalue
eigenvalue spectrum
spectrum and
and spatiotempospatiotemporal
modes
with
those
obtained
via
compressed
ral modes with those obtained via compressed DMDc
DMDc
and
and compressed
compressed sensing
sensing DMDc.
DMDc.
99

B known
B known

C-type 1
C-type 1
C-type 2
C-type 2
C-type 3
C-type 3
C-type 1
C-type 1

BB
unknown
unknown C-type 2
C-type 2
C-type
C-type33

DMDc
c-DMDc
cs-DMDc
DMDc
c-DMDc
cs-DMDc
3.631e-13
3.627e-13
3.627e-13
3.631e-13
3.627e-13
3.627e-13
3.631 e-13
3.443 e-13
3.443 e-13
3.631 e-13
3.443 e-13
3.443 e-13
3.631e-13
4.210e-13
4.210e-13
3.631e-13
4.210e-13
4.210e-13
1.758e-16
4.840e-17
9.558e-17
1.758e-16 5.396e-13
4.840e-17 5.926e-13
9.558e-17
4.481e-13
4.481e-13
5.396e-13
5.926e-13
1.758e-16
4.731e-17
2.845e-16
1.758e-16
4.731e-17
2.845e-16
4.481e-13
4.159e-13
3.872e-13
4.481e-13
4.159e-13
3.872e-13
1.758e-16
3.022e-16
2.785e-16
1.758e-16
3.022e-16
2.785e-16
4.481e-13
4.322e-13
5.691e-13
4.481e-13
4.322e-13
5.691e-13

ˆ kF and
Table
1: 1:
Normalized
kB B̂k using DMDc, compressed DMDc, and compressed sensing
Table
Normalizederror
errorofofk kΦ − Φ̂k
F and kB − B̂k22 using DMDc, compressed DMDc, and compressed sensing
DMDc.
Three
types
ofof
compression
(C-type1),
1),Gaussian
Gaussianrandom
randomprojections
projections
DMDc.
Three
types
compressionare
areshown:
shown:uniform
uniform random
random projections
projections (C-type
(C-type
2)
and
single
pixel
measurements
(C-type
3).
When
B
is
unknown,
the
error
of
the
estimated
B̂
is
given
(C-type 2) and single pixel measurements (C-type 3). When B unknown, the error of the estimated B̂ is given inin
thethe
ˆ
upper
triangle,
and
the
error
cases, B
B isisaalinear
linearcombination
combinationofofcolumns
columns
upper
triangle,
and
the
errorofof Φ̂isisgiven
givenininthe
thelower
lowertriangle.
triangle. In all cases,
ofof
P.P.

(a) 0.01

DCT
matrix
thesparsifying
sparsifyingbasis
basisfor
forcompressed
compressed
thethe
DCT
matrix
asas
the
sensing.
With
spatialmodes
modesininP,P,ititisispossible
possibleto
to lift
lift
sensing.
With
thethe
spatial
low-dimensional
statex̃kx̃ktotoa ahigh-dimensional
high-dimensionalstate
state
thethe
low-dimensional
state
. The lifted state is shown in Fig. 3 (b) for the
xk x=k =
Px̃Px̃
k . kThe lifted state is shown in Fig. 3 (b) for the
first
31
snapshots,
fromt t==0 0totot t==3.3. Note
Note that
that the
the
first 31 snapshots, from
high-dimensional
actuation
vector
B
is
chosen
to
be
in
high-dimensional actuation vector B is chosen to be in
the
span
of
the
columns
of
P,
i.e.
B
=
P
B̃,
so
that
the
the span of the columns of P, i.e. B = PB̃, so that the
actuation
only
exciteslow-dimensional
low-dimensionaldynamics;
dynamics; other
other
actuation
only
excites
types of actuation will be examined further in Fig. 6.
types of actuation will be examined further in Fig. 6.
To investigate the proposed cDMDc algorithm, we now
To investigate the proposed cDMDc algorithm, we now
consider compressed measurements y given by Eq. (12b).
consider compressed measurements y given by Eq. (12b).
Specifically, the compression matrix C can be built with
Specifically,
the compression matrix C can be built with
entries drawn from uniform (C-type 1), Gaussian (C-type
entries
drawn
from(C-type
uniform3)(C-type
1), Gaussian
2) or Bernoulli
distributions.
Fig. (C-type
4 shows
2) the
or Bernoulli
(C-type
3)
distributions.
Fig.
4 shows
DMDc mode reconstruction using compressed
DMDc
the(Path
DMDc
mode
reconstruction
using
compressed
1 in
Fig. 2)
and compressed
sensing
DMDc DMDc
(Path 2
(Path
1 in2)
Fig.
compressed
sensing DMDc (Path
2
in Fig.
for2)p and
= 128
Gaussian measurements.
In both
in cases,
Fig. 2)itfor
p = 128 that
Gaussian
measurements. actuation
In both
is assumed
the high-dimensional
cases,
it is
assumed
that the and
high-dimensional
actuation
input
vector
B is known,
the reconstructed
modes
input
vectorreproduce
B is known,
reconstructed
faithfully
the and
true the
coherent
structuresmodes
of the
faithfully
reproduce
the true
coherent
structures
ofcomthe
underlying
system (i.e.,
the DMD
modes
are a linear
underlying
system
(i.e., the of
DMD
are a linear
combination of
the columns
P). modes
The results
remain
unchanged
B is unknown
be reconbination
of when
the columns
of P). and
The must
resultsalso
remain
unstructed.
changed
when B is unknown and must also be reconstructed.
Table 1 shows the error in the reconstructed eigen-

B

0

-0.01
0.1 0

(b)

256

512

768

1 0

(c)

-0.1
0.1 0

256

512

768

256

512

768

1024

True
DMDc
C-DMDc
1024
CS-DMDc

2 0
-0.1
0

1024

i
Figure 4: Reconstruction of (a) actuation vector B, (b)-(c)
Figure 4: Reconstruction of (a) actuation vector B, (b)-(c)
modes from DMDc, compressed DMDc and compressed
modes from DMDc, compressed DMDc and compressed
sensing DMDc. For compressive DMD, 128 Gaussian ransensing DMDc. For compressive DMD, 128 Gaussian random
measurements are collected from 1024 state dimendom measurements are collected from 1024 state dimensions.
sions.

functions
Φ̂ andthe
estimated
actuation
vector B̂ (when
unTable
1 shows
error in
the reconstructed
eigenknown) ˆfor
DMDc,
compressed
DMDc,
and
functions
and
estimated
actuation
vector
B̂compressed
(when unsensing
against
theand
truecompressed
values. In
known)
forDMDc,
DMDc,compared
compressed
DMDc,
addition,
we investigate
effectthe
of true
different
sensing
sensing
DMDc,
compared the
against
values.
In
strategies
discussed
above.
In
all
cases,
the
reconstrucaddition, we investigate the effect of different sensing
tion error
is small, above.
as no noise
is added
simulated
strategies
discussed
In all
cases, to
thethe
reconstrucdata.
tion error is small, as no noise is added to the simulated
data.
10

10

Compressed sensing DMDc is able to uncover the
underlying
dynamics
andDMDc
spatio-temporal
fromthe
Compressed
sensing
is able tomodes
uncover
noiseless subsampled
relaxing the requirement
of
underlying
dynamics data,
and spatio-temporal
modes from
high-dimensional
measurements.
However,
in realis- of
noiseless
subsampled
data, relaxing
the requirement
tic experimental conditions,
measurement
noise in
willrealisalhigh-dimensional
measurements.
However,
ways
be present conditions,
and is known
to effect DMD
computic
experimental
measurement
noise
will altations.be Figure
5 shows
performance
DMDc
and
ways
present
and is the
known
to effectofDMD
compucompressed
DMDc
for varying
levels of measurement
tations.
Figure
5 shows
the performance
of DMDc and
noise, averaged
overfor
100varying
differentlevels
noise of
realizations
in
compressed
DMDc
measurement
each
case.
Note
that
compressed
sensing
DMDc
and
comnoise, averaged over 100 different noise realizations in
pressed
DMDc
identical eigenvalues,
as bothand
metheach
case.
Notehave
that compressed
sensing DMDc
comods compute the DMD spectrum from the same compressed DMDc have identical eigenvalues, as both methods compute the DMD spectrum from the same com-

(a)
(a)
⌘ ⇠ 10%
⌘ ⇠ 10%
(b)
(b)
⌘ ⇠ 25%
⌘ ⇠ 25%
(c)
(c)
⌘ ⇠ 50%
⌘ ⇠ 50%

B known
B known

B unknown
B unknown

0.2
0.2
0
0
-0.2
-0.2

0.2
0.2
True
True
DMDc
DMDc
cDMDc
cDMDc

0
0
-0.2
-0.2

0.1 0.3 0.5 0.7 0.9 1.1
0.2 0.1 0.3 0.5 0.7 0.9 1.1
0.2

0.2
0.2

0
0

0
0

-0.2
-0.2

-0.2
-0.2

0.1 0.3 0.5 0.7 0.9 1.1
0.2 0.1 0.3 0.5 0.7 0.9 1.1
0.2 C
C
0
0

-0.2
-0.2

0.1 0.3 0.5 0.7 0.9 1.1
0.1 0.3 0.5 0.7 0.9 1.1

0.9
0.9

0.1 0.3 0.5 0.7 0.9 1.1
0.2 0.1 0.3 0.5 0.7 0.9 1.1
0.2

0
0
-0.2
-0.2

0.1 0.3 0.5 0.7 0.9 1.1
0.1 0.3 0.5 0.7 0.9 1.1
0.1 0.3 0.5 0.7 0.9 1.1
0.1 0.3 0.5 0.7 0.9 1.1
Figure
5:
Noise
dependency
of
DMD
spectrum
for
the
true
system,
DMDc
andcDMDc
cDMDcbased
basedon
on100
100different
differentnoise
noise
Figure 5: Noise dependency of DMD spectrum for the true system, DMDc and
Figure 5: Noise
dependency
oftoDMD
spectrum
for
the ηtrue
system,
DMDc
and σcDMDc=based
on 100
different
noise
realizations.
Rows
correspond
different
noise
levels
∈
{0.1,
0.25,
0.5}
with
η
max(σ
),
where
σ
denotes
noise
i
i
realizations. Rows correspond to different noise levels ⌘ 2 {0.1, 0.25, 0.5} with noise = ⌘ max( i ), where i denotes
realizations.
Rows correspond
to different noise levels ⌘ 2 {0.1, 0.25, 0.5} with noise = ⌘ max( i ), where i denotes
the
standard
the
standarddeviation
deviation of
of each
each spatial
spatial measurement.
measurement.
the standard deviation of each spatial measurement.

BB22span{P
span{P11, ,PP22}}

B
Brandom
random

B
/2
B2
/ span{P
span{P11,,P
P22}}

B̂B̂

||

True
True
DMDc
DMDc
C-DMDc
C-DMDc
CS-DMDc
CS-DMDc

1,2 | |
1,2

Figure
6: Estimated
B̂B̂and
mode
magnitudes
B for
DMDc, compressed
DMDc
and
compressed
Figure
EstimatedB̂
andmode
modemagnitudes
magnitudesfor
fordifferent
differentchoices
choicesof
of
compressed
DMDc
and
Figure
6:6:
Estimated
and
for
different
choices
of B
B for
for DMDc,
DMDc,
compressed
DMDc
andcompressed
compressed
sensing
DMDc.
For
compressive
DMDc,
the
same
number
of
measurements
is
employed
as
for
Fig.
4.
sensing
DMDc.For
Forcompressive
compressiveDMDc,
DMDc,the
thesame
same number
number of
of measurements
measurements is
sensing
DMDc.
is employed
employedas
asfor
forFig.
Fig.4.4.
pressed
data.
For
moderate
noise
levels,
the
presseddata.
data.For
Formoderate
moderatenoise
noiselevels,
levels,the
thecompressed
compressed
pressed
compressed
DMDc
algorithm
provides
reasonably
accurate
eigenvalDMDc
algorithm
provides
reasonably
accurate
eigenvalDMDc
algorithm
provides
reasonably
accurate
eigenvalues,
although
they
are
less
accurate
than
those
predicted
ues, although they are less accurate than those predicted
ues,
although
theynoise
are less accurate
than those 10%,
predicted
by
byDMDc.
DMDc.As
Asthe
the noiseintensity
intensityisisincreased
increasedto
to 10%,both
both
by
DMDc.
As
the
noise
intensity
is
increased
to
greater
the
theDMDc
DMDcand
andcDMDc
cDMDceigenvalues
eigenvalues have
have significant
significant ererthan
50%, both
the DMDccomplex
and cDMDc eigenvalues
begin
rors,
rors,and
andininsome
somecases,
cases, complexconjugate
conjugate eigenvalues
eigenvalues
toare
deviate, and in some
cases, complex
conjugate eigenaremiscomputed
miscomputedas
aspurely
purelyreal
real eigenvalues.
eigenvalues. Although
Although
values
are of
miscomputed
as purelybyreal
eigenvalues.
Althe
effect
noise
is
exacerbated
compression,
the effect of noise is exacerbated by compression, sensisensithough
the
effect
of
noise
is
exacerbated
by
compression,
tivity
of
DMD
eigenvalues
with
measurement
noise
is
tivity of DMD eigenvalues with measurement noise is aa
sensitivity
of DMD
eigenvalues
with measurement
noise
known
and
been
studied
knownissue,
issue,
andhas
has
beenextensively
extensively
studiedand
andcharcharisacterized
a
known
issue,
and
has
been
extensively
studied
and
[4,
23,
32].
There
are
a
number
of
algorithacterized [4, 23, 32]. There are a number of algorithcharacterized
[4,
31,
23].
There
are
a
number
of
algorithmic
extensions
that
improve
the
eigenvalue
prediction
mic extensions that improve the eigenvalue prediction
mic
extensions
that improve
thetotal
eigenvalue
prediction
with
noise,
using
least
[32],
with
noise, including
including
using the
the
total
least squares
squares
[32],
aaforward-backward
symmetrizing
algorithm
[23],
with
noise,
including
using
the
total
least
squares
[31],
forward-backward symmetrizing algorithm [23], subsubDMD
DMD
or
aspace
forward-backward
symmetrizing
algorithm
[23],DMD
subspace
DMD[68],
[68],Bayesian
Bayesian
DMD[69],
[69],
oran
anoptimal
optimal
DMD
based
on
projection
methods
Each
these
space
DMD
[68], Bayesian
DMD
[67], or[2].
an optimal
based
onvariable
variable
projection
methods
[2].
Each of
ofDMD
these
methods
may
the
based
on variable
projectioncombined
methods with
[2].
of these
methods
maybe
beeffectively
effectively
combined
withEach
theproposed
proposed
cDMDc
to
more
methods
may
be effectively
combined
with eigenvalues.
the
proposed
cDMDcarchitecture
architecture
toyield
yield
moreaccurate
accurate
eigenvalues.
cDMDc architecture to yield more accurate eigenvalues.

4.1.1
of Actuation
4.1.1 Effect
Effect
4.1.1
Effect of
of Actuation
Actuation
Finally,
Finally,we
weinvestigate
investigatethe
theperformance
performance of
of cDMDc
cDMDc for
for difdifFinally,
we
investigate
the
performance
of
cDMDc
foraddifferent
choices
of
the
actuation
vector
B
in
Fig.
6.
ferent choices of the actuation vector B in Fig. 6. In
In
addition
to
the
cases
presented
so
far,
where
B
is
in
the
ferent
choices
of
the
actuation
vector
B
in
Fig.
6.
In
addition to the cases presented so far, where B is in the
subspace
P,
we
explore
when
randition
to of
the
presented
so
far, where
the
subspace
of
P,cases
we now
now
explore scenarios
scenarios
whenBB
Bisis
is in
randomly
or
in
complementary
subspace
of P, we
scenariossubspace
when B of
is
randomlygenerated
generated
ornow
inthe
theexplore
complementary
subspace
of P.
P.
For
and
of
domly
generated
or in theB,
complementary
subspace
of P.
For all
all types
types of
of actuation
actuation
B,
and regardless
regardless
of whether
whether
or
not
B
is
known,
compressed
DMDc
and
compressed
For
all Btypes
of actuation
B, and
regardless
of whether
or not
is known,
compressed
DMDc
and compressed
sensing
both
accurately
modes
or
not BDMDc
is known,
DMDc the
andtrue
compressed
sensing
DMDc
both compressed
accurately identify
identify
the
true
modes
and
,
which
is
also
consistent
with
DMDc.
When
B
sensing
both
accurately
identify
the true
modes
11 and DMDc
22 , which
is also
consistent
with DMDc.
When
B
is
unknown,
cDMDc
and
csDMDc
both
accurately
idenφ
and φ2 , which
is also
withaccurately
DMDc. When
is1unknown,
cDMDc
andconsistent
csDMDc both
iden-B
tify
B,
of
subspace
itit belongs
to,
as
is
unknown,
cDMDc
and
csDMDc
both accurately
tify
B,regardless
regardless
of the
the
subspace
belongs
to, as
as long
longidenas
the
columns
of
B
are
sparse.
The
algorithms
do
not
actify
B, regardless
subspace
belongs to,
long
the columns
of B of
arethe
sparse.
The italgorithms
doasnot
ac-as
curately
capture
the
B
matrix
itit is
(i.e.,
curately
capture
the
B sparse.
matrix when
when
is not
not sparse
sparse
(i.e.,
the
columns
of B
are
The algorithms
do not
acrandom
actuation
vector);
however,
in
this
case,
DMDc
random
actuation
vector);
however,
in
this
case,
DMDc
curately capture the B matrix when it is not sparse (i.e.,
also
actuation
vector.
alsomisidentifies
misidentifies
the
actuation
vector. in this case, DMDc
random
actuationthe
vector);
however,
also misidentifies the actuation vector.
11
11
11

Measurement window

(a)

t=0

(b)

t = 11.5

(c)

t = 21.5
q2

20
0

0
-20

Pitch velocity

(d)

Pitch angle [deg]

q1

5

-5
0

10

20
30
Convective Time, t

40

50

Figure 7: Simulation of a pitching airfoil: (a)-(c) Instantaneous vorticity fields at different times t = 0, 11.5, 21.5 and
Figure 7: Simulation of a pitching airfoil: (a)-(c) Instantaneous vorticity fields at different times t = 0, 11.5, 21.5 and
showing the measurement window, (d) time series of the actuation inputs specifying the pitching of the plate. Black
showing the measurement window, (d) time series of the actuation inputs specifying the pitching of the plate. Black
filled circles depict the time instants of the shown vorticity plots in (a)-(c).
filled circles depict the time instants of the shown vorticity plots in (a)-(c).

4.2 Fluid flow past a pitching plate

maneuver described
described in
in [55].
[56]. Six
Six rapid
rapid pitch
pitch maneumaneupitch maneuver
performed, and
and then
then the
the data
data is
is symmetrized
symmetrized by
by
vers are performed,
mirror image
image of
of these
these 251
251 snapshots
snapshotswith
with
concatenating a mirror
signed vorticity,
vorticity, in
in an
an attempt
attempt to
toidentify
identifyununthe opposite signed
Figure 77 illustrates
illustrates the
the vorticity
vorticityfield
fieldat
atdifdifbiased modes. Figure
= 0,
0, 11.5,
11.5, 21.5.
21.5. We
We focus
focus on
on aa downstream
downstream
ferent times, t =
window of
of the
the size
size 399
399 ×
⇥ 141
141 to
to mimic
mimic aa
measurement window
The actuation
actuation input
input takes
takes 4.6
4.6 convective
convective
PIV window. The
reach the
the observation
observation window,
window, and
and the
the acactime units to reach
that used
used as
as input
input for
for cDMDc
cDMDc and
and DMDc,
DMDc,
tuation input, that
accordingly. Similar
Similar small
small amplitude
amplitude pitching
pitching
is shifted accordingly.
been shown
shown to
to be
be well-approximated
well-approximated with
with
motions have been
[15].
linear models [18].

In
In the
the second example, we apply the proposed compressive
sive DMD with control algorithm to model the vorticity
field
field downstream of a pitching plate at Reynolds number
ber Re = 100. This pitching airfoil has been studied previously
viously in the context of reduced-order models for flow
control
31, 30].
67]. The flow is simulated uscontrol [15,
[66, 17,
22, 22,
18, 14,
ing
ing the immersed boundary projection method (IBPM)21
method
⇥ 159 on a
method [22,
[66, 67]
22] with a grid resolution of 799 ×
domain
domain of size 10 ⇥
× 2, nondimensionalized by the plate
length
length L. The flow is simulated with a time-step of ∆tt =
0.01
0.01 dimensionless convective time units, nondimensionalized
alized by the length L and free-stream velocity U . 251
snapshots
snapshots are sampled at a rate of ∆tt = 0.1. In this example,
ple, the airfoil is rapidly pitched up and down between
±5
±5◦ at irregular intervals in time, using the canonical

eigenvalues of
of Ã
Ã and
and Ã
ÃY
given by
by DMDc
DMDc and
and cDcDThe eigenvalues
Y given
shown in
in Fig.
Fig. 8.
8. The
The cloud
cloud of
of cDMDc
cDMDc eigeneigenMDc are shown
generated from
from an
an ensemble
ensemble of
of 50
50 realizations
realizations
values is generated
compressed measurements
measurements (i.e.,
(i.e., pp =
= 0.1n)
0.1n)
using 10% compressed
Gaussian random
random projections.
projections. We
We use
use rr =
= 99 in
in
based on Gaussian

1 Code available at https://github.com/cwrowley/ibpm
2 Code available at https://github.com/cwrowley/ibpm

12
12

(a)

0

0

-0.2

1b
2b
3b

0.2

0.4

0.6

0.8
Re( )

1

3a
2a

2

0
1b

2b
3b
4b

-4
1.2

1a

0

-2

4b

-0.4

4a

4

3a
2a
1a

0.2
Im( )

(b)

4a

DMDc
cDMDc

Im(!)

0.4

-0.04

-0.02

0
Re(!)

0.02

0.04

Figure 8: Comparison of the spectrum in (a) discrete time and (b) continuous time, both obtained using DMDc and
Figure 8: Comparison of the spectrum in (a) discrete time and (b) continuous time, both obtained using DMDc and
cDMDc with 10% Gaussian random measurements. The first r = 9 modes are shown, ordered by their real part.
cDMDc with 10% Gaussian random measurements. The first r = 9 modes are shown, ordered by their real part.

this study to identify the spatiotemporal coherent DMD
this study to identify the spatiotemporal coherent DMD
modes, and they capture 96% of the total energy, based
modes, and they capture 96% of the total energy, based
on the singular values. These nine modes are ordered
on the singular values. These nine modes are ordered
by their decay rate (i.e., the real part) and the higherby their decay rate (i.e., the real part) and the higherorder modes have larger frequencies of oscillation. The
order modes have larger frequencies of oscillation. The
discrete-time eigenvalues are located close to the unit cirdiscrete-time eigenvalues are located close to the unit circle, and they are slightly stable, as seen in the continuouscle, and they are slightly stable, as seen in the continuoustime plot. The DMDc eigenvalues appear to have a decay
time plot. The DMDc eigenvalues appear to have a decay
rate that is too small, although the dynamics of transients
rate that is too small, although the dynamics of transients
are captured quite well. Nevertheless, the cDMDc eigenare captured quite well. Nevertheless, the cDMDc eigenvalues agree well with the DMDc eigenvalues, which is
values agree well with the DMDc eigenvalues, which is
the goal.
the goal.
In Fig. 9, we present the spatial structure of the DMDc
In Fig. 9, we present the spatial structure of the DMDc
mode pairs and the B matrix. Similar DMD modes have
mode pairs and the B matrix. Similar DMD modes have
been observed in the flow past a cylinder [20]. As the tembeen observed in the flow past a cylinder [20]. As the temporal frequency associated with an eigenvalue increases,
poral frequency associated with an eigenvalue increases,
the modes are characterized by higher spatial wavenumthe modes are characterized by higher spatial wavenumbers, which is characteristic of bluff-body flows. The
bers, which is characteristic of bluff-body flows. The
accuracy of mode reconstruction is similar for Gaussian
accuracy of mode reconstruction is similar for Gaussian
random projections and single-pixel measurements, as
random projections and single-pixel measurements, as
indicated in Fig. 10. However, single pixel measurements
indicated in Fig. 10. However, single pixel measurements
may be less expensive and more realistic in real-world
may be less expensive and more realistic in real-world
applications, and sampling 10% of the original measureapplications, and sampling 10% of the original measurements results in a reconstruction accuracy of 90% in this
ments results in a reconstruction accuracy of 90% in this
example.
example.
Taking the DMDc results as the reference, Fig. 11 shows
Taking the DMDc results as the reference, Fig. 11 shows
the error of the eigenvalues with increasing number of
the error of the eigenvalues with increasing number of
measurements used in cDMDc. Note that the eigenvalmeasurements used in cDMDc. Note that the eigenvalues are the same for compressed DMDc and compressed
ues are the same for compressed DMDc and compressed
sensing DMDc, as shown in Algorithm 4 . The eigenvalue
sensing DMDc, as shown in Algorithm 4 . The eigenvalue
λ0 corresponding to zero frequency is estimated with the
0 corresponding to zero frequency is estimated with the

greatest accuracy from the fewest measurements, pregreatest accuracy from the fewest measurements, presumably because this mode contains the most energy in
sumably because this mode contains the most energy in
the flow. The error in eigenvalues associated with other
the flow. The error in eigenvalues associated with other
modes decreases logarithmically with increasing commodes decreases logarithmically with increasing compression ratio. Single pixel measurements have similar
pression ratio. Single pixel measurements have similar
performance compared with Gaussian random measureperformance compared with Gaussian random measurements for small compression ratios. When p = n, the
ments for small compression ratios. When p = n, the
error goes to zero for single pixel measurement, since the
error goes to zero for single pixel measurement, since the
measurement is an invertible permutation of the identity
measurement is an invertible permutation of the identity
matrix, and cDMDc is equivalent to DMDc in this case.
matrix, and cDMDc is equivalent to DMDc in this case.
Successful reconstruction using compressed sensing
Successful reconstruction using compressed sensing
relies on the sparsity of the state in some transform barelies on the sparsity of the state in some transform basis. Indeed, both Φ and B are sparse in the DCT basis.
sis. Indeed, both
and B are sparse in the DCT basis.
In general, the DMDc modes are more sparse than the
In general, the DMDc modes are more sparse than the
actuation matrix, and we choose K = 300 to ensure a
actuation matrix, and we choose K = 300 to ensure a
good reconstruction of the modes and B achieving an ergood reconstruction of the modes and B achieving an error of about 1%. In particular, we use the CoSaMP algoror of about 1%. In particular, we use the CoSaMP algorithm [54] to perform the l1 -minimization with 10 iterarithm [55] to perform the l1 -minimization with 10 iterations and the desired sparsity
of K = 300. The L2 ertions and the desired sparsity of K = 300. The L2 errors between the compressive DMDc and DMDc modes
rors between the compressive DMDc and DMDc modes
are shown in Fig. 10. The convergence of error versus
are shown in Fig. 10. The convergence of error versus
compression ratio in compressed sensing DMDc is much
compression ratio in compressed sensing DMDc is much
slower than that in compressed DMDc, especially for the
slower than that in compressed DMDc, especially for the
zero frequency mode.
zero frequency mode.
The necessary number of measurements is given from
The necessary number of measurements is given from
theory to be p ∼ 4K log10 (n/K) ≈ 2728, which cortheory to be p ⇠ 4K log10 (n/K) ⇡ 2728, which corresponds to a compression ratio of 5%. This matches
responds to a compression ratio of 5%. This matches
the observation that the error curves plateau at around
the observation that the error curves plateau at around
p/n = 0.05. Overall, the error for compressed measurep/n = 0.05. Overall, the error for compressed measurements using Gaussian random projections is comparable
ments using Gaussian random projections is comparable
with using single pixel measurements.
with using single pixel measurements.
13
13

DMDc

compressed DMDc

compressed sensing DMDc

Mode 0:

Mode 1:
Real

Imag
Mode 2:

Mode 3:

Mode 4:

B1 :

B2 :

Figure 9: Estimation of mode pairs (real and imaginary parts) and B using DMDc, compressed DMDc and compressed
Figure 9: Estimation of mode pairs (real and imaginary parts) and B using DMDc, compressed DMDc and compressed
sensing DMDc. The cDMDc modes and actuation matrix B (two columns, B and B ) are constructed using p Gaussian
sensing DMDc. The cDMDc modes and actuation matrix B (two columns, B11 and B22) are constructed using p Gaussian
random projections, where p/n = 0.1 (i.e., 10% compression ratio). Note that the results of single pixel measurements
random projections, where p/n = 0.1 (i.e., 10% compression ratio). Note that the results of single pixel measurements
are very similar (not shown).
are very similar (not shown).
14

14

(a)
(a) 100
0

(b)
(b) 100
0

10-3
10-3

10-3
10-3

i i 2

ˆ ˆk2k

csDMDc
csDMDc

i

cDMDc
cDMDc

10-10
10-10
0.001
0.001

10-5
10-5

k ki

10-5
10-5

10
10-1
10-1

k ki

i

i i

ˆ ˆk2k2

10
10-1
10-1

0.01
0.01

0.1
0.1

10-10
10-10
0.001
0.001

0.5 1
0.5 1

0.01

0.1

0.5 1

0.01
0.1
Compression ratio p/n
Compression
ratio
p/n 0.5 1
Compression ratio p/n
Compression ratio p/n
Figure 10: Error of estimated modes ˆˆ i and actuation matrix B̂ for increasing compression ratio using compressed
Figure
Error
estimatedmodes
modesφ̂ i and
and actuation
B̂
ratio
using
compressed
Figure
10:10:
Error
ofofestimated
matrix
B̂ for
for increasing
increasing
compression
ratio
using
compressed
DMDc
and
compressed
sensing
DMDci for
(a)actuation
Gaussianmatrix
random
projections
andcompression
(b) single pixel
measurements.
The
DMDc
and
compressed
sensing
DMDc
for
(a)
Gaussian
random
projections
and
(b)
single
pixel
measurements.
TheThe
DMDc
and
compressed
sensing
DMDc
(a)asGaussian
random
projections
(b) single pixel measurements.
DMDc
modes
and actuation
matrix
are for
used
the reference,
denoted
by andand
B, respectively.
DMDc
modes
and
actuationmatrix
matrixare
areused
usedas
as the
the reference,
reference, denoted
DMDc
modes
and
actuation
denotedby
by φand
andB,
B,respectively.
respectively.
(a) 1000
(b) 1000
(a) 10-1
(b) 10-1

10
10-1

10-3
10-3

10-3
10-3

10-10
10-10
0.001
0.001

i

10-5
10-5

| |i

i

10-5
10-5

| |i

i

ˆiˆ| i |

ˆiˆ| |

10
10-1

0.01
0.01

0.1
0.1

0.5 1
0.5 1

10-10
10-10
0.001
0.001

0.01
0.01

0.1
0.1

0.5 1
0.5 1

Compression ratio, p/n
Compression ratio, p/n
Compression ratio, p/n
Compression ratio, p/n
Figure 11: Error of eigenvalues for increasing compression ratio using (a) Gaussian random projections, and (b) single
Figure 11: Error of eigenvalues for increasing compression ratio using (a) Gaussian random projections, and (b) single
pixel measurements.
The DMDc
eigenvalues
are used as the reference
i . Gaussian random projections, and (b) single
Figure
Error of eigenvalues
foreigenvalues
increasing compression
using (a)
pixel11:
measurements.
The DMDc
are used as theratio
reference
i.
pixel measurements. The DMDc eigenvalues are used as the reference λi .

5

Conclusions
55 Conclusions

the proposed method. In addition, we have released our
the proposed method. In addition, we have released our
entire
code base todecompositions
promote reproducible
research with
and reaccurate
are achieved
entire
codemodal
base to promote reproducible
research and
re-surduce the barrier to implement these methods.
prisingly
few
measurements,
showcasing
the
efficacy
of
duce the barrier to implement these methods.
are a number
important
extensions
futureour
theThere
proposed
method.ofIn
addition,
we have and
released
There are a number of important extensions and future
directions
that
arise
out of this
work. First, it will be
inentire codethat
base
to promote
reproducible
directions
arise
out of this
work. First,research
it will beand
in- reteresting
to
further
investigate the
relationship
between
duce
the
barrier
to
implement
these
methods.
teresting to further investigate the relationship between
the controllable and observable subspaces and full-state
theThere
controllable
and observable
subspaces
and full-state
number
of important
extensions
future
recovery are
via acompressed
sensing. The
goal is aand
generrecovery
via
compressed
sensing.
The
goal
isitawill
generdirections
that
arise
out of this
work. of
First,
be inalized
theory
that
combines
the notion
controllability
alized
theory
that combines
the notion
of controllability
teresting
to further
investigate
the relationship
between
and
observability,
based
on the structure
of the A, B,
and
and
observability,and
based
on the structure
of theand
A, B,
and
thematrices,
controllable
observable
subspaces
full-state
C
and the notion
of sparse
signal recovery,
via
C
matrices,
and
the
notion
of
sparse
signal
recovery,
via
recovery
via of
compressed
sensing.
The P.
goal
is a also
generthe
structure
the low-rank
embedding
It may
the
structure
ofthat
the combines
low-rank embedding
It may also
alized
theory
notion P.
of
controllability
be
possible
to extend
results tothe
nonlinear
estimation
[65,
be possible to extend results to nonlinear estimation [65,
andand
observability,
onthrough
the structure
of the A,to B,
66]
control [36,based
42, 58]
the connection
66] and control [36, 42, 58] through the connection to
and Koopman
C matrices,
and theInnotion
of sparse
signal recovthe
operator.
addition,
the methods
dethe Koopman operator. In addition, the methods described
directlyofapplicable
to experimental
ery, viahere
the are
structure
the low-rank
embeddingmeaP. It
scribed here are directly applicable to experimental measurements
[64],
as theytodoextend
not require
access
to a model
may
also
be
possible
results
to
nonlinear
estisurements [64], as they do not require access to a model
of
the system.
Itand
maycontrol
be possible
to significantly
reduce
mation
[64,
65]
[57,
41,
35]
through
the
conof the system. It may be possible to significantly reduce

In Conclusions
summary, we have presented a unified framework
In summary, we have presented a unified framework
for compressive system identification based on the dyfor compressive system identification based on the dynamic mode decomposition. First, we describe the preIn namic
summary,
have presented
unified
framework
modewe
decomposition.
First,a we
describe
the previously developed compressed sensing DMD (csDMD)
forviously
compressive
system
identification
based
on(csDMD)
the dydeveloped
compressed
sensing
DMD
and DMD with control (DMDc) algorithms in a common
and DMD
control (DMDc)
algorithms
in a common
namic
mode with
decomposition.
First,
we describe
the premathematical framework, providing algorithmic implemathematical
framework, providing
algorithmic
impleviously
developed
sensing
(csDMD)
mentation
details.compressed
Next, we show
howDMD
it is possible
to
mentation
details.
Next,
we
show
how
it
is
possible
to
and
DMD with
control (DMDc)
algorithms
in a common
construct
reduced-order
models
from compressed
meaconstruct reduced-order
models
from
compressedimplemeamathematical
framework,
providing
algorithmic
surements and
then reconstruct
full-state
modes corresurements
and
then
reconstruct
full-state
modes
corrementation
westates
show via
howcompressed
it is possible
to
spondingdetails.
to the Next,
reduced
senssponding
to the reduced
statesfrom
via compressed
compressed measensconstruct
reduced-order
models
ing. This
lifting procedure
adds interpretability
to othing. Thisand
lifting
procedure
adds interpretability
to othsurements
then
reconstruct
modes
erwise black-box
models.
The full-state
compressed
DMDcorrewith
erwise
black-box
models.
The
compressed
DMD
with
sponding
to the reduced
states viaoncompressed
control algorithm
is demonstrated
two examplesenssyscontrol algorithm is demonstrated on two example sysing.
Thisincluding
lifting procedure
adds interpretability
to othtems,
a high-dimensional
discretized simulatems, including a high-dimensional discretized simulaerwise
black-box
models.
The compressed
with
tion of
fluid flow
past a pitching
airfoil. InDMD
both cases,
tion of fluid flow past a pitching airfoil. In both cases,
accurate
modal is
decompositions
with syssurcontrol
algorithm
demonstratedare
on achieved
two example
accurate modal decompositions are achieved with surprisingly
few
measurements,
showcasing
the
efficacy
of
tems,
including
a high-dimensional
discretized
simulaprisingly
few measurements,
showcasing
the efficacy
of

tion of fluid flow past a pitching airfoil. In both cases,

15
15

15

nection to the Koopman operator. In addition, the methods described here are directly applicable to experimental measurements [63], as they do not require access to a
model of the system. It may be possible to significantly
reduce the required spatial resolution in experiments, improving the effective bandwidth and enabling the characterization of faster flow phenomena. A sparsifying POD
basis may be obtained first with non-time-resolved measurements at full resolution. It is then possible to collect many fewer spatial measurements at much higher
temporal resolution, identify a reduced-order model, and
characterize the full-state modes in the offline library.
This also suggests that it may be possible to combine
space and time compressed sensing strategies for DMD.

[2] T. A SKHAM AND J. N. K UTZ, Variable projection methods for an optimized dynamic mode decomposition, arXiv
preprint arXiv:1704.02343, (2017).
[3] S. B AGHERI, Koopman-mode decomposition of the cylinder wake, Journal of Fluid Mechanics, 726 (2013),
pp. 596–623.
[4] S. B AGHERI, Effects of weak noise on oscillating flows:
linking quality factor, Floquet modes, and Koopman spectrum, Physics of Fluids, 26 (2014), p. 094104.
[5] Z. B AI , S. L. B RUNTON , B. W. B RUNTON , J. N.
K UTZ , E. K AISER , A. S POHN , AND B. R. N OACK,
Data-driven methods in fluid dynamics: Sparse classification from experimental data, in Invited chapter for
Whither Turbulence and Big Data in the 21st Century, 2015.

The growing intersection of dynamical systems, machine learning, and advanced optimization are driving
tremendous innovations in the characterization and control of complex systems [25, 43]. Although data is becoming increasingly abundant, there remain applications
such as feedback flow control, where real-time measurements are costly and control decisions must be made with
low latency to ensure robust performance. In these applications, techniques that strategically select sensor data
into the most relevant information will enable higher performance in more sophisticated flow control applications.
It is likely that the methods developed here may be combined with principled sensor selection methods to promote enhanced sparsity based on learned structures and
patterns [11, 5, 36, 47]. Moreover, it may also be possible
to enforce known physics or symmetries in the regression
procedure, as in [46], to improve model performance and
accelerate learning from data.

[6] Z. B AI , T. W IMALAJEEWA , Z. B ERGER , G. WANG ,
M. G LAUSER , AND P. K. VARSHNEY, Lowdimensional approach for reconstruction of airfoil
data via compressive sensing, AIAA Journal, 53 (2014),
pp. 920–933.
[7] R. G. B ARANIUK, Compressive sensing, IEEE Signal
Processing Magazine, 24 (2007), pp. 118–120.
[8] P. B ENNER , S. G UGERCIN , AND K. W ILLCOX, A survey of projection-based model reduction methods for parametric dynamical systems, SIAM Review, 57 (2015),
pp. 483–531.
[9] E. B ERGER , M. S ASTUBA , D. V OGT, B. J UNG , AND
H. B. A MOR, Estimation of perturbations in robotic behavior using dynamic mode decomposition, Journal of
Advanced Robotics, 29 (2015), pp. 331–343.

Acknowledgements

[10] D. A. B ISTRIAN AND L. M. N AVON, Randomized
dynamic mode decomposition for non-intrusive reduced
order modelling, International Journal for Numerical
Methods in Engineering, 112 (2017), pp. 3–25.

SLB acknowledges support from the Air Force Office of
Scientific Research (AFOSR FA9550-16-1-0650; Program
Manager: Dr. Douglas R. Smith) and the Army Research
Office (ARO W911NF-17-1-0118; Program Manager: Dr.
Matthew Munson). We would also like to thank several people for valuable discussions about sparse sensing, dynamic mode decomposition, and modeling highdimensional fluids, including Bing Brunton, Krithika
Manohar, Igor Mezić, Bernd Noack, and Sam Taira.

[11] B. W. B RUNTON , S. L. B RUNTON , J. L. P ROCTOR ,
AND J. N. K UTZ , Sparse sensor placement optimization
for classification, SIAM Journal on Applied Mathematics, 76 (2016), pp. 2099–2122.
[12] B. W. B RUNTON , L. A. J OHNSON , J. G. O JEMANN ,
AND J. N. K UTZ , Extracting spatial–temporal coherent
patterns in large-scale neural recordings using dynamic
mode decomposition, Journal of Neuroscience Methods, 258 (2016), pp. 1–15.

References
[1] H. A RBABI AND I. M EZI Ć, Ergodic theory, dynamic mode decomposition and computation of spectral properties of the koopman operator, arXiv preprint
arXiv:1611.06664, (2016).

[13] S. L. B RUNTON , B. W. B RUNTON , J. L. P ROCTOR ,
E. K AISER , AND J. N. K UTZ, Chaos as an intermittently forced linear system, Nature Communications, 8
(2017).
16

[14] S. L. B RUNTON , S. T. D AWSON , AND C. W. R OWLEY , State-space model identification and feedback control of unsteady aerodynamic forces, J. Fluid & Struc., 50
(2014), pp. 253–270.

[26] N. B. E RICHSON , S. L. B RUNTON , AND J. N.
K UTZ, Randomized dynamic mode decomposition,
arXiv:1702.02912, (2017).
[27] J. G ROSEK AND J. N. K UTZ, Dynamic mode decomposition for real-time background/foreground separation in
video, arXiv preprint arXiv:1404.7592, (2014).

[15] S. L. B RUNTON AND B. R. N OACK, Closed-loop turbulence control: Progress and challenges, Applied Mechanics Reviews, 67 (2015), pp. 050801–1–050801–48.

[28] F. G UENIAT, L. M ATHELIN , AND L. PASTUR, A dynamic mode decomposition approach for large and arbitrarily sampled systems, Physics of Fluids, 27 (2015),
p. 025113.

[16] S. L. B RUNTON , J. L. P ROCTOR , AND J. N. K UTZ,
Discovering governing equations from data by sparse
identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences, 113
(2016), pp. 3932–3937.

[29] R. H ECKEL AND H. B ÖLCSKEI, Identification of sparse
linear operators, arXiv preprint arXiv:1209.5187v2,
(2013).

[17] S. L. B RUNTON , J. L. P ROCTOR , J. H. T U , AND J. N.
K UTZ, Compressed sensing and dynamic mode decomposition, Journal of Computational Dynamics, 2 (2015),
pp. 165–191.

[30] M. S. H EMATI , S. T. D AWSON , AND C. W. R OWLEY , Parameter-varying aerodynamics models for aggressive pitching-response prediction, AIAA Journal,
(2016), pp. 1–9.

[18] S. L. B RUNTON , C. W. R OWLEY, AND D. R.
W ILLIAMS, Reduced-order unsteady aerodynamic models at low Reynolds numbers, Journal of Fluid Mechanics, 724 (2013), pp. 203–233.

[31] M. S. H EMATI , C. W. R OWLEY, E. A. D EEM , AND
L. N. C ATTAFESTA, De-biasing the dynamic mode
decomposition for applied Koopman spectral analysis,
arXiv:1502.03854, (2015).

[19] E. J. C ANDÈS , J. R OMBERG , AND T. TAO, Robust
uncertainty principles: exact signal reconstruction from
highly incomplete frequency information, IEEE Transactions on Information Theory, 52 (2006), pp. 489–509.

[32] P. J. H OLMES , J. L. L UMLEY, G. B ERKOOZ , AND
C. W. R OWLEY, Turbulence, coherent structures, dynamical systems and symmetry, Cambridge Monographs in Mechanics, Cambridge University Press,
Cambridge, England, 2nd ed., 2012.

[20] K. K. C HEN , J. H. T U , AND C. W. R OWLEY, Variants of dynamic mode decomposition: Boundary condition, Koopman, and Fourier analyses, Journal of Nonlinear Science, 22 (2012), pp. 887–915.

[33] M. R. J OVANOVI Ć , P. J. S CHMID , AND J. W.
N ICHOLS, Sparsity-promoting dynamic mode decomposition, Physics of Fluids, 26 (2014), p. 024103.

[21] Y. C HEN , Y. G U , AND A. O. H. III, Sparse lms for
system identification, IEEE International Conference
on Acoustics, Speech and Signal Processing (ICASSP
2009), (2009), p. 10700918.

[34] J. N. J UANG AND R. S. PAPPA, An eigensystem realization algorithm for modal parameter identification and
model reduction, Journal of Guidance, Control, and
Dynamics, 8 (1985), pp. 620–627.

[22] T. C OLONIUS AND K. TAIRA, A fast immersed boundary method using a nullspace approach and multi-domain
far-field boundary conditions, Computer Methods in
Applied Mechanics and Engineering, 197 (2008),
pp. 2131–2146.

[35] E. K AISER , J. N. K UTZ , AND S. L. B RUNTON, Datadriven discovery of Koopman eigenfunctions for control,
arXiv preprint arXiv:1707.01146, (2017).

[23] S. T. D AWSON , M. S. H EMATI , M. O. W ILLIAMS ,
AND C. W. R OWLEY , Characterizing and correcting for
the effect of sensor noise in the dynamic mode decomposition, Experiments in Fluids, 57 (2016), pp. 1–19.

[36] E. K AISER , M. M ORZYNSKI , G. D AVILLER , J. N.
K UTZ , B. W. B RUNTON , AND S. L. B RUNTON, Sparsity enabled cluster reduced-order models for control,
arXiv preprint arXiv:1701.00038, (2016).

[24] D. L. D ONOHO, Compressed sensing, IEEE Transactions on Information Theory, 52 (2006), pp. 1289–
1306.

[37] E. K AISER , B. R. N OACK , L. C ORDIER , A. S POHN ,
M. S EGOND , M. A BEL , G. D AVILLER , J. O STH ,
S. K RAJNOVIC , AND R. K. N IVEN, Cluster-based
reduced-order modelling of a mixing layer, J. Fluid
Mech., 754 (2014), pp. 365–414.

[25] T. D URIEZ , S. L. B RUNTON , AND B. R. N OACK, Machine Learning Control: Taming Nonlinear Dynamics
and Turbulence, Springer, 2016.
17

[38] N. K ALOUPTSIDIS , G. M ILEOUNIS , B. B ABADI , AND
V. TAROKH, Adaptive algorithms for sparse system identification, Signal Processing, (2011), pp. 1910–1919.

[51] I. M EZI Ć, Analysis of fluid flows via spectral properties
of the Koopman operator, Ann. Rev. Fluid Mech., 45
(2013), pp. 357–378.

[39] B. O. K OOPMAN, Hamiltonian systems and transformation in Hilbert space, Proceedings of the National
Academy of Sciences, 17 (1931), pp. 315–318.

[52] B. C. M OORE, Principal component analysis in linear
systems: Controllability, observability, and model reduction, IEEE Transactions on Automatic Control, AC26 (1981), pp. 17–32.

[40] Y. K OPSINIS , K. S LAVAKIS , AND S. T HEODORIDIS,
Online sparse system identification and signal reconstruction using projections onto weighted l1 balls, arXiv
preprint arXiv:1004.3040, (2010).

[53] M. N. AND D. C., Nonlinear system identification using compressed sensing, IEEE 2012 Conference Record
of the Forty Sixth Asilomar Conference on Signals, Systems and Computers (ASILOMAR), (2012),
p. 13415721.

[41] M. K ORDA AND I. M EZI Ć, Linear predictors for nonlinear dynamical systems: Koopman operator meets
model predictive control, arXiv:1611.03537, (2016).

[54] D. N EEDELL AND J. A. T ROPP, CoSaMP: iterative
signal recovery from incomplete and inaccurate samples,
Communications of the ACM, 53 (2010), pp. 93–100.

[42] B. K RAMER , P. G ROVER , P. B OUFOUNOS ,
M. B ENOSMAN , AND S. N ABI, Sparse sensing and
dmd based identification of flow regimes and bifurcations
in complex flows, arXiv preprint arXiv:1510.02831,
(2015).

[55] M. V. OL, A. A LTMAN , J. D. E LDREDGE , D. J. G AR MANN , AND Y. L IAN, Résumé of the AIAA FDTC
low Reynolds number discussion group’s canonical cases.
AIAA Paper 2010-1085, 48th Aerospace Sciences
Meeting, January 2010.

[43] J. N. K UTZ , S. L. B RUNTON , B. W. B RUNTON ,
AND J. L. P ROCTOR , Dynamic Mode Decomposition:
Data-Driven Modeling of Complex Systems, Society for
Industrial and Applied Mathematics, Philadelphia,
PA, 2016.

[56] J. L. P ROCTOR , S. L. B RUNTON , AND J. N. K UTZ,
Dynamic mode decomposition with control, SIAM J. Applied Dynamical Systems, 15 (2016), pp. 142–161.

[44] H. L EUNG, System identification using chaos with application to equalization of a chaotic modulation system,
IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 45 (1998), pp. 314–
320.

[57] J. L. P ROCTOR , S. L. B RUNTON , AND J. N. K UTZ,
Generalizing Koopman theory to allow for inputs and
control, arXiv:1602.07647, (2016).
[58] J. L. P ROCTOR AND P. A. E CKHOFF, Discovering dynamic patterns from infectious disease data using dynamic mode decomposition, International health, 7
(2015), pp. 139–145.

[45] L. L JUNG, System identification, Wiley Online Library,
1999.
[46] J.-C. L OISEAU AND S. L. B RUNTON, Constrained sparse Galerkin regression, arXiv preprint
arXiv:1611.03271, (2016).

[59] C. W. R OWLEY, Model reduction for fluids using balanced proper orthogonal decomposition., International
Journal of Bifurcation and Chaos, 15 (2005), pp. 997–
1013.

[47] K. M., B. W. B RUNTON , J. N. K UTZ , AND S. L.
B RUNTON, Data-driven sparse sensor placement, To appear in IEEE Control Systems Magazine, (2017).

[60] C. W. R OWLEY, I. M EZI Ć , S. B AGHERI , P. S CHLATTER , AND D. S. H ENNINGSON, Spectral analysis
of nonlinear flows, Journal of Fluid Mechanics, 641
(2009), pp. 115–127.

[48] Z. M A , S. A HUJA , AND C. W. R OWLEY, Reduced order models for control of fluids using the eigensystem realization algorithm, Theor. and Comp. Fluid Dyn., 25
(2011), pp. 233–247.

[61] P. J. S CHMID, Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Mechanics,
656 (2010), pp. 5–28.

[49] J. M ANN AND J. N. K UTZ, Dynamic mode decomposition for financial trading strategies, Quantitative Finance, (2016), pp. 1–13.

[62] P. J. S CHMID , K. E. M EYER , AND O. P UST, Dynamic
mode decomposition and proper orthogonal decomposition of flow in a lid-driven cylindrical cavity. PIV090186, 8th International Symposium on Particle Image Velocimetry, 2009.

[50] I. M EZI Ć, Spectral properties of dynamical systems,
model reduction and decompositions, Nonlinear Dynamics, 41 (2005), pp. 309–325.
18

[63] A. S. S HARMA , I. M EZI Ć , AND B. J. M C K EON, Correspondence between koopman mode decomposition, resolvent mode decomposition, and invariant solutions of
the navier-stokes equations, Physical Review Fluids, 1
(2016), p. 032402.
[64] A. S URANA, Koopman operator based observer synthesis for control-affine nonlinear systems, in 55th IEEE
Conference on Decision and Control (CDC, 2016,
pp. 6492–6499.
[65] A. S URANA AND A. B ANASZUK, Linear observer synthesis for nonlinear systems using koopman operator
framework, IFAC-PapersOnLine, 49 (2016), pp. 716–
723.
[66] K. TAIRA AND T. C OLONIUS, The immersed boundary
method: a projection approach., Journal of Computational Physics, 225 (2007), pp. 2118–2137.
[67] N. TAKEISHI , Y. K AWAHARA , Y. TABEI , AND
T. YAIRI, Bayesian dynamic mode decomposition,
Twenty-Sixth International Joint Conference on Artificial Intelligence, (2017).
[68] N. TAKEISHI , Y. K AWAHARA , AND T. YAIRI, Subspace dynamic mode decomposition for stochastic Koopman analysis, arXiv:1705.04908, (2017).
[69] J. H. T U , C. W. R OWLEY, J. N. K UTZ , AND J. K.
S HANG, Spectral analysis of fluid flows using subNyquist-rate PIV data, Experiments in Fluids, 55
(2014), pp. 1–13.
[70] J. H. T U , C. W. R OWLEY, D. M. L UCHTENBURG ,
S. L. B RUNTON , AND J. N. K UTZ, On dynamic
mode decomposition: Theory and applications, Journal
of Computational Dynamics, 1 (2014), pp. 391–421.
[71] M. B. WAKIN , B. M. S ANANDAJI , AND T. L. V IN CENT , On the observability of linear systems from random, compressive measurements, in 49th IEEE Conference on Decision and Control (CDC), Dec 2010,
pp. 4447–4454.
[72] K. W ILLCOX AND J. P ERAIRE, Balanced model reduction via the proper orthogonal decomposition, AIAA
Journal, 40 (2002), pp. 2323–2330.
[73] J. Y. G U AND S. M EI, l0 norm constraint lms algorithm
for sparse system identification, IEEE Signal Processing
Letters, (2009), pp. 774–777.

19

