# Dynamic mode decomposition
for compressive system identification

Zhe Bai1∗, Eurika Kaiser1, Joshua L. Proctor3, J. Nathan Kutz2, Steven L.
Brunton1

1 Department of Mechanical Engineering, University of Washington, Seattle, WA
98195, United States
2 Department of Applied Mathematics, University of Washington, Seattle, WA
98195, United States
3 Institute of Disease Modeling, Bellevue, WA 98004, United States

###### 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 111Code is publicly available at:
https://github.com/zhbai/cDMDc. 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

## Nomenclature

$\mathbf{A},\mathbf{\tilde{A}}$ State transition matrix for $\mathbf{x}$ and
$\tilde{\mathbf{x}}$ $\mathbf{A_{Y}},\mathbf{\tilde{A}_{Y}}$ State transition
matrix for $\mathbf{y}$ $\mathbf{B},\mathbf{\tilde{B}},\mathbf{\hat{B}}$
Actuation matrix $\hat{\mathbf{B}}$ Estimate of $\mathbf{B}$ $\mathbf{b}$
Vector of DMD amplitudes $\mathbf{C}$ Output measurement matrix $\mathbf{F}$
Discrete cosine transform (DCT) $\mathbf{G}$ Augmented matrix of $\mathbf{A}$
and $\mathbf{B}$ $K$ Number of nonzero DCT coeffs. $\mathbf{P}$ Projection
matrix acting on $\mathbf{x}$ $\mathbf{q}$ Spatial coordinate $r,\tilde{r}$
Rank of truncated SVD of $\mathbf{X},\boldsymbol{\Omega}$ $\mathbf{s}$ DCT
coefficients of $\mathbf{x}$ $t$ Time $\Delta t$ Time step $t_{k}$ $k$th
discrete time step $\mathbf{U},\mathbf{\hat{U}},\mathbf{\tilde{U}}$ Left
singular vectors (POD modes) of
$\mathbf{X},\mathbf{X}^{\prime},\boldsymbol{\Omega}$
$\mathbf{U_{Y}},\mathbf{\hat{U}_{Y}},\mathbf{\tilde{U}_{Y}}$ Left singular
vectors (POD modes) of
$\mathbf{Y},\mathbf{Y}^{\prime},\boldsymbol{\Omega_{\mathbf{Y}}}$ $\mathbf{u}$
Control input, $\mathbf{u}\in\mathbb{R}^{q}$
$\mathbf{V},\mathbf{\hat{V}},\mathbf{\tilde{V}}$ Right singular vectors of
$\mathbf{X},\mathbf{X}^{\prime},\boldsymbol{\Omega}$
$\mathbf{V_{Y}},\mathbf{\hat{V}_{Y}},\mathbf{\tilde{V}_{Y}}$ Right singular
vectors of $\mathbf{Y},\mathbf{Y}^{\prime},\boldsymbol{\Omega_{\mathbf{Y}}}$
$\mathbf{W},\mathbf{W_{Y}}$ Eigenvectors of
$\mathbf{\tilde{A}},\mathbf{\tilde{A}_{Y}}$ $\mathbf{X}$ Data matrix of
states, $\mathbf{X}\in\mathbb{R}^{n\times m}$ $\mathbf{X}^{\prime}$ Shifted
data matrix, $\mathbf{X}^{\prime}\in\mathbb{R}^{n\times m}$
$\tilde{\mathbf{X}}$ Data matrix of reduced states,
$\tilde{\mathbf{X}}\in\mathbb{R}^{r\times m}$ $\xi$ Spatial variable
$\mathbf{x}$ State vector, $\mathbf{x}\in\mathbb{R}^{n}$ $\tilde{\mathbf{x}}$
Low-rank state vector, $\tilde{\mathbf{x}}\in\mathbb{R}^{r}$ $\mathbf{Y}$ Data
matrix of measurements, $\mathbf{Y}\in\mathbb{R}^{p\times m}$
$\mathbf{Y}^{\prime}$ Shifted matrix of measurements,
$\mathbf{Y}^{\prime}\in\mathbb{R}^{p\times m}$ $\mathbf{y}$ Output vector,
$\mathbf{y}\in\mathbb{R}^{p}$ $\eta$ Noise magnitude ${\lambda}$ DMD
eigenvalue of $\mathbf{A}$ $\nu$ Wavenumber
$\boldsymbol{\Lambda},\boldsymbol{\Lambda}_{\mathbf{Y}}$ Matrix of DMD
eigenvalues of $\mathbf{A},\mathbf{A_{Y}}$ $\boldsymbol{\Omega}$ State and
input snapshot matrices
$[\mathbf{X}^{T}\enspace\boldsymbol{\Upsilon}^{T}]^{T}$
$\boldsymbol{\Omega_{\mathbf{Y}}}$ Output and input snapshot matrices
$[\mathbf{Y}^{T}\enspace\boldsymbol{\Upsilon}^{T}]^{T}$ $\omega$ Continuous-
time DMD eigenvalue of $\mathbf{A}$
$\boldsymbol{\Phi},\boldsymbol{\Phi_{\mathbf{Y}}}$ Matrix of DMD (DMDc) modes
of $\mathbf{A},\mathbf{A_{Y}}$ $\boldsymbol{\Phi_{\mathbf{S}}}$ Matrix of
sparse representation of $\boldsymbol{\Phi}$ $\hat{\boldsymbol{\Phi}}$
Compressed sensing estimate of $\boldsymbol{\Phi}$ $\boldsymbol{\Psi}$
Orthonormal basis (e.g. Fourier or POD)
$\boldsymbol{\Sigma},\boldsymbol{\hat{\Sigma}},\boldsymbol{\tilde{\Sigma}}$
Matrix of singular values of
$\mathbf{X},\mathbf{X}^{\prime},\boldsymbol{\Omega}$
$\boldsymbol{\Sigma_{\mathbf{Y}}},\boldsymbol{\hat{\Sigma}_{\mathbf{Y}}},\boldsymbol{\tilde{\Sigma}_{\mathbf{Y}}}$
Matrix of singular values of
$\mathbf{Y},\mathbf{Y}^{\prime},\boldsymbol{\Omega_{\mathbf{Y}}}$
$\boldsymbol{\Theta}$ Product of measurement matrix and sparsifying basis,
$\boldsymbol{\Theta}=\mathbf{C}\boldsymbol{\Psi}$ $\boldsymbol{\Upsilon}$ Data
matrix of control inputs

##  1 Introduction

Dynamic mode decomposition (DMD) is a dimensionality reduction technique
introduced by Schmid [61] in the fluid dynamics community to decompose high-
dimensional 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, oscillate 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].

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

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 Jovanović 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 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 reduced-order models on downsampled
spatial measurements of a high-dimensional system, and then reconstruct the
full-state 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.

Algorithm 1 Exact DMD [70]

1:Data matrix $\mathbf{X}$, shifted data matrix $\mathbf{X}^{\prime}$, and
target rank $r$.

2:DMD spectrum $\boldsymbol{\Lambda}$ and modes $\boldsymbol{\Phi}$.

3:procedure DMD($\mathbf{X},\mathbf{X}^{\prime},r$)

4: $[\mathbf{U},\boldsymbol{\Sigma},\mathbf{V}]\leftarrow$ SVD($\mathbf{X},r$)
$\triangleright$ Truncated $r$-rank SVD of $\mathbf{X}$.

5:
$\mathbf{\tilde{A}}\leftarrow\mathbf{U}^{*}\mathbf{X}^{\prime}\mathbf{V}\boldsymbol{\Sigma}^{-1}$
$\triangleright$ Low-rank approximation of $\mathbf{A}$.

6: $[\mathbf{W},\boldsymbol{\Lambda}]\leftarrow\text{EIG}(\mathbf{\tilde{A}})$
$\triangleright$ Eigen-decomposition of $\mathbf{\tilde{A}}$.

7:
$\boldsymbol{\Phi}\leftarrow\mathbf{X}^{\prime}\mathbf{V}\boldsymbol{\Sigma}^{-1}\mathbf{W}$
$\triangleright$ DMD modes of $\mathbf{A}$.

8:end procedure

Note: If ${\lambda}_{i}=0$, then
$\boldsymbol{\phi}_{i}=\mathbf{U}\mathbf{w}_{i}$ for step $5$. In the original
DMD algorithm [62] all modes are computed as
$\boldsymbol{\phi}_{i}=\mathbf{U}\mathbf{w}_{i}$.

##  2 Background: The Dynamic Mode Decomposition

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
${\bf X}$ into the rank $r$ approximation

$\displaystyle\mathbf{X}\approx\boldsymbol{\Phi}\boldsymbol{\Lambda}\mathbf{b}$
(1)

where the columns of $\boldsymbol{\Phi}$ are the DMD modes (spatial
structures), the elements of the diagonal matrix $\boldsymbol{\Lambda}$ are
the corresponding DMD eigenvalues (with angular frequency
$\lambda_{i}=\exp(\omega_{i}\Delta t)$), and the vector $\mathbf{b}$
determines the weighting of each of the $r$ modes. Thus, each spatial mode in
$\boldsymbol{\Phi}$ is associated with a single temporal eigenvalue of
$\boldsymbol{\Lambda}$.

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 future-state prediction [33, 2] and control for input-output
systems [56], even when using only a small number of measurements [17].

The exact DMD algorithm formulates the decomposition as a least-squares
regression. Specifically, a series of snapshots
$\mathbf{x}_{k}\in\mathbb{R}^{n}$ sampled at discrete instances in time
$t_{k}$, $k=0,\ldots,m$ are arranged into the data matrix

$\displaystyle\mathbf{X}$
$\displaystyle=\begin{bmatrix}\vline&\vline&&\vline\\\
\mathbf{x}_{0}&\mathbf{x}_{1}&\ldots&\mathbf{x}_{m-1}\\\
\vline&\vline&&\vline\end{bmatrix}$ (2)

and the time-shifted matrix

$\displaystyle\mathbf{X}^{\prime}$
$\displaystyle=\begin{bmatrix}\vline&\vline&&\vline\\\
\mathbf{x}_{1}&\mathbf{x}_{2}&\ldots&\mathbf{x}_{m}\\\
\vline&\vline&&\vline\end{bmatrix}.$ (3)

Typically $n\gg m$, i.e. there are many more spatial measurements available
than temporal. The vector $\mathbf{x}_{k}$ 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 $\mathbf{A}$, chosen to minimize
$\|\mathbf{x}_{k+1}-\mathbf{A}\mathbf{x}_{k}\|_{2}$ over the
${k=0,2,3,\cdots,{m-1}}$ snapshots, so that
$\mathbf{X}^{\prime}\approx\mathbf{A}\mathbf{X}$. Computationally, the matrix
$\mathbf{A}$ is obtained as

$\displaystyle\mathbf{A}=\mathbf{X}^{\prime}\mathbf{X}^{\dagger}$ (4)

where $\mathbf{A}\in\mathbb{R}^{n\times n}$ and $\mathbf{X}^{\dagger}$ is the
Moore-Penrose pseudo-inverse of $\mathbf{X}$. The dominant eigenvectors of
$\mathbf{A}$ are the dynamic modes $\boldsymbol{\Phi}$, and the associated
eigenvalues determine how these modes behave in time.

In practice, the high-dimensional $\mathbf{A}$ is not computed directly.
Instead, the SVD can be used to first project to a low-rank subspace and then
compute the matrix $\mathbf{\tilde{A}}=\mathbf{U}^{*}\mathbf{A}\mathbf{U}$
which has many of the same eigenvalues as $\mathbf{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].

Algorithm 2 DMD with control [56]

1:Data matrices $\mathbf{X}$, $\mathbf{X}^{\prime}$, input snapshot matrix
$\boldsymbol{\Upsilon}$, target rank $r$ of $\mathbf{X}$ or
$\mathbf{X}^{\prime}$ and $\tilde{r}$ of $\boldsymbol{\Omega}$.

2:Actuation matrix $\mathbf{B}$.

3:Spectrum $\boldsymbol{\Lambda}$, modes $\boldsymbol{\Phi}$, [ and actuation
matrix $\mathbf{\hat{B}}$]. $\triangleright$ Optional outputs in brackets
[$\cdot$].

4:procedure
DMDc($\mathbf{X},\mathbf{X}^{\prime},\boldsymbol{\Upsilon},r,\tilde{r},[\mathbf{B}]$)

5: if $\mathbf{B}$ is known then

6: $[\boldsymbol{\Lambda},\boldsymbol{\Phi}]\leftarrow$
DMD($\mathbf{X},\mathbf{X}^{\prime}-\mathbf{B}\boldsymbol{\Upsilon},r$)
$\triangleright$ Perform _DMD_ (Algorithm 1) adjusted

7: $\triangleright$ for known actuation.

8: else

9: $\boldsymbol{\Omega}\leftarrow\begin{bmatrix}\mathbf{X}\\\
\boldsymbol{\Upsilon}\end{bmatrix}$ $\triangleright$ Matrix of the state and
input snapshots.

10:
$[\mathbf{\tilde{U}},\boldsymbol{\tilde{\Sigma}},\mathbf{\tilde{V}}]\leftarrow$
SVD($\boldsymbol{\Omega},\tilde{r}$) $\triangleright$ Truncated
$\tilde{r}$-rank SVD of $\boldsymbol{\Omega}$.

11:
$\mathbf{\tilde{U}}_{1},\;\mathbf{\tilde{U}}_{2}\leftarrow\mathbf{\tilde{U}}$
$\triangleright$ Split $\mathbf{\tilde{U}}$ into two components.

12: $[\mathbf{\hat{U}},\boldsymbol{\hat{\Sigma}},\mathbf{\hat{V}}]\leftarrow$
SVD($\mathbf{X}^{\prime},r$) $\triangleright$ Truncated $r$-rank SVD of
$\mathbf{X}^{\prime}$.

13:
$\mathbf{\tilde{A}}\leftarrow\mathbf{\hat{U}}^{*}\mathbf{X}^{\prime}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}\mathbf{\tilde{U}}^{*}_{1}\mathbf{\hat{U}}$.
$\triangleright$ Low-rank approximation of $\mathbf{A}$.

14:
$\mathbf{\tilde{B}}\leftarrow\mathbf{\hat{U}}^{*}\mathbf{X}^{\prime}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}\mathbf{\tilde{U}}^{*}_{2}$
$\triangleright$ Estimate reduced actuation matrix $\mathbf{\tilde{B}}$.

15:
$\mathbf{\hat{B}}\leftarrow\mathbf{X}^{\prime}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}\mathbf{\tilde{U}}^{*}_{2}$
$\triangleright$ Estimate actuation matrix $\hat{\mathbf{B}}$.

16:
$[\mathbf{W},\boldsymbol{\Lambda}]\leftarrow\text{EIG}(\mathbf{\tilde{A}})$
$\triangleright$ Eigendecomposition of $\mathbf{\tilde{A}}$.

17:
$\boldsymbol{\Phi}\leftarrow\mathbf{X}^{\prime}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}\mathbf{\tilde{U}}^{*}_{1}\mathbf{\hat{U}}\mathbf{W}$
$\triangleright$ DMD modes of $\mathbf{A}$.

18: end if

19:end procedure

Note: If ${\lambda}_{i}=0$, then
$\boldsymbol{\phi}_{i}=\mathbf{\tilde{U}}_{1}\mathbf{\tilde{U}}^{*}_{1}\mathbf{\hat{U}}\mathbf{w}_{i}\text{
for step }13$.

###  2.1 Dynamic mode decomposition with control

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 $\mathbf{u}_{k}\in\mathbb{R}^{q}$

$\displaystyle\mathbf{x}_{k+1}=\mathbf{A}\mathbf{x}_{k}+\mathbf{B}\mathbf{u}_{k}$
(5)

where $\mathbf{B}\in\mathbb{R}^{n\times q}$ . The matrix form of the actuation
is

$\displaystyle\boldsymbol{\Upsilon}=\begin{bmatrix}\vline&\vline&&\vline\\\
\mathbf{u}_{0}&\mathbf{u}_{1}&\ldots&\mathbf{u}_{m-1}\\\
\vline&\vline&&\vline\end{bmatrix}.$ (6)

The system can be written in matrix form as:

$\displaystyle\mathbf{X}^{\prime}=\mathbf{A}\mathbf{X}+\mathbf{B}\boldsymbol{\Upsilon}$
(7)

where each column $\mathbf{u}_{k}$ of the input snapshot matrix
$\boldsymbol{\Upsilon}$ is the input at each snapshot in Eq. (6) and the data
matrices $\mathbf{X},\mathbf{X}^{\prime}$ are formulated in the same way as in
Eq. (2).

A least-squares regression algorithm can once again be used to determine the
matrix $\mathbf{A}$ and its associated DMD modes and eigenvalues. Two
distinguishing cases can be considered, when $\mathbf{B}$ is known and when
$\mathbf{B}$ is unknown. When the input matrix $\mathbf{B}$ is known, or can
be well-estimated, the output is a simple linear combination of states and
inputs. The DMD modes can then be obtained following the procedure of the
exact DMD algorithm 1 with $\mathbf{X}^{\prime}$ replaced with
$\mathbf{X}^{\prime}-\mathbf{B}\boldsymbol{\Upsilon}$. For an unknown
$\mathbf{B}$ it is possible to compute the DMD modes and an approximation to
the matrix $\mathbf{B}$ via regression [56]. For this case, an augmented
matrix containing both the state and input snapshots is constructed

$\displaystyle\boldsymbol{\Omega}={\begin{bmatrix}\mathbf{X}\\\
\boldsymbol{\Upsilon}\end{bmatrix}}$ (8)

along with an augmented matrix containing the two unknown system matrices

$\displaystyle\mathbf{G}=\begin{bmatrix}\mathbf{A}&\mathbf{B}\end{bmatrix}\,.$
(9)

The regression problem is then formulated as

$\displaystyle\mathbf{X}^{\prime}=\begin{bmatrix}\mathbf{A}&\mathbf{B}\end{bmatrix}\begin{bmatrix}\mathbf{X}\\\
\boldsymbol{\Upsilon}\end{bmatrix}=\mathbf{G}\boldsymbol{\Omega}$ (10)

where $\boldsymbol{\Omega}\in\mathbb{R}^{(n+q)\times m}$ is the combination of
the state and control snapshots. The DMDc algorithm with an unknown
$\mathbf{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 $\mathbf{A}$ and $\mathbf{B}$
through low-rank structures. Much like the DMD algorithm, DMDc relies on
least-squares regression of the data to build a linear model $\mathbf{A}$ of
the state dynamics that is disambiguated from the discovered actuation matrix
$\mathbf{B}$.

Algorithm 3 Compressive DMD [17]

1:Measurements $\mathbf{Y},\mathbf{Y}^{\prime}$, measurement matrix
$\mathbf{C}$, sparsifying basis $\boldsymbol{\Psi}$, and target rank $r$.

2:$\mathbf{X},\mathbf{X}^{\prime}$.

3:cDMD spectrum $\boldsymbol{\Lambda}$ and modes $\hat{\boldsymbol{\Phi}}$.

4:procedure
cDMD($\mathbf{Y},\mathbf{Y}^{\prime},\mathbf{C},r,[\mathbf{X},\mathbf{X}^{\prime}]$).

5:
$[\mathbf{U_{Y}},\boldsymbol{\Sigma_{\mathbf{Y}}},\mathbf{V_{Y}}]\leftarrow$
SVD($\mathbf{Y},r$) $\triangleright$ Truncated $r$-rank SVD of $\mathbf{Y}$.

6:
$\mathbf{\tilde{A}_{Y}}\leftarrow\mathbf{U}^{*}_{\mathbf{Y}}\mathbf{Y}^{\prime}\mathbf{V_{Y}}\boldsymbol{\Sigma}^{-1}_{\mathbf{Y}}$
$\triangleright$ Low-rank approximation of $\mathbf{A_{Y}}$.

7:
$[\mathbf{W_{Y}},\boldsymbol{\Lambda}_{\mathbf{Y}}]\leftarrow\text{EIG}(\mathbf{\tilde{A}_{Y}})$
$\triangleright$ Eigendecomposition of $\mathbf{\tilde{A}_{Y}}$.

8: if $\mathbf{X}$ is known then $\triangleright$ Perform _compressed DMD_.

9:
$\hat{\boldsymbol{\Phi}}\leftarrow\mathbf{X}^{\prime}\mathbf{V_{Y}}\boldsymbol{\Sigma}^{-1}_{\mathbf{Y}}\mathbf{W_{Y}}$
$\triangleright$ Estimate DMD modes of $\mathbf{A}$.

10: else$\triangleright$ Perform _compressed sensing DMD_.

11:
$\boldsymbol{\Phi_{\mathbf{Y}}}\leftarrow\mathbf{Y}^{\prime}\mathbf{V_{Y}}\boldsymbol{\Sigma}^{-1}_{\mathbf{Y}}\mathbf{W_{Y}}$
$\triangleright$ DMD modes of $\mathbf{A_{Y}}$.

12: $\boldsymbol{\Phi_{\mathbf{S}}}\leftarrow$ Compressed
Sensing($\boldsymbol{\Phi_{\mathbf{Y}}},\mathbf{C},\boldsymbol{\Psi}$)
$\triangleright$ Perform $l_{1}$ minimization on
$\boldsymbol{\phi_{\mathbf{Y}}}_{i}$ to solve for
$\boldsymbol{\phi_{\mathbf{S}}}_{i}$.

13:
$\hat{\boldsymbol{\Phi}}\leftarrow\boldsymbol{\Psi}\boldsymbol{\Phi_{\mathbf{S}}}$
$\triangleright$ Estimate DMD modes of $\mathbf{A}$.

14: end if

15:end procedure

Note: If ${\lambda}_{i}=0$, then
$\boldsymbol{\phi}_{i}=\mathbf{U}\mathbf{w}_{\mathbf{Y},i},\boldsymbol{\phi}_{\mathbf{Y},i}=\mathbf{U_{Y}}\mathbf{w}_{\mathbf{Y},i}$
for steps $6$ and $8$.

###  2.2 Compressed sensing and dynamic mode decomposition

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 $\mathbf{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].

Consider compressed or subsampled data $\mathbf{Y}$ given by

$\displaystyle\mathbf{Y}=\mathbf{C}\mathbf{X},$ (11)

where $\mathbf{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 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.

\begin{overpic}[width=216.81pt]{fig01.pdf} \small\put(20.0,54.0){$\mathbf{X}$,
$\mathbf{X}^{\prime}$} \put(20.0,36.0){$\mathbf{C}$}
\put(20.0,16.0){$\mathbf{Y}$, $\mathbf{Y}^{\prime}$} \put(46.0,57.0){DMD}
\put(43.0,51.0){$\mathbf{X}^{\prime}=\mathbf{A}\mathbf{X}$}
\put(72.5,54.0){$\boldsymbol{\Lambda}$, $\boldsymbol{\Phi}$}
\put(70.0,16.0){$\boldsymbol{\Lambda}_{\mathbf{Y}}$,
$\boldsymbol{\Phi_{\mathbf{Y}}}$} \put(45.0,40.0){Path 1:}
\put(36.0,36.0){Compressed DMD} \put(34.5,6.7){Path 2: Compressed sensing DMD}
\end{overpic} Figure 1: Schematic of DMD and compressive DMD. Path $1$ shows
compressed DMD and path $2$ shows compressed sensing DMD ( [17]).

In addition to the the data matrices $\mathbf{X}$ and $\mathbf{X}^{\prime}$,
it is also now required to specify a measurements matrix $\mathbf{C}$. These
three matrices together are required to execute Algorithm 3. In practice, we
often consider point measurements so that the rows of $\mathbf{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.

##  3 Compressive system identification

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.

We now consider a general input–output system with high-dimensional state
$\mathbf{x}$:

$\displaystyle\mathbf{x}_{k+1}$ $\displaystyle=$
$\displaystyle\mathbf{A}\mathbf{x}_{k}+\mathbf{B}\mathbf{u}_{k}$ (12a)
$\displaystyle\mathbf{y}_{k}$ $\displaystyle=$
$\displaystyle\mathbf{C}\mathbf{x}_{k}.$ (12b)

As in DMD, the goal is to obtain the leading eigendecomposition of
$\mathbf{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
$\mathbf{y}$ and disambiguate internal state dynamics from the effect of
actuation $\mathbf{u}$. Writing Eq. (12) in terms of data matrices yields:

$\displaystyle\mathbf{X}^{\prime}$ $\displaystyle=$
$\displaystyle\mathbf{A}\mathbf{X}+\mathbf{B}\boldsymbol{\Upsilon}$ (13a)
$\displaystyle\mathbf{Y}$ $\displaystyle=$
$\displaystyle\mathbf{C}\mathbf{X}.$ (13b)

Under certain conditions it is possible to apply DMDc to the compressed data
$\mathbf{Y}$ and then recover full-state DMD modes via compressed sensing. As
in the compressive DMD algorithm [17], if full-state data $\mathbf{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.

\begin{overpic}[width=216.81pt]{fig01} \small\put(14.5,54.0){$\mathbf{X}$,
$\mathbf{X}^{\prime}$, $\boldsymbol{\Upsilon}$, $[\mathbf{B}]$}
\put(20.0,36.0){$\mathbf{C}$} \put(14.5,16.0){$\mathbf{Y}$,
$\mathbf{Y}^{\prime}$, $\boldsymbol{\Upsilon}$, $[\mathbf{B}]$}
\put(45.0,57.0){DMDc}
\put(37.7,51.0){$\mathbf{X}^{\prime}=\mathbf{A}\mathbf{X}+\mathbf{B}\boldsymbol{\Upsilon}$}
\put(69.5,54.0){$\boldsymbol{\Lambda}$,
$\boldsymbol{\Phi},[\mathbf{\hat{B}}]$}
\put(65.0,16.0){$\boldsymbol{\Lambda}_{\mathbf{Y}}$,
$\boldsymbol{\Phi_{\mathbf{Y}}},[\mathbf{\hat{B}}_{\mathbf{Y}}]$}
\put(45.0,41.5){Path $1$:} \put(35.0,37.5){Compressed DMDc}
\put(33.0,6.8){Path $2$: Compressed sensing DMDc} \end{overpic} Figure 2:
Schematic of compressive DMD control. Path $1$ shows compressed DMDc and path
$2$ shows compressed sensing DMDc. Algorithm 4 Compressive DMD with control

1:Measurement snapshot matrices $\mathbf{Y},\mathbf{Y}^{\prime}$, sensing
matrix $\mathbf{C}$, input snapshot matrix $\boldsymbol{\Upsilon}$, and target
rank $r,\tilde{r}$.

2:Full-state data matrices $\mathbf{X},\mathbf{X}^{\prime}$ and actuation
matrix $\mathbf{B}$.

3:cDMDc spectrum $\boldsymbol{\Lambda}$ and modes
$\hat{\boldsymbol{\Phi}},[\mathbf{\hat{B}}]$.

4:procedure
cDMDc($\mathbf{Y},\mathbf{Y}^{\prime},\mathbf{C},r,\tilde{r},[\mathbf{X},\mathbf{X}^{\prime},\mathbf{B}]$)

5: if $\mathbf{X}$ is known then

6: if $\mathbf{B}$ is known then

7: $\boldsymbol{\Lambda},\hat{\boldsymbol{\Phi}}\leftarrow$
cDMD($\mathbf{Y},\mathbf{Y}^{\prime}-\mathbf{C}\mathbf{B}\boldsymbol{\Upsilon},\mathbf{C},\mathbf{X},\mathbf{X}^{\prime},r$)
$\triangleright$ Perform _compressed DMD_.

8: else

9:
$\boldsymbol{\Lambda}_{\mathbf{Y}},\mathbf{\hat{U}_{Y}},\mathbf{\tilde{U}}_{\mathbf{Y},1},\mathbf{\tilde{U}}_{\mathbf{Y},2},\boldsymbol{\tilde{\Sigma}}_{\mathbf{Y}},\mathbf{\tilde{V}_{Y}},\mathbf{W_{Y}}$

10: $\leftarrow$
DMDc($\mathbf{Y},\mathbf{Y}^{\prime},\boldsymbol{\Upsilon},r,\tilde{r},\mathbf{X},\mathbf{X}^{\prime}-\mathbf{B}\boldsymbol{\Upsilon}$)
$\triangleright$ Perform _DMDc_.

11:
$\hat{\boldsymbol{\Phi}}\leftarrow\mathbf{X}^{\prime}\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}^{*}_{\mathbf{Y},1}\mathbf{\hat{U}_{Y}}\mathbf{W_{Y}}$
$\triangleright$ Estimate DMD modes of $\mathbf{A}$.

12:
$\mathbf{\hat{B}}\leftarrow\mathbf{X}^{\prime}\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}^{*}_{\mathbf{Y},2}$
$\triangleright$ Estimate actuation matrix $\hat{\mathbf{B}}$.

13: end if

14: else

15: if $\mathbf{B}$ is known then

16: $\boldsymbol{\Lambda},\hat{\boldsymbol{\Phi}}\leftarrow$
cDMD($\mathbf{Y},\mathbf{Y}^{\prime}-\mathbf{C}\mathbf{B}\boldsymbol{\Upsilon},\mathbf{C},r$)
$\triangleright$ Perform _compressed sensing DMD_.

17: else

18:
$\boldsymbol{\Lambda}_{\mathbf{Y}},\boldsymbol{\Phi_{\mathbf{Y}}},\mathbf{\hat{B}}_{\mathbf{Y}}\leftarrow$
DMDc($\mathbf{Y},\mathbf{Y}^{\prime},\boldsymbol{\Upsilon},r,\tilde{r}$)
$\triangleright$ Perform _DMDc_.

19: $\hat{\boldsymbol{\Phi}}\leftarrow$ Compressed
Sensing($\boldsymbol{\Phi_{\mathbf{Y}}},\boldsymbol{\Psi}$) $\triangleright$
Estimate DMD modes of $\mathbf{A}$.

20: $\mathbf{\hat{B}}\leftarrow$ Compressed
Sensing($\mathbf{\hat{B}}_{\mathbf{Y}},\boldsymbol{\Psi}$) $\triangleright$
Estimate actuation matrix $\hat{\mathbf{B}}$.

21: end if

22: end if

23:end procedure

Note: If $\lambda_{i}=0$, then
$\boldsymbol{\phi}_{i}=\mathbf{\tilde{U}}\mathbf{\tilde{U}}^{*}_{\mathbf{Y},1}\mathbf{\hat{U}_{Y}}\mathbf{w}_{\mathbf{Y},i}\text{
for step }7.$

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

######  Assumption 1.

The measurement matrix $\mathbf{C}$ preserves the temporal information in
$\boldsymbol{\Omega}$ so that
$\mathbf{\tilde{V}_{Y}}\mathbf{\tilde{V}_{Y}}^{*}\mathbf{\tilde{V}}=\mathbf{\tilde{V}}$.
This requires the columns of $\mathbf{\tilde{V}}$ to be in the column space of
$\mathbf{\tilde{V}_{Y}}$ and will only be approximately satisfied with
measurement noise.

######  Assumption 2.

The columns of the full-state output matrix $\mathbf{X}^{\prime}$ are in the
subspace of the upper left singular vectors $\mathbf{\tilde{U}}_{1}$ of the
full-state augmented matrix $\boldsymbol{\Omega}$ so that
$\mathbf{\tilde{U}}_{1}\mathbf{\tilde{U}}^{*}_{1}\mathbf{X}^{\prime}\approx\mathbf{X}^{\prime}$.

######  Lemma 1.

If Assumption 1 holds, we have an identity
$\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}=\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},1}^{*}\mathbf{C}\mathbf{\tilde{U}}_{1}$,
similar to that derived in [17]:

###### Proof.

$\displaystyle\mathbf{Y}$ $\displaystyle=\mathbf{C}\mathbf{X}$
$\displaystyle\mathbf{\tilde{U}}_{\mathbf{Y},1}\boldsymbol{\tilde{\Sigma}_{\mathbf{Y}}}\mathbf{\tilde{V}_{Y}}^{*}$
$\displaystyle=\mathbf{C}\mathbf{\tilde{U}}_{1}\boldsymbol{\tilde{\Sigma}}\mathbf{\tilde{V}}^{*}$
$\displaystyle\mathbf{\tilde{V}}^{*}_{\mathbf{Y}}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}$
$\displaystyle=\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},1}^{*}\mathbf{C}\mathbf{\tilde{U}}_{1}$
$\displaystyle\mathbf{\tilde{V}_{Y}}\mathbf{\tilde{V}_{Y}}^{*}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}$
$\displaystyle=\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},1}^{*}\mathbf{C}\mathbf{\tilde{U}}_{1}$
$\displaystyle\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}$
$\displaystyle=\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},1}^{*}\mathbf{C}\mathbf{\tilde{U}}_{1}.$

∎

The next theorem uses the following definitions of $\mathbf{A}$ and
$\mathbf{A_{Y}}$ from the DMDc and cDMDc algorithms:

$\displaystyle\mathbf{A}$
$\displaystyle=\mathbf{X}^{\prime}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}\mathbf{\tilde{U}}^{*}_{1}$
$\displaystyle\mathbf{A_{Y}}$
$\displaystyle=\mathbf{Y}^{\prime}\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},1}^{*}.$

######  Theorem 1.

If Assumption 1 and 2 hold, then

$\displaystyle\mathbf{C}\mathbf{A}\mathbf{X}^{\prime}=\mathbf{A_{Y}}\mathbf{C}\mathbf{X}^{\prime}.$
(14)

###### Proof.

Using Lemma 1, we have

$\displaystyle\mathbf{C}\mathbf{A}\mathbf{X}^{\prime}$
$\displaystyle=\mathbf{C}(\mathbf{X}^{\prime}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}\mathbf{\tilde{U}}^{*}_{1})\mathbf{X}^{\prime}$
$\displaystyle=\mathbf{Y}^{\prime}(\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},1}^{*}\mathbf{C}\mathbf{\tilde{U}}_{1})\mathbf{\tilde{U}}^{*}_{1}\mathbf{X}^{\prime}$
$\displaystyle=\mathbf{A_{Y}}\mathbf{C}\mathbf{\tilde{U}}_{1}\mathbf{\tilde{U}}^{*}_{1}\mathbf{X}^{\prime}$
$\displaystyle=\mathbf{A_{Y}}\mathbf{C}\mathbf{X}^{\prime}.$

∎

Thus, the compression matrix $\mathbf{C}$ commutes with the dynamics in
$\mathbf{A}$, when applied to data $\mathbf{X}^{\prime}$. We use this to
obtain the central result: compressed dynamic modes of the full data are
dynamic modes of the compressed data.

######  Theorem 2.

If Assumptions 1 and 2 hold, then compressing a full-state DMDc eigenvector
$\boldsymbol{\phi}$ yields a DMDc eigenvector of the compressed data with the
same eigenvalue:

$\displaystyle\mathbf{A_{Y}}\mathbf{C}\boldsymbol{\phi}=\lambda\mathbf{C}\boldsymbol{\phi}.$
(15)

###### Proof.

$\displaystyle\mathbf{A_{Y}}\mathbf{C}\boldsymbol{\phi}$
$\displaystyle=\mathbf{A_{Y}}\mathbf{C}\mathbf{X}^{\prime}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}\mathbf{\tilde{U}}_{1}^{*}\mathbf{\hat{U}}\mathbf{w}$
$\displaystyle=\mathbf{C}\mathbf{A}\mathbf{X}^{\prime}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}\mathbf{\tilde{U}}_{1}^{*}\mathbf{\hat{U}}\mathbf{w}$
$\displaystyle=\mathbf{C}\mathbf{A}\boldsymbol{\phi}$
$\displaystyle=\lambda\mathbf{C}\boldsymbol{\phi}.$

∎

If $\mathbf{C}$ is chosen poorly so that $\boldsymbol{\phi}$ is in its
nullspace, then Theorem 2 applies trivially. This theorem does not guarantee
that every eigenvector of $\mathbf{A_{Y}}$ is a compressed eigenvector of
$\mathbf{A}$, although under reasonable assumptions the dominant eigenvalues
of $\mathbf{A_{Y}}$ will approximate those of $\mathbf{A}$, as shown in [17].
We then have

$\displaystyle\boldsymbol{\phi}_{\mathbf{Y}}=\mathbf{C}\boldsymbol{\phi}_{\mathbf{X}}=\mathbf{C}\boldsymbol{\Psi}\boldsymbol{\phi}_{\mathbf{S}},$
(16)

where $\boldsymbol{\phi}_{\mathbf{S}}$ is the sparse representation of
$\boldsymbol{\phi}_{\mathbf{X}}$ in a universal basis $\boldsymbol{\Psi}$,
such as a Fourier or wavelet basis. Thus, it is possible to recover
$\boldsymbol{\phi}_{\mathbf{S}}$, and hence $\boldsymbol{\phi}_{\mathbf{X}}$,
from compressed DMD modes $\boldsymbol{\phi}_{\mathbf{Y}}$, given enough
incoherent measurements [17].

##### Compressed recovery of the actuation matrix.

We now establish a similar relationship between the full-state actuation
matrix $\mathbf{B}$ and the compressed matrix $\mathbf{B}_{\mathbf{Y}}$.

######  Assumption 3.

The columns of $\mathbf{\tilde{U}}_{\mathbf{Y},2}$ are spanned by those of
$\mathbf{\tilde{U}}_{2}$, so
$\mathbf{\tilde{U}}_{2}\mathbf{\tilde{U}}_{2}^{*}\mathbf{\tilde{U}}_{\mathbf{Y},2}=\mathbf{\tilde{U}}_{\mathbf{Y},2}$
and
$\mathbf{\tilde{U}}_{\mathbf{Y},2}^{*}\mathbf{\tilde{U}}_{2}\mathbf{\tilde{U}}_{2}^{*}=\mathbf{\tilde{U}}_{\mathbf{Y},2}^{*}$.

######  Lemma 2.

We have an identity
$\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}=\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},2}^{*}\mathbf{\tilde{U}}_{2}$.

###### Proof.

Expanding $\boldsymbol{\Upsilon}$ in the SVD bases of
$\boldsymbol{\Omega}_{\mathbf{Y}}$ and $\boldsymbol{\Omega}$, we find:

$\displaystyle\mathbf{\tilde{U}}_{\mathbf{Y},2}\boldsymbol{\tilde{\Sigma}_{\mathbf{Y}}}\mathbf{\tilde{V}_{Y}}^{*}$
$\displaystyle=\mathbf{\tilde{U}}_{2}\boldsymbol{\tilde{\Sigma}}\mathbf{\tilde{V}}^{*}$
$\displaystyle\mathbf{\tilde{V}}^{*}_{\mathbf{Y}}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}$
$\displaystyle=\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},2}^{*}\mathbf{\tilde{U}}_{2}$
$\displaystyle\mathbf{\tilde{V}_{Y}}\mathbf{\tilde{V}_{Y}}^{*}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}$
$\displaystyle=\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},2}^{*}\mathbf{\tilde{U}}_{2}$
$\displaystyle\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}$
$\displaystyle=\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},2}^{*}\mathbf{\tilde{U}}_{2}.$

∎

The next theorem uses the following definitions of $\mathbf{B}$ and
$\mathbf{B}_{\mathbf{Y}}$ from the DMDc and cDMDc algorithms:

$\displaystyle\mathbf{B}$
$\displaystyle=\mathbf{X}^{\prime}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}\mathbf{\tilde{U}}^{*}_{2}$
$\displaystyle\mathbf{B}_{\mathbf{Y}}$
$\displaystyle=\mathbf{Y}^{\prime}\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},2}^{*}.$

######  Theorem 3.

If Assumptions 1 and 3 hold, then

$\displaystyle\mathbf{C}\mathbf{B}=\mathbf{B}_{\mathbf{Y}}.$ (17)

###### Proof.

Using Lemma 2, we have

$\displaystyle\mathbf{C}\mathbf{B}$
$\displaystyle=\mathbf{C}\mathbf{X}^{\prime}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}\mathbf{\tilde{U}}^{*}_{2}$
$\displaystyle=\mathbf{Y}^{\prime}\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},2}^{*}\mathbf{\tilde{U}}_{2}\mathbf{\tilde{U}}^{*}_{2}$
$\displaystyle=\mathbf{Y}^{\prime}\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},2}^{*}$
$\displaystyle=\mathbf{B}_{\mathbf{Y}}.$

∎

Therefore, we can first compute the DMD modes $\boldsymbol{\Phi_{\mathbf{Y}}}$
and the compressed actuation matrix $\mathbf{B}_{\mathbf{Y}}$ and then
reconstruct the full-state $\boldsymbol{\Phi}$ and $\mathbf{B}$ through
compressed sensing.

###  3.1 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].

######  Theorem 4.

If Assumptions 1 and 2 hold, then

$\displaystyle\mathbf{C}\mathbf{A}\mathbf{B}=\mathbf{A_{Y}}\mathbf{C}\mathbf{B}.$
(18)

###### Proof.

Using Lemma 1, we have

$\displaystyle\mathbf{C}\mathbf{A}\mathbf{B}$
$\displaystyle=\mathbf{C}(\mathbf{X}^{\prime}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}\mathbf{\tilde{U}}^{*}_{1})\mathbf{B}$
$\displaystyle=\mathbf{Y}^{\prime}(\mathbf{\tilde{V}_{Y}}\boldsymbol{\tilde{\Sigma}}^{-1}_{\mathbf{Y}}\mathbf{\tilde{U}}_{\mathbf{Y},1}^{*}\mathbf{C}\mathbf{\tilde{U}}_{1})\mathbf{\tilde{U}}^{*}_{1}\mathbf{B}$
$\displaystyle=\mathbf{A_{Y}}\mathbf{C}{\mathbf{\tilde{U}}_{1}\mathbf{\tilde{U}}^{*}_{1}\mathbf{X}^{\prime}}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}\mathbf{\tilde{U}}^{*}_{2}$
$\displaystyle=\mathbf{A_{Y}}\mathbf{C}\mathbf{X}^{\prime}\mathbf{\tilde{V}}\boldsymbol{\tilde{\Sigma}}^{-1}\mathbf{\tilde{U}}^{*}_{2}$
$\displaystyle=\mathbf{A_{Y}}\mathbf{C}\mathbf{B}.$

∎

######  Corollary 1.

Theorem 4 can be expanded to further steps in dynamics, for
$k\in\mathbb{Z}^{+}$, such as:

$\displaystyle\mathbf{C}\mathbf{A}^{k}\mathbf{B}=\mathbf{A_{Y}}^{k}\mathbf{C}\mathbf{B}=\mathbf{A_{Y}}^{k}\mathbf{B}_{\mathbf{Y}}.$
(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.

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

$\displaystyle{\mathcal{C}}=\begin{bmatrix}\mathbf{B}&\mathbf{A}\mathbf{B}&\cdots\mathbf{A}^{n-1}\mathbf{B}\end{bmatrix}$
(20)

we have the following relationship:

$\displaystyle\mathbf{C}\,{\mathcal{C}}$
$\displaystyle=\mathbf{C}\begin{bmatrix}\mathbf{B}&\mathbf{A}\mathbf{B}&\cdots\mathbf{A}^{n-1}\mathbf{B}\end{bmatrix}$
(21a)
$\displaystyle=\begin{bmatrix}\mathbf{B}_{\mathbf{Y}}&\mathbf{A}_{\mathbf{Y}}\mathbf{B}_{\mathbf{Y}}&\cdots&\mathbf{A}^{n-1}_{\mathbf{Y}}\mathbf{B}_{\mathbf{Y}}\end{bmatrix}=\mathcal{C}_{\mathbf{Y}}.$
(21b)

Therefore, if the full-state systems is controllable, i.e. the controllability
matrix ${\mathcal{C}}$ is full rank, and the compressed measurement matrix
$\mathbf{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].

###  3.2 Computational considerations

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

  1. 1.

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

  2. 2.

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

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

In Algorithm 4, these four scenarios are discussed. When both $\mathbf{X}$ and
$\mathbf{B}$ are known, the spectrum $\boldsymbol{\Lambda}_{\mathbf{Y}}$ and
modes $\boldsymbol{\Phi}$ are obtained by performing compressed DMD on the
pre-compressed data $\mathbf{Y}$ and the shifted matrix correcting for the
effect of control
$\mathbf{Y}^{\prime}-\mathbf{C}\mathbf{B}\boldsymbol{\Upsilon}$. When
$\mathbf{X}$ is known and $\mathbf{B}$ is unknown, DMDc is computed on
$\mathbf{Y},\mathbf{Y}^{\prime}$ and the dynamic modes $\boldsymbol{\Phi}$ and
actuation matrix $\mathbf{B}$ are reconstructed as a linear combination of the
full-state data $\mathbf{X}^{\prime}$. When $\mathbf{X}$ is only partially
known and $\mathbf{B}$ is known, compressed sensing DMD is performed on
$\mathbf{Y}$ and
$\mathbf{Y}^{\prime}-\mathbf{C}\mathbf{B}\boldsymbol{\Upsilon}$ and the full-
state modes $\boldsymbol{\Phi}$ are reconstructed from the compressed
$\boldsymbol{\Phi_{\mathbf{Y}}}$. When both $\mathbf{X}$ and $\mathbf{B}$ are
unavailable, DMDc is computed on $\mathbf{Y},\mathbf{Y}^{\prime}$ and the
dynamic modes $\boldsymbol{\Phi}$ and actuation matrix $\mathbf{B}$ are
reconstructed using compressed sensing.

Figure 3: Two-dimensional system with known dynamics: (a) phase portrait
(color denotes the progression of time), (b) x-t diagram of inflated high-
dimensional system, $\mathbf{X}\in\mathbb{R}^{1024\times 301}$ (first 31
snapshots shown), (c) two orthogonal modes of $\mathbf{P}$ and (d) DCT
coefficients of the two modes of $\mathbf{P}$.

##  4 Results

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

###  4.1 Stochastically forced linear system

This experiment is designed to test the compressive DMD with control framework
on an example where the low-rank dynamics are known. Thus, it is possible to
directly compare the true eigenvalue spectrum and spatiotemporal modes with
those obtained via compressed DMDc and compressed sensing DMDc.

The state matrix $\mathbf{\tilde{A}}$ and input matrix $\mathbf{\tilde{B}}$
are designed to yield a stable, controllable system:

$\displaystyle\begin{bmatrix}\tilde{x}^{1}\\\
\tilde{x}^{2}\end{bmatrix}_{k+1}$
$\displaystyle=\underbrace{\begin{bmatrix}0.9&0.2\\\
-0.1&0.9\end{bmatrix}}_{\mathbf{\tilde{A}}}\begin{bmatrix}\tilde{x}^{1}\\\
\tilde{x}^{2}\end{bmatrix}_{k}+\underbrace{\begin{bmatrix}0.1\\\
0.01\end{bmatrix}}_{\mathbf{\tilde{B}}}u_{k}.$ (22)

The dynamics are excited via Gaussian random input excitation in $u_{k}$,
starting from an initial condition
$\mathbf{x}_{0}=\begin{bmatrix}0.25&0.25\end{bmatrix}^{T}$. The system is
integrated from $t=0$ to $t=30$ with a time-step of $\Delta t=0.1$, resulting
in $301$ snapshots. A sample trajectory of the low-dimensional system is shown
in Fig. 3 (a).

| DMDc | c-DMDc | cs-DMDc
---|---|---|---
$\mathbf{B}$ known |  $\mathbf{C}$-type $1$ | 3.631e-13 | 3.627e-13 | 3.627e-13
$\mathbf{C}$-type $2$ | 3.631 e-13 | 3.443 e-13 | 3.443 e-13
$\mathbf{C}$-type $3$ | 3.631e-13 | 4.210e-13 | 4.210e-13
$\mathbf{B}$ unknown |  $\mathbf{C}$-type $1$ |  4.481e-131.758e-16 |  5.396e-134.840e-17 |  5.926e-139.558e-17
$\mathbf{C}$-type $2$ |  4.481e-131.758e-16 |  4.159e-134.731e-17 |  3.872e-132.845e-16
$\mathbf{C}$-type $3$ |  4.481e-131.758e-16 |  4.322e-133.022e-16 |  5.691e-132.785e-16
Table 1: Normalized error of
$\|\boldsymbol{\Phi}-\hat{\boldsymbol{\Phi}}\|_{F}$ and
$\|\mathbf{B}-\hat{\mathbf{B}}\|_{2}$ using DMDc, compressed DMDc, and
compressed sensing DMDc. Three types of compression are shown: uniform random
projections ($\mathbf{C}$-type $1$), Gaussian random projections
($\mathbf{C}$-type $2$) and single pixel measurements ($\mathbf{C}$-type $3$).
When $\mathbf{B}$ is unknown, the error of the estimated $\hat{\mathbf{B}}$ is
given in the upper triangle, and the error of $\hat{\boldsymbol{\Phi}}$ is
given in the lower triangle. In all cases, $\mathbf{B}$ is a linear
combination of columns of $\mathbf{P}$.

To inflate the state dimension, we associate each of the two states
$\tilde{x}^{1}$ and $\tilde{x}^{2}$ with a high-dimensional mode that is
sparse in the spatial wavenumber domain. These modes, shown in Fig. 3 (c), are
given by the columns of $\mathbf{P}\in\mathbb{R}^{1024\times 2}$, where each
column is constructed to have $K=4$ non-zero elements in the discrete cosine
transform (DCT) basis. The nonzero DCT coefficients of each mode have the same
wave number with different magnitudes, as shown in Fig. 3 (d). Thus, it is
possible to use the DCT matrix as the sparsifying basis for compressed
sensing. With the spatial modes in $\mathbf{P}$, it is possible to lift the
low-dimensional state $\mathbf{\tilde{x}}_{k}$ to a high-dimensional state
$\mathbf{x}_{k}=\mathbf{P}\mathbf{\tilde{x}}_{k}$. The lifted state is shown
in Fig. 3 (b) for the first $31$ snapshots, from $t=0$ to $t=3$. Note that the
high-dimensional actuation vector $\mathbf{B}$ is chosen to be in the span of
the columns of $\mathbf{P}$, i.e. $\mathbf{B}=\mathbf{P}\tilde{\mathbf{B}}$,
so that the actuation only excites low-dimensional dynamics; other types of
actuation will be examined further in Fig. 6.

To investigate the proposed cDMDc algorithm, we now consider compressed
measurements $\mathbf{y}$ given by Eq. (12b). Specifically, the compression
matrix $\mathbf{C}$ can be built with entries drawn from uniform
($\mathbf{C}$-type 1), Gaussian ($\mathbf{C}$-type 2) or Bernoulli
($\mathbf{C}$-type 3) distributions. Fig. 4 shows the DMDc mode reconstruction
using compressed DMDc (Path 1 in Fig. 2) and compressed sensing DMDc (Path 2
in Fig. 2) for $p=128$ Gaussian measurements. In both cases, it is assumed
that the high-dimensional actuation input vector ${\mathbf{B}}$ is known, and
the reconstructed modes faithfully reproduce the true coherent structures of
the underlying system (i.e., the DMD modes are a linear combination of the
columns of $\mathbf{P}$). The results remain unchanged when $\mathbf{B}$ is
unknown and must also be reconstructed.

Figure 4: Reconstruction of (a) actuation vector $\mathbf{B}$, (b)-(c) modes
from DMDc, compressed DMDc and compressed sensing DMDc. For compressive DMD,
$128$ Gaussian random measurements are collected from $1024$ state dimensions.
Figure 5:  Noise dependency of DMD spectrum for the true system, DMDc and
cDMDc based on $100$ different noise realizations. Rows correspond to
different noise levels $\eta\in\\{0.1,0.25,0.5\\}$ with
$\sigma_{noise}=\eta\,\mathrm{max}(\sigma_{i})$, where $\sigma_{i}$ denotes
the standard deviation of each spatial measurement. Figure 6: Estimated
$\hat{\mathbf{B}}$ and mode magnitudes for different choices of $\mathbf{B}$
for DMDc, compressed DMDc and compressed sensing DMDc. For compressive DMDc,
the same number of measurements is employed as for Fig. 4.

Table 1 shows the error in the reconstructed eigenfunctions
$\hat{\boldsymbol{\Phi}}$ and estimated actuation vector $\hat{\mathbf{B}}$
(when unknown) for DMDc, compressed DMDc, and compressed sensing DMDc,
compared against the true values. In addition, we investigate the effect of
different sensing strategies discussed above. In all cases, the reconstruction
error is small, as no noise is added to the simulated data.

Compressed sensing DMDc is able to uncover the underlying dynamics and spatio-
temporal modes from noiseless subsampled data, relaxing the requirement of
high-dimensional measurements. However, in realistic experimental conditions,
measurement noise will always be present and is known to effect DMD
computations. Figure 5 shows the performance of DMDc and compressed DMDc for
varying levels of measurement noise, averaged over 100 different noise
realizations in each case. Note that compressed sensing DMDc and compressed
DMDc have identical eigenvalues, as both methods compute the DMD spectrum from
the same compressed data. For moderate noise levels, the compressed DMDc
algorithm provides reasonably accurate eigenvalues, although they are less
accurate than those predicted by DMDc. As the noise intensity is increased to
greater than 50%, both the DMDc and cDMDc eigenvalues begin to deviate, and in
some cases, complex conjugate eigenvalues are miscomputed as purely real
eigenvalues. Although the effect of noise is exacerbated by compression,
sensitivity of DMD eigenvalues with measurement noise is a known issue, and
has been extensively studied and characterized [4, 31, 23]. There are a number
of algorithmic extensions that improve the eigenvalue prediction with noise,
including using the total least squares [31], a forward-backward symmetrizing
algorithm [23], subspace DMD [68], Bayesian DMD [67], or an optimal DMD based
on variable projection methods [2]. Each of these methods may be effectively
combined with the proposed cDMDc architecture to yield more accurate
eigenvalues.

####  4.1.1 Effect of Actuation

Finally, we investigate the performance of cDMDc for different choices of the
actuation vector $\mathbf{B}$ in Fig. 6. In addition to the cases presented so
far, where $\mathbf{B}$ is in the subspace of $\mathbf{P}$, we now explore
scenarios when $\mathbf{B}$ is randomly generated or in the complementary
subspace of $\mathbf{P}$. For all types of actuation $\mathbf{B}$, and
regardless of whether or not $\mathbf{B}$ is known, compressed DMDc and
compressed sensing DMDc both accurately identify the true modes
$\boldsymbol{\phi}_{1}$ and $\boldsymbol{\phi}_{2}$, which is also consistent
with DMDc. When $\mathbf{B}$ is unknown, cDMDc and csDMDc both accurately
identify $\mathbf{B}$, regardless of the subspace it belongs to, as long as
the columns of $\mathbf{B}$ are sparse. The algorithms do not accurately
capture the $\mathbf{B}$ matrix when it is not sparse (i.e., random actuation
vector); however, in this case, DMDc also misidentifies the actuation vector.

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 filled circles depict the time instants of the shown vorticity plots in
(a)-(c).

###  4.2 Fluid flow past a pitching plate

In the second example, we apply the proposed compressive DMD with control
algorithm to model the vorticity field downstream of a pitching plate at
Reynolds number $Re=100$. This pitching airfoil has been studied previously in
the context of reduced-order models for flow control [66, 22, 18, 14, 30]. The
flow is simulated using the immersed boundary projection method (IBPM)222Code
available at https://github.com/cwrowley/ibpm method [66, 22] with a grid
resolution of $799\times 159$ on a domain of size $10\times 2$,
nondimensionalized by the plate length $L$. The flow is simulated with a time-
step of $\Delta t=0.01$ dimensionless convective time units,
nondimensionalized by the length $L$ and free-stream velocity $U$. $251$
snapshots are sampled at a rate of $\Delta t=0.1$. In this example, the
airfoil is rapidly pitched up and down between $\pm 5^{\circ}$ at irregular
intervals in time, using the canonical pitch maneuver described in [55]. Six
rapid pitch maneuvers are performed, and then the data is symmetrized by
concatenating a mirror image of these $251$ snapshots with the opposite signed
vorticity, in an attempt to identify unbiased modes. Figure 7 illustrates the
vorticity field at different times, $t=0,11.5,21.5$. We focus on a downstream
measurement window of the size $399\times 141$ to mimic a PIV window. The
actuation input takes $4.6$ convective time units to reach the observation
window, and the actuation input, that used as input for cDMDc and DMDc, is
shifted accordingly. Similar small amplitude pitching motions have been shown
to be well-approximated with linear models [18].

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.

The eigenvalues of $\mathbf{\tilde{A}}$ and $\mathbf{\tilde{A}_{Y}}$ given by
DMDc and cDMDc are shown in Fig. 8. The cloud of cDMDc eigenvalues is
generated from an ensemble of $50$ realizations using $10\%$ compressed
measurements (i.e., $p=0.1n$) based on Gaussian random projections. We use
$r=9$ in this study to identify the spatiotemporal coherent DMD modes, and
they capture $96\%$ of the total energy, based on the singular values. These
nine modes are ordered by their decay rate (i.e., the real part) and the
higher-order modes have larger frequencies of oscillation. The discrete-time
eigenvalues are located close to the unit circle, and they are slightly
stable, as seen in the continuous-time plot. The DMDc eigenvalues appear to
have a decay rate that is too small, although the dynamics of transients are
captured quite well. Nevertheless, the cDMDc eigenvalues agree well with the
DMDc eigenvalues, which is the goal.

In Fig. 9, we present the spatial structure of the DMDc mode pairs and the
$\mathbf{B}$ matrix. Similar DMD modes have been observed in the flow past a
cylinder [20]. As the temporal frequency associated with an eigenvalue
increases, the modes are characterized by higher spatial wavenumbers, which is
characteristic of bluff-body flows. The accuracy of mode reconstruction is
similar for Gaussian random projections and single-pixel measurements, as
indicated in Fig. 11. However, single pixel measurements may be less expensive
and more realistic in real-world applications, and sampling $10\%$ of the
original measurements results in a reconstruction accuracy of $90\%$ in this
example.

Figure 9: Estimation of mode pairs (real and imaginary parts) and $\mathbf{B}$
using DMDc, compressed DMDc and compressed sensing DMDc. The cDMDc modes and
actuation matrix $\mathbf{B}$ (two columns, $\mathbf{B}_{1}$ and
$\mathbf{B}_{2}$) 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 are very similar (not shown).

Taking the DMDc results as the reference, Fig. 11 shows the error of the
eigenvalues with increasing number of measurements used in cDMDc. Note that
the eigenvalues are the same for compressed DMDc and compressed sensing DMDc,
as shown in Algorithm 4 . The eigenvalue ${\lambda}_{0}$ corresponding to zero
frequency is estimated with the greatest accuracy from the fewest
measurements, presumably because this mode contains the most energy in the
flow. The error in eigenvalues associated with other modes decreases
logarithmically with increasing compression ratio. Single pixel measurements
have similar performance compared with Gaussian random measurements for small
compression ratios. When $p=n$, the error goes to zero for single pixel
measurement, since the measurement is an invertible permutation of the
identity matrix, and cDMDc is equivalent to DMDc in this case.

Successful reconstruction using compressed sensing relies on the sparsity of
the state in some transform basis. Indeed, both $\boldsymbol{\Phi}$ and
$\mathbf{B}$ are sparse in the DCT basis. In general, the DMDc modes are more
sparse than the actuation matrix, and we choose $K=300$ to ensure a good
reconstruction of the modes and $\mathbf{B}$ achieving an error of about
$1\%$. In particular, we use the CoSaMP algorithm [54] to perform the
$l_{1}$-minimization with $10$ iterations and the desired sparsity of $K=300$.
The $L^{2}$ errors between the compressive DMDc and DMDc modes are shown in
Fig. 11. The convergence of error versus compression ratio in compressed
sensing DMDc is much slower than that in compressed DMDc, especially for the
zero frequency mode.

Figure 10: Error of estimated modes $\hat{\boldsymbol{\phi}}_{i}$ and
actuation matrix $\hat{\mathbf{B}}$ for increasing compression ratio using
compressed DMDc and compressed sensing DMDc for (a) Gaussian random
projections and (b) single pixel measurements. The DMDc modes and actuation
matrix are used as the reference, denoted by $\boldsymbol{\phi}$ and
$\mathbf{B}$, respectively.

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 $\lambda_{i}$.

The necessary number of measurements is given from theory to be $p\sim
4K\log_{10}(n/K)\approx 2728$, which corresponds to a compression ratio of
$5\%$. This matches the observation that the error curves plateau at around
$p/n=0.05$. Overall, the error for compressed measurements using Gaussian
random projections is comparable with using single pixel measurements.

##  5 Conclusions

In summary, we have presented a unified framework for compressive system
identification based on the dynamic mode decomposition. First, we describe the
previously developed compressed sensing DMD (csDMD) and DMD with control
(DMDc) algorithms in a common mathematical framework, providing algorithmic
implementation details. Next, we show how it is possible to construct reduced-
order models from compressed measurements and then reconstruct full-state
modes corresponding to the reduced states via compressed sensing. This lifting
procedure adds interpretability to otherwise black-box models. The compressed
DMD with control algorithm is demonstrated on two example systems, including a
high-dimensional discretized simulation of fluid flow past a pitching airfoil.
In both cases, accurate modal decompositions are achieved with surprisingly
few measurements, showcasing the efficacy of the proposed method. In addition,
we have released our entire code base to promote reproducible research and
reduce the barrier to implement these methods.

There are a number of important extensions and future directions that arise
out of this work. First, it will be interesting to further investigate the
relationship between the controllable and observable subspaces and full-state
recovery via compressed sensing. The goal is a generalized theory that
combines the notion of controllability and observability, based on the
structure of the $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ matrices, and
the notion of sparse signal recovery, via the structure of the low-rank
embedding $\mathbf{P}$. It may also be possible to extend results to nonlinear
estimation [64, 65] and control [57, 41, 35] through the connection 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.

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.

## Acknowledgements

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 high-dimensional
fluids, including Bing Brunton, Krithika Manohar, Igor Mezić, Bernd Noack, and
Sam Taira.

## References

  * [1] H. Arbabi and I. Mezić, Ergodic theory, dynamic mode decomposition and computation of spectral properties of the koopman operator, arXiv preprint arXiv:1611.06664, (2016).
  * [2] T. Askham and J. N. Kutz, Variable projection methods for an optimized dynamic mode decomposition, arXiv preprint arXiv:1704.02343, (2017).
  * [3] S. Bagheri, Koopman-mode decomposition of the cylinder wake, Journal of Fluid Mechanics, 726 (2013), pp. 596–623.
  * [4] S. Bagheri, Effects of weak noise on oscillating flows: linking quality factor, Floquet modes, and Koopman spectrum, Physics of Fluids, 26 (2014), p. 094104.
  * [5] Z. Bai, S. L. Brunton, B. W. Brunton, J. N. Kutz, E. Kaiser, A. Spohn, and B. R. Noack, 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.
  * [6] Z. Bai, T. Wimalajeewa, Z. Berger, G. Wang, M. Glauser, and P. K. Varshney, Low-dimensional approach for reconstruction of airfoil data via compressive sensing, AIAA Journal, 53 (2014), pp. 920–933.
  * [7] R. G. Baraniuk, Compressive sensing, IEEE Signal Processing Magazine, 24 (2007), pp. 118–120.
  * [8] P. Benner, S. Gugercin, and K. Willcox, A survey of projection-based model reduction methods for parametric dynamical systems, SIAM Review, 57 (2015), pp. 483–531.
  * [9] E. Berger, M. Sastuba, D. Vogt, B. Jung, and H. B. Amor, Estimation of perturbations in robotic behavior using dynamic mode decomposition, Journal of Advanced Robotics, 29 (2015), pp. 331–343.
  * [10] D. A. Bistrian and L. M. Navon, Randomized dynamic mode decomposition for non-intrusive reduced order modelling, International Journal for Numerical Methods in Engineering, 112 (2017), pp. 3–25.
  * [11] B. W. Brunton, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Sparse sensor placement optimization for classification, SIAM Journal on Applied Mathematics, 76 (2016), pp. 2099–2122.
  * [12] B. W. Brunton, L. A. Johnson, J. G. Ojemann, and J. N. Kutz, Extracting spatial–temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition, Journal of Neuroscience Methods, 258 (2016), pp. 1–15.
  * [13] S. L. Brunton, B. W. Brunton, J. L. Proctor, E. Kaiser, and J. N. Kutz, Chaos as an intermittently forced linear system, Nature Communications, 8 (2017).
  * [14] S. L. Brunton, S. T. Dawson, and C. W. Rowley, State-space model identification and feedback control of unsteady aerodynamic forces, J. Fluid & Struc., 50 (2014), pp. 253–270.
  * [15] S. L. Brunton and B. R. Noack, Closed-loop turbulence control: Progress and challenges, Applied Mechanics Reviews, 67 (2015), pp. 050801–1–050801–48.
  * [16] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences, 113 (2016), pp. 3932–3937.
  * [17] S. L. Brunton, J. L. Proctor, J. H. Tu, and J. N. Kutz, Compressed sensing and dynamic mode decomposition, Journal of Computational Dynamics, 2 (2015), pp. 165–191.
  * [18] S. L. Brunton, C. W. Rowley, and D. R. Williams, Reduced-order unsteady aerodynamic models at low Reynolds numbers, Journal of Fluid Mechanics, 724 (2013), pp. 203–233.
  * [19] E. J. Candès, J. Romberg, and T. Tao, Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information, IEEE Transactions on Information Theory, 52 (2006), pp. 489–509.
  * [20] K. K. Chen, J. H. Tu, and C. W. Rowley, Variants of dynamic mode decomposition: Boundary condition, Koopman, and Fourier analyses, Journal of Nonlinear Science, 22 (2012), pp. 887–915.
  * [21] Y. Chen, Y. Gu, and A. O. H. III, Sparse lms for system identification, IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP 2009), (2009), p. 10700918.
  * [22] T. Colonius 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.
  * [23] S. T. Dawson, M. S. Hemati, M. O. Williams, and C. W. Rowley, Characterizing and correcting for the effect of sensor noise in the dynamic mode decomposition, Experiments in Fluids, 57 (2016), pp. 1–19.
  * [24] D. L. Donoho, Compressed sensing, IEEE Transactions on Information Theory, 52 (2006), pp. 1289–1306.
  * [25] T. Duriez, S. L. Brunton, and B. R. Noack, Machine Learning Control: Taming Nonlinear Dynamics and Turbulence, Springer, 2016.
  * [26] N. B. Erichson, S. L. Brunton, and J. N. Kutz, Randomized dynamic mode decomposition, arXiv:1702.02912, (2017).
  * [27] J. Grosek and J. N. Kutz, Dynamic mode decomposition for real-time background/foreground separation in video, arXiv preprint arXiv:1404.7592, (2014).
  * [28] F. Gueniat, L. Mathelin, and L. Pastur, A dynamic mode decomposition approach for large and arbitrarily sampled systems, Physics of Fluids, 27 (2015), p. 025113.
  * [29] R. Heckel and H. Bölcskei, Identification of sparse linear operators, arXiv preprint arXiv:1209.5187v2, (2013).
  * [30] M. S. Hemati, S. T. Dawson, and C. W. Rowley, Parameter-varying aerodynamics models for aggressive pitching-response prediction, AIAA Journal, (2016), pp. 1–9.
  * [31] M. S. Hemati, C. W. Rowley, E. A. Deem, and L. N. Cattafesta, De-biasing the dynamic mode decomposition for applied Koopman spectral analysis, arXiv:1502.03854, (2015).
  * [32] P. J. Holmes, J. L. Lumley, G. Berkooz, and C. W. Rowley, Turbulence, coherent structures, dynamical systems and symmetry, Cambridge Monographs in Mechanics, Cambridge University Press, Cambridge, England, 2nd ed., 2012.
  * [33] M. R. Jovanović, P. J. Schmid, and J. W. Nichols, Sparsity-promoting dynamic mode decomposition, Physics of Fluids, 26 (2014), p. 024103.
  * [34] J. N. Juang and R. S. Pappa, An eigensystem realization algorithm for modal parameter identification and model reduction, Journal of Guidance, Control, and Dynamics, 8 (1985), pp. 620–627.
  * [35] E. Kaiser, J. N. Kutz, and S. L. Brunton, Data-driven discovery of Koopman eigenfunctions for control, arXiv preprint arXiv:1707.01146, (2017).
  * [36] E. Kaiser, M. Morzynski, G. Daviller, J. N. Kutz, B. W. Brunton, and S. L. Brunton, Sparsity enabled cluster reduced-order models for control, arXiv preprint arXiv:1701.00038, (2016).
  * [37] E. Kaiser, B. R. Noack, L. Cordier, A. Spohn, M. Segond, M. Abel, G. Daviller, J. Osth, S. Krajnovic, and R. K. Niven, Cluster-based reduced-order modelling of a mixing layer, J. Fluid Mech., 754 (2014), pp. 365–414.
  * [38] N. Kalouptsidis, G. Mileounis, B. Babadi, and V. Tarokh, Adaptive algorithms for sparse system identification, Signal Processing, (2011), pp. 1910–1919.
  * [39] B. O. Koopman, Hamiltonian systems and transformation in Hilbert space, Proceedings of the National Academy of Sciences, 17 (1931), pp. 315–318.
  * [40] Y. Kopsinis, K. Slavakis, and S. Theodoridis, Online sparse system identification and signal reconstruction using projections onto weighted $l_{1}$ balls, arXiv preprint arXiv:1004.3040, (2010).
  * [41] M. Korda and I. Mezić, Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control, arXiv:1611.03537, (2016).
  * [42] B. Kramer, P. Grover, P. Boufounos, M. Benosman, and S. Nabi, Sparse sensing and dmd based identification of flow regimes and bifurcations in complex flows, arXiv preprint arXiv:1510.02831, (2015).
  * [43] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2016.
  * [44] H. Leung, 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.
  * [45] L. Ljung, System identification, Wiley Online Library, 1999.
  * [46] J.-C. Loiseau and S. L. Brunton, Constrained sparse Galerkin regression, arXiv preprint arXiv:1611.03271, (2016).
  * [47] K. M., B. W. Brunton, J. N. Kutz, and S. L. Brunton, Data-driven sparse sensor placement, To appear in IEEE Control Systems Magazine, (2017).
  * [48] Z. Ma, S. Ahuja, and C. W. Rowley, Reduced order models for control of fluids using the eigensystem realization algorithm, Theor. and Comp. Fluid Dyn., 25 (2011), pp. 233–247.
  * [49] J. Mann and J. N. Kutz, Dynamic mode decomposition for financial trading strategies, Quantitative Finance, (2016), pp. 1–13.
  * [50] I. Mezić, Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dynamics, 41 (2005), pp. 309–325.
  * [51] I. Mezić, Analysis of fluid flows via spectral properties of the Koopman operator, Ann. Rev. Fluid Mech., 45 (2013), pp. 357–378.
  * [52] B. C. Moore, Principal component analysis in linear systems: Controllability, observability, and model reduction, IEEE Transactions on Automatic Control, AC-26 (1981), pp. 17–32.
  * [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.
  * [54] D. Needell and J. A. Tropp, CoSaMP: iterative signal recovery from incomplete and inaccurate samples, Communications of the ACM, 53 (2010), pp. 93–100.
  * [55] M. V. OL, A. Altman, J. D. Eldredge, D. J. Garmann, and Y. Lian, Résumé of the AIAA FDTC low Reynolds number discussion group’s canonical cases.  AIAA Paper 2010-1085, 48th Aerospace Sciences Meeting, January 2010.
  * [56] J. L. Proctor, S. L. Brunton, and J. N. Kutz, Dynamic mode decomposition with control, SIAM J. Applied Dynamical Systems, 15 (2016), pp. 142–161.
  * [57] J. L. Proctor, S. L. Brunton, and J. N. Kutz, Generalizing Koopman theory to allow for inputs and control, arXiv:1602.07647, (2016).
  * [58] J. L. Proctor and P. A. Eckhoff, Discovering dynamic patterns from infectious disease data using dynamic mode decomposition, International health, 7 (2015), pp. 139–145.
  * [59] C. W. Rowley, Model reduction for fluids using balanced proper orthogonal decomposition., International Journal of Bifurcation and Chaos, 15 (2005), pp. 997–1013.
  * [60] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson, Spectral analysis of nonlinear flows, Journal of Fluid Mechanics, 641 (2009), pp. 115–127.
  * [61] P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Mechanics, 656 (2010), pp. 5–28.
  * [62] P. J. Schmid, K. E. Meyer, and O. Pust, Dynamic mode decomposition and proper orthogonal decomposition of flow in a lid-driven cylindrical cavity.  PIV09-0186, 8th International Symposium on Particle Image Velocimetry, 2009.
  * [63] A. S. Sharma, I. Mezić, and B. J. McKeon, 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. Surana, 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. Surana and A. Banaszuk, Linear observer synthesis for nonlinear systems using koopman operator framework, IFAC-PapersOnLine, 49 (2016), pp. 716–723.
  * [66] K. Taira and T. Colonius, The immersed boundary method: a projection approach., Journal of Computational Physics, 225 (2007), pp. 2118–2137.
  * [67] N. Takeishi, Y. Kawahara, Y. Tabei, and T. Yairi, Bayesian dynamic mode decomposition, Twenty-Sixth International Joint Conference on Artificial Intelligence, (2017).
  * [68] N. Takeishi, Y. Kawahara, and T. Yairi, Subspace dynamic mode decomposition for stochastic Koopman analysis, arXiv:1705.04908, (2017).
  * [69] J. H. Tu, C. W. Rowley, J. N. Kutz, and J. K. Shang, Spectral analysis of fluid flows using sub-Nyquist-rate PIV data, Experiments in Fluids, 55 (2014), pp. 1–13.
  * [70] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz, On dynamic mode decomposition: Theory and applications, Journal of Computational Dynamics, 1 (2014), pp. 391–421.
  * [71] M. B. Wakin, B. M. Sanandaji, and T. L. Vincent, 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. Willcox and J. Peraire, Balanced model reduction via the proper orthogonal decomposition, AIAA Journal, 40 (2002), pp. 2323–2330.
  * [73] J. Y. Gu and S. Mei, $l_{0}$ norm constraint lms algorithm for sparse system identification, IEEE Signal Processing Letters, (2009), pp. 774–777.