###### Abstract

A data-driven analysis method known as dynamic mode decomposition (DMD)
approximates the linear Koopman operator on projected space. In the spirit of
Johnson-Lindenstrauss Lemma, we will use random projection to estimate the DMD
modes in reduced dimensional space. In practical applications, snapshots are
in high dimensional observable space and the DMD operator matrix is massive.
Hence, computing DMD with the full spectrum is infeasible, so our main
computational goal is estimating the eigenvalue and eigenvectors of the DMD
operator in a projected domain. We will generalize the current algorithm to
estimate a projected DMD operator. We focus on a powerful and a simple random
projection algorithm that will reduce the computational and storage cost.
While clearly, a random projection simplifies the algorithmic complexity of a
detailed optimal projection, as we will show, generally the results can be
excellent nonetheless, and quality understood through a well-developed theory
of random projections. We will demonstrate that modes can be calculated for a
low cost by the projected data with sufficient dimension.

###### keywords:

Koopman Operator, Dynamic Mode Decomposition(DMD), Johnson-Lindenstrauss
Lemma, Random Projection, Data-driven method.

xx 1 5 Received: date; Accepted: date; Published: date Randomized Projection
Learning Method for Dynamic Mode Decomposition Sudam Surasinghe 1 and Erik M.
Bollt 2 Sudam Surasinghe, Ioannis Kevrekidis and Erik M. Bollt

##  1 Introduction

Modeling real-world phenomena in physical sciences to engineering,
videography, economics is limited due to computational costs. Real-world
systems require dynamic nonlinear modeling and as Dynamic Mode Decomposition
(DMD)Schmid (2010); Rowley et al. (2009); Tu et al. (2014) is an emerging tool
in this area which can use the data directly rather than intermediary
differential equations. However, the data in real-world applications is
enormous. The DMD algorithm’s reliance on a Singular Value Decomposition(SVD)
appears to be a limiting factor due to storing and calculating the SVD of such
a matrix. One can notice that SVD calculation of snapshot matrices in big
projects would require the use of supercomputers or days of computation. In
order to allow processing on small scale computers and in a shorter time
frame, we propose to develop an algorithm based on a randomized
projectionDasgupta (2000) which is used to reduce the dimension of observable
space in DMD which we call rDMD. In order to utilize and carefully analyzed
the rDMD, we will use the Johnson-Lindenstrauss lemma. It is clear that a
random projection is simple as compared to a detailed optimal projection
method, but our analysis and examples demonstrate nonetheless the quality and
efficiency.

Strong theoretical support from Johnson-Lindenstrauss(JL) Lemma Johnson et al.
makes the random projection method reliable and has extensive utilization in
the field of data science. The JL lemma says that if data points lie in
sufficiently high dimensional space, then those data points may be projected
into a sufficiently low dimensional space while approximately preserving the
distance of the data points. Furthermore, the projection can be done just by a
random matrix which makes algorithms based on JL lemma both past and simple.
Hence, this tool is more powerful and adopted heavily in data science. JL
lemma based Random projection and SVD based projection can be used to project
$N$ dimensional data into lower dimension $L<<N$. Data matrix $X_{N\times M}$
can be projected by random projection into lower dimension ($L$) subspace as
$X_{L}:=RX$ where $R$ is a random matrix with unit length. Hence, the random
projection is very simple because it relies only on matrix multiplication.
Also computational complexity is $\mathcal{O}(MLN)$ while SVD has
computational complexity $\mathcal{O}(NM^{2})$ when $M<N$ Bingham and Mannila
(2001). We will use the random projection to project high dimensional snapshot
matrices into a manageable low dimensional space. In theoretical perspective,
the dimension of input and output spaces in the Koopman operator can be
reduced by the random projection method thus reducing the storage and
computational cost of the DMD algorithm.

Both DMD and rDMD are grounded in theory through the application of the
Koopman operator (Rowley.). They are numerical methods to estimate the linear
Koopman operator that identifies spatial and temporal patterns from a
dynamical system. Theoretical support of the Koopman operator theory makes the
these algorithms strong. Our new randomized DMD (rDMD) algorithm targets to
address issues that arise in SVD based existing DMD methods by reducing the
dimensionality of the matrix just by using matrix multiplication. rDMD can
achieve very accurate results with low-dimensional data embedded in high-
dimensional observable space.

This paper will summarize the Koopman operator theory, existing DMD
algorithms, and random projection theory in section (2). Then we will discuss
our proposed randomized DMD algorithm in section (3), and finally, in section
(4), we will provide examples that support our approach.

##  2 Dynamic Mode Decomposition and Background Theory

The focus of this paper is to approximate the eigenpairs of the Koopman
operator based on the randomized dynamic mode decomposition method. We will
first review the underlining theory about the Koopman operator.

###  2.1 Koopman Operator

Consider a discrete-time dynamical system

$\displaystyle x_{n+1}=S(x_{n}),$ (1)

where $S:\mathcal{M}\to\mathcal{M}$ and $\mathcal{M}$ is a finite dimensional
manifold. (If we have a differential equation or continuous time dynamical
system, flow map can be considered.) The variable $x$ is often recognized as a
state variable and $\mathcal{M}$ as phase space. The associated Koopman
operator is described as the evaluation of observable functions (Fig. 1)
$\psi:\mathcal{M}\to\mathbb{R}$ in function space $\mathcal{F}$. Instead of
analysing the individual trajectories in phase space, the Koopman operator
operates on the observations Koopman (1931); Bollt (2021); Tu et al. (2014);
Rowley et al. (2009).

Figure 1: This figure shows the behavior of Koopman operator $\mathcal{K}$ in
observable space $\mathcal{F}$ associated with a dynamical system $S$. The
Koopman operator evaluate the observable $\psi$ at downstream or future
$x^{\prime}=x_{n+1}=s(x_{n})$.

[Koopman operator Koopman (1931)] The Koopman operator $\mathcal{K}$ for a map
$S$ is defined as following composition,

$\displaystyle\mathcal{K}$ $\displaystyle:\mathcal{F}\rightarrow\mathcal{F}$
$\displaystyle\psi\mapsto\mathcal{K}[\psi]$ $\displaystyle=\psi\circ S$ (2)

on the function space $\mathcal{F}$. It is straight forward to prove Koopman
(1931),

$\displaystyle\mathcal{K}[a\psi_{1}+b\psi_{2}]=a(\psi_{1}\circ
S)+b(\psi_{2}\circ S)=a\mathcal{K}[\psi_{1}]+b\mathcal{K}[\psi_{2}]$ (3)

for $\psi_{1},\psi_{2}\in\mathcal{F}$ and $a,b\in\mathbb{C}$, and therefore
the Koopman operator is linear on $\mathcal{F}$. This is interesting and
important property of the operator because the associated map $S$ most
probably will be non-linear. Even though the operator is associated with a map
that evolves in finite dimensional space, $\mathcal{F}$ the function space
which the operator acts on could possibly be an infinite dimensional. This is
the trade-off between cost for the linearity Bollt (2021).

Spectral analysis of the Koopman operator can be used to decompose the
dynamics which becomes the key success in the DMD. Assuming the spectrum of
Koopman operator $\mathcal{K}$ is given by

$\displaystyle\mathcal{K}\psi_{i}(x)=\lambda_{i}\psi_{i}(x)~{}~{}~{}i=1,\ 2,\
3,\ \dots$ (4)

then vector-valued observables
$\boldsymbol{g}:\mathcal{M}\to\mathbb{R}^{N}$(or $\mathbb{C}^{N}$) can be
represented by

$\displaystyle\boldsymbol{g(x)}=\sum_{i=1}^{\infty}\psi_{i}(\boldsymbol{x})\boldsymbol{\phi_{i}},$
(5)

where $\boldsymbol{\phi_{i}}\in\mathbb{R}^{N}$(or $\mathbb{C}^{N}$) are the
vector coefficients of the expansion and called “Koopman modes”(here we
assumed that components of $\boldsymbol{g}$ lie within the span of the
eigenfunctions of $\mathcal{K}$). Note that the observable value at time $n+1$
is given by

$\displaystyle\boldsymbol{g(x_{n+1})}=\sum_{i=1}^{\infty}\lambda_{i}^{n}\psi_{i}(\boldsymbol{x_{0}})\boldsymbol{\phi_{i}}.$
(6)

This decomposition can be used to separate the spacial and time components of
the dynamical system and can be used to isolate the specific dynamics.

###  2.2 Dynamic Mode Decomposition

The dynamic mode decomposition is a data-driven method to estimate the Koopman
modes from numerical or experimental data Rowley et al. (2009). Suppose
dynamics are govern by eq. (1) for any state $\boldsymbol{x}$ and vector
valued measurements are given by observable
$\boldsymbol{g(x)}\in\mathbb{R}^{N}$. For a given set of data
$X=[\boldsymbol{g(x_{0})}\ \boldsymbol{g(x_{1})}\ \dots\
\boldsymbol{g(x_{M-1})}],\ Y=[\boldsymbol{y_{0}}\ \boldsymbol{y_{1}}\ \dots\
\boldsymbol{y_{M-1}}]$ where $\boldsymbol{y_{i}}=\boldsymbol{g(s(x_{i}))}$,
the Koopman modes and eigenvalues of the Koopman operator can be estimated
through solving the least-squares problem

$\displaystyle\mathbb{K}=\argmin_{K}||KX-Y||_{F}^{2}=\argmin_{K}\sum_{i=0}^{M-1}||K\boldsymbol{g(x_{i})}-\boldsymbol{y_{i}}||_{2}^{2}$
(7)

and $\mathbb{K}=YX^{\dagger}$ (here $X^{\dagger}$ is the pseudo-inverse of
$X$) is defined as the "Exact DMD" operator Tu et al. (2014). The eigenvalue
($\hat{\lambda}$) of $\mathbb{K}$ is an approximation of an eigenvalue
($\lambda$) of $\mathcal{K}$; the corresponding right
eigenvector($\hat{\phi}$) is called the DMD mode and approximates the Koopman
mode ($\phi$). Then the observable value $\boldsymbol{g(x(t))}$ at time $t$
can be modeled as

$\displaystyle\boldsymbol{g(x(t))}=\sum_{i=1}^{r}\psi_{i}(\boldsymbol{x_{0}})\hat{\boldsymbol{\phi}}_{i}\hat{\lambda}_{i}^{t}$
(8)

where $r$ is the number of selected DMD modes and demonstrates the finite
dimensional approximation for vector-valued observable $\boldsymbol{g}$ under
the Koopman operator. Based on this decomposition, data matrices can be
expressed as

$\displaystyle X_{N\times M}$ $\displaystyle=\Phi_{N\times r}T_{r\times M}$
(9) $\displaystyle Y_{N\times M}$ $\displaystyle=\Phi_{N\times
r}\Lambda_{r\times r}T_{r\times M}$

where
$\Phi=[\psi_{1}(x_{0})\boldsymbol{\phi}_{1}~{}~{}\psi_{2}(x_{0})\boldsymbol{\phi}_{2}~{}\dots~{}\psi_{r}(x_{0})\boldsymbol{\phi}_{r}]$,
$T$ is a Vandermonde matrix with $T_{ij}=\lambda_{i}^{j-1}$ for $i=1,2,\dots
r$, $j=1,2,\dots,M$ and
$\Lambda=diag\\{\lambda_{1},~{}\lambda_{2},\dots,\lambda_{r}\\}$. Note that
with the above decomposition $\mathbb{K}=YX^{\dagger}=\Phi\Lambda
TT^{\dagger}\Phi^{\dagger}$. We will suppose $\mathbb{K}$ has distinct
eigenvalues $\lambda_{i}$, columns of $X$ are linearly independent and $r\leq
M$. In practical applications, we are expected to fully understand the data
set by relatively few ($r<<M$) modes. This can be considered as one of the
dimension reduction steps of the algorithm. Additionally, dimension of columns
of the data matrix need to be reduced.

In practice, the columns of data matrix $X$ (and $Y$) are constructed by the
snapshot matrices of spatial observable data. More often, those snapshots lie
in high dimension space $\mathbb{R}^{N}$ ($N>>1$ and roughly
$\mathcal{O}(10^{15})$ to $\mathcal{O}(10^{10})$), but the number of snapshots
or time steps ($M$) are small and often it is $\mathcal{O}(10^{3})$ to
$\mathcal{O}(10^{1})$ Chen et al. (2012). Hence, computing the spectrum of
matrix $\mathbb{K}$ is infeasible even though most of the eigenvalues will be
zero. Therefore we can project our data matrices $X,Y$ into a low dimensional
space $R^{L}$ with $r\leq L\leq M<<N$, therefore need to estimate the spectrum
of $\mathbb{K}$ based on the computation on projected space. Our proposed rDMD
method is focused on this dimension reduction step.

In the next section (Sec. (3)) we will discuss more details about the
calculation. Note that current methods are based on the singular value
decomposition of the data matrix $X$ to construct a projection and our
proposed algorithm is based on the random projection method to project data
into a low dimensional space.

###  2.3 Random projection

The random projection method is based on the Johnson-Lindenstrauss lemma which
is applied by many data analysis methods.

[Johnson-Lindenstrauss LemmaJohnson et al. ] For any $0<\epsilon<1$ and any
integer $M>1,$ let $L$ be a positive integer such that $L\geq L_{0}$ with
$L_{0}=\frac{C\ln M}{\epsilon^{2}},$ where $C$ is a suitable constant
($C\approx 8$ in practice,$C=2$ is good enough). Then for any set $X$ of $M$
data points in $\mathbb{R}^{N}$, there exists a map
$f:\mathbb{R}^{N}\to\mathbb{R}^{L}$ such that for all $x_{1},x_{2}\in X$,

$(1-\epsilon)||x_{1}-x_{2}||^{2}\leq||f(x_{1})-f(x_{2})||^{2}\leq(1+\epsilon)||x_{1}-x_{2}||^{2}.$

[Random Projection Dasgupta (2000)] For any $0<\epsilon,\delta<\frac{1}{2}$
and positive integer $N$, there exists a random matrix of $B$ of size $L\times
N$ such that for $L\geq L_{0}$ with
$L_{0}=\frac{C\ln(1/\delta)}{\epsilon^{2}}$. and for any unit-length vector
$x\in R^{N}$

$Pr\\{|||Bx||^{2}-1|>\epsilon\\}\leq\delta$

or

$Pr\\{|||Bx||^{2}-1|>\epsilon\\}\leq e^{-CL\epsilon^{2}}$

A low rank approximation for both $X,Y$ can be found using the random
projection method. Notice that both these matrices have $M$ points from $N$
dimensional observable space and therefore we can use random projection matrix
$B$ of size $L\times N$ with the $L\geq\frac{C\ln N}{\epsilon^{2}}$ which
provides $\epsilon-$ isometry to $\mathbb{R}^{N}$. (See Fig. 2 for details).

##  3 Randomized Dynamic Mode Decomposition

In this section we will generalize currently used DMD algorithms and then we
will discuss our proposed randomized DMD algorithm.

###  3.1 DMD on projected space

As we mentioned in section 2.2, computational and storage cost of DMD can be
reduced by projecting data into a low dimensional observable space. Let
$P\in\mathbb{R}^{L\times N}$ be any rank $L$ projection matrix, then dimension
of data matrices $X,\ Y\in\mathbb{R}^{N\times M}$ can be reduced to $L\times
M$ by the projection $X_{L}=PX,\ Y_{L}=PY$. The DMD operator on projected
space (see Fig. (2)) is given by,

$\displaystyle\hat{\mathbb{K}}=\argmin_{K\in\mathbb{R}^{L\times
L}}||KX_{L}-Y_{L}||_{F}^{2}=\argmin_{K\in\mathbb{R}^{L\times L}}||KPX-
PY||_{F}^{2}$ (10)

and $\hat{\mathbb{K}}=PY(PX)^{\dagger}$. Therefor

$\displaystyle\hat{\mathbb{K}}=PY(PX)^{\dagger}=PYX^{\dagger}P^{\dagger}=P\mathbb{K}P^{\dagger}$
(11)

where $\mathbb{K}=YX^{\dagger}$ is the DMD operator on original space.

Some eigenpairs $(\lambda,\phi)$ of $\mathbb{K}$ can be obtain by
$(\lambda_{L},\phi_{L})$ of projected DMD $\hat{\mathbb{K}}$ with
$\lambda=\lambda_{L}$ and $\phi=P^{\dagger}\phi_{L}$.

###### Proof.

Let $(\lambda_{L},\phi_{L})$ be an eigenpair of $\hat{\mathbb{K}}$. Then
$\hat{\mathbb{K}}\phi_{L}=\lambda_{L}\phi_{L}$ and by eq.(11),
$P\mathbb{K}P^{\dagger}\phi_{L}=\lambda_{L}\phi_{L}$. Now let
$P^{\dagger}\phi_{L}=\phi$, then $\phi_{L}=P\phi$ because $PP^{\dagger}=I$.
Hence $P\mathbb{K}P^{\dagger}\phi_{L}=P\mathbb{K}\phi=\lambda_{L}P\phi$ and
$P(\mathbb{K}\phi-\lambda_{L}\phi)=0$. Since
$\mathbb{K}\phi-\lambda_{L}\phi=0$ is a solution to the above equation,
$\lambda_{L}$ is an eigenvalue and the corresponding eigenvector is
$\phi=P^{\dagger}\phi_{L}$ of $\mathbb{K}$. ∎

In other words, we can lift up the dimension of eigenvectors in projected
space by $P^{\dagger}$ to obtain an eigenvector in original data space.
However to avoid the direct calculation of the pseudo-inverse of projection
matrix , we can calculate the eigenvector in output space $Y_{L}$ of the DMD
operator and lift up the vector into the original output space $Y$. We can
easily show that $\hat{\phi}=Y(PX)^{\dagger}\phi_{L}$ is an eigenvector of
$\mathbb{K}$ for corresponding non-zero eigenvalues.

$\displaystyle\mathbb{K}\hat{\phi}$
$\displaystyle=\mathbb{K}Y(PX)^{\dagger}\phi_{L}=\mathbb{K}YX^{\dagger}P^{\dagger}\phi_{L}$
(12) $\displaystyle=\mathbb{K}\mathbb{K}P^{\dagger}\phi_{L}$
$\displaystyle=\mathbb{K}\lambda P^{\dagger}\phi_{L}$ $\displaystyle=\lambda
YX^{\dagger}P^{\dagger}\phi_{L}$ $\displaystyle=\lambda
Y(PX)^{\dagger}\phi_{L}$ $\displaystyle=\lambda\hat{\phi}$

Also notice, $P\hat{\phi}=PY(PX)^{\dagger}\phi_{L}=\hat{\mathbb{K}}\phi_{L}$
and therefor $\hat{\phi}$ estimate the eigenvector on output space $Y$.
Detailed view of this lifting operator is shown by Fig. (3). It provides the
relationship of the lifting operator with the DMD operator acted on any
general observable vector $z=\boldsymbol{g(x(t))}\in\mathbb{R}^{N}$.

Next, the focus move to the spatial-temporal decomposition of projected data
matrices by spectrum of the DMD operator. Note that observable value
$\boldsymbol{g(x(t))}$ at time $t$ can be modeled as
$\boldsymbol{g(x(t))}=\sum_{i=1}^{r}\psi_{i}(\boldsymbol{x_{0}})P^{\dagger}\boldsymbol{(\phi_{L})}_{i}\hat{\lambda}_{i}^{t}$
and similar to the eq. (9), data can be decomposed as

$\displaystyle X_{N\times M}$ $\displaystyle=P_{N\times
L}^{\dagger}\tilde{\Phi}_{L\times r}T_{r\times M}$ (13) $\displaystyle
Y_{N\times M}$ $\displaystyle=P_{N\times L}^{\dagger}\tilde{\Phi}_{L\times
r}\Lambda_{r\times r}T_{r\times M}.$

This decomposition leads to
$\mathbb{K}=YX^{\dagger}=P^{\dagger}\tilde{\Phi}\Lambda
TT^{\dagger}\tilde{\Phi}^{\dagger}P$ and if $r\leq L$ all the non-zero
eigenvalues and corresponding eigenvectors of $\mathbb{K}$ can be constructed
by the projected DMD operator. Further, eq. (13) can be use to isolate the
spatial profile of interesting dynamical compotes such as attractors, periodic
behaviors, etc.

Based on the choice of the projection matrix we will have alternative ways to
estimated the spectrum of the DMD operator.

######  Remark 3.1 (Projection by SVD).

Commonly used projection matrix is based on SVD of the input matrix $X=U\Sigma
V^{*}$ and projection matrix is chosen to be $P=U^{*}$, here ∗ represents the
conjugate transpose of a matrix. Using eq. 11 and SVD of X, the operator on
projected space can be formulated as $\hat{\mathbb{K}}=U^{*}YV\Sigma^{-1}$.

######  Remark 3.2 (Standard DMD and Exact DMD).

Let eigenpair of a SVD based $\hat{\mathbb{K}}=U^{*}YV\Sigma^{-1}$ be given by
$(\lambda,\phi_{L})$. In standard DMD (Ref. Schmid Paper& Tu paper) use the
eigenvector $P^{\dagger}\phi_{L}=U\phi_{L}$ to estimate eigenvectors of
$\mathbb{K}$. On the other hand, in exact DMD(ref. Tu paper) this eigenvector
is estimated by $Y(PX)^{\dagger}\phi_{L}=YV\Sigma^{-1}\phi_{L}$.

######  Remark 3.3.

QR decomposition based projection methods on both input and output data $[X\
Y]$ can be used Tu et al. (2014).

In this paper we are proposing a simple random projection based method to
estimate the spectrum of the DMD operator.

Figure 2: DMD operator on projected space. This figure shows the relationship
between DMD operator on original space and DMD on projected space. The
operator $\mathbb{\hat{K}}$ on projected space is defined in Eq. (10) and can
be calculated by Eq. (10).

###  3.2 Randomized Methods

Our suggested randomized Dynamic Mode Decomposition(rDMD) is based on the
random projection applied to the theory of DMD on projected space. We can
reduce the dimension of data matrix $X,\ Y$ in DMD by using a random
projection matrix $R_{L\times N}$. In other words we will construct a
projection matrix $P$ discussed in sec. 3.1 as a random matrix $R$ whose
columns have unit length and entries are selected independently and
identically from a probability distribution. Therefore, the rDMD matrix on the
projected space is given by $\hat{\mathbb{K}}=RY(RX)^{\dagger}$, and if an
eigenpair of $\hat{\mathbb{K}}$ is given by $(\lambda,\phi_{L})$, then the
eigenpair of $\mathbb{K}$ is given $(\lambda,Y(RX)^{\dagger}\phi_{L})$.
Algorithm 1 represents the major steps needed to estimate the eigenvalues and
corresponding eigenvectors of the DMD operator with random the projection
method.

The calculation of the projection matrix of a standard or exact DMD algorithm
based on SVD of snapshot matrix $X$ is needed to store a full high resolution
data matrix which leads to memory issues. Our proposed rDMD algorithm can
avoid these storage issues, because low dimensional matrices $X_{L},\ Y_{L}$
obtained by matrix multiplications only need to store one row and one column
of each matrix at a time. Additionally, this algorithm reduces the
computational cost since we only need to calculate pseudo-inverse of
comparatively lower dimensional matrix. Choice of the distribution of $R$ can
further reduce the computational cost Achlioptas (2003).

Data: $X,\ Y\in\mathbb{R}^{N\times M}$

$\epsilon$ $L_{0}=\frac{C\ln M}{\epsilon^{2}}$;

Choose $L$ such that $L\geq L_{0}$;

Construct a random matrix $R=\frac{1}{\sqrt{L}}(r_{ij})\in\mathbb{R}^{L\times
N}$ such that $r_{ij}~{}N(0,1)$;

Calculate $X_{L}:=RX,\ Y_{L}:=RY$ ;

Calculate $\hat{\mathbb{K}}=Y_{L}X_{L}^{\dagger}$;

[$\lambda$ $\ \Phi_{L}$]=eigs($\hat{\mathbb{K}}$);

Result: $diag(\Lambda)$, $YX_{L}^{\dagger}\Phi_{L}$

Algorithm 1 Randomized DMD(rDMD)

One time step forecasting error for any given snapshot by using rDMD algorithm
can be bounded by using the JL theory.

[Error Bound] Let $z=\boldsymbol{g(x(t))},\
z^{\prime}=\boldsymbol{g(x(t+1))}\in\mathbb{R}^{N}$. Error bound of estimating
$z^{\prime}$ by using the rDMD as $\hat{z}^{\prime}=YX^{\dagger}R^{\dagger}Rz$
is given by

$\displaystyle
E(z^{\prime};L):=||z^{\prime}-\hat{z}^{\prime}||\leq\frac{||Rz^{\prime}-\hat{\mathbb{K}}Rz||}{1-\epsilon}:=UB$
(14)

with at least the probability of $\mathcal{O}(1/M^{2})$ for any $0<\epsilon<1$
with $L>\frac{C\log(M)}{\epsilon^{2}}$.

######  Proof 3.4.

Since $\hat{\mathbb{K}}=RYX^{\dagger}R^{\dagger}$, the rDMD acts on the
projected vector which can be rearranged as
$\hat{\mathbb{K}}Rz=R\hat{z}^{\prime}$. Therefore,

$\displaystyle||Rz^{\prime}-\hat{\mathbb{K}}Rz||=||Rz^{\prime}-R\hat{z}^{\prime}||.$

Now we can apply the JL theory to attain the desired error bound.

$\displaystyle(1-\epsilon)||z^{\prime}-\hat{z}^{\prime}||\leq||Rz^{\prime}-R\hat{z}^{\prime}||=||Rz^{\prime}-\hat{\mathbb{K}}Rz||.$

Hence
$||z^{\prime}-\hat{z}^{\prime}||\leq\frac{||Rz^{\prime}-\hat{\mathbb{K}}Rz||}{1-\epsilon}$.

Figure 3:  The figure shows the projecting operator $P$ and DMD related
lifting operator $Y(PX)^{\dagger}=\mathbb{K}P^{\dagger}$ which should be used
in DMD algorithms. Instead of using $P^{\dagger}$ as the lifting operator,
$Y(PX)^{\dagger}$ can be used for efficient calculations. Moreover, notice
that $\hat{z}^{\prime}_{L}:=P\hat{z}^{\prime}=\hat{\mathbb{K}}z_{L}$.

##  4 Results and Discussion

In this section we will demonstrate the theory of rDMD with a few examples.
The first two examples consider the computation for known dynamics and
demonstrate the error analysis. The final example will demonstrate application
in the field of oceanography and isolate the interesting features by rDMD and
compere with the resulting modes with the exact DMD results.

###  4.1 Logistic Map

We will first consider a dataset of 300 snapshots from a logistic map,

$\displaystyle x_{n+1}=ax_{n}(1-x_{n})$

with $a=3.56994$. In this case all the initial conditions will converge to a
period-256 orbit. Therefore rank of the snapshot matrix with relatively high
samples should be 256. We forecast the data by using the rDMD method and then
analyzed the error of the prediction and compared it with the theoretical
upper bound. With $N=5000$ initial conditions and $M=300$ samples, the
dimension $L$ of the projecting space can be chosen as
$L\geq\frac{C\ln(300)}{\epsilon^{2}}\approx\frac{34.22}{\epsilon^{2}}$ when
$C=6$. rDMD with projection into 50 dimensional space can accurately forecast
the time series data. (Fig. (4) shows the original vs predicted data for one
trajectory.) Furthermore, fig. (5) demonstrates the bound of the error of
forecast explained in eq. (14) and how error relates with the distortion
parameter $\epsilon$ (fig. (5) (a)) and dimension of the projected space
((fig. (5) (b))). Since the rank of the snapshot matrix is $256$, any $L\geq
256$ will perform very accurately. This example validates the error bound we
discussed in Eq. (14) and error of the prediction depends on the error
exhibited by the projected DMD operator and the distortion parameter
($\epsilon$ or the projected dimension) from the JL theory.

Figure 4:  (Left) shows the predicted data using the rDMD algorithm with
projected dimension $L=50$ compared to the original data from logistic map
$x_{n+1}=3.56994x_{n}(1-x_{n})$ for initial condition $x(0)=0.4967$. Further,
(Right) figure shows the estimated error and theoretical upper bounds (Eq. 14)
for the some projected dimension $L$ and this example validates the
theoretical bound.

(a)

(b)

Figure 5:  (a),(b) shows the prediction error of logistic map
$x_{n+1}=3.56994x_{n}(1-x_{n})$ by rDMD and its theoretical upper bounds(Eq.
14). Figure (a) represents the error with respect to the distortion $\epsilon$
and (b) shows the error with dimenstion of the projected space that will
guarantee the bound for this example.

###  4.2 Toy Example: Demonstrates the Variable Separation and isolating
dynamics.

To demonstrate the variable separation and to isolate the spatial structures
based on the time dynamics, we consider a toy example(motivated by the example
in Kutz et al. (2016)) ,

$\displaystyle
z(x,t)=\sum_{j=1}^{20}j\sech(0.1x+j)e^{i\gamma_{j}t}=\sum_{j=1}^{20}\Phi_{j}(x)T_{j}(t)$
(15)

where $\gamma_{j}$’s are constants, and let $\Phi_{j}(x)=j\sech(0.1x+j)$ and
$T_{j}(t)=e^{i\gamma_{j}t}$.

Figure 6: original dataset constructed by eq. (15).(Left) shows the
$\real(z(x,t))$ plot for all $x,t$ values at $20000\times 5000$ grid points.
(middle) represents the time series plot for two initial conditions. (Right)
provides the snapshots of few different time points. Our goal is separate and
isolate spatial variables $\Phi_{j}(x)=j\sech(0.1x+j)$ and $T_{j}(t)=e^{ijt}$
from the given data constructed by $z(x,t)$.

Comparing this Eq. (15) with decomposition Eq. (8), the rDMD algorithm is
expected to isolate 20 periodic modes by rDMD algorithm. The data set
(snapshot matrix) for this problem is constructed by $N=20000$ spatial grid
points and $M=5001$ temporal grid points with $\gamma_{j}=j$ (see fig. (6)).
As discussed in the previous section, if $L\geq 20$ then those expected modes
can be isolated and there exist eigenvalues $\lambda_{j}$ of rDMD operator
such that

$\displaystyle\omega_{j}=\ln(\lambda_{j})=ji$

for $j=1,\ 2,\ \dots,\ 20$ (see fig. (7)). Furthermore, we are expecting
corresponding rDMD modes equal to spatial variables of the model such that

$\displaystyle b_{j}(x_{0})\phi_{j}(x)=j\sech(0.1x+j)=\Phi_{j}(x).$

As expected, we noticed that calculated modes have negligible error when the
dimension of projected space $L\geq r=20$. Figure 7 shows the absolute error
of eigenvalues and DMD modes. (All the modes behave similarly and here we
present mode 10 for demonstration purpose.) by SVD based exact DMD method and
random projection based rDMD method. Notice that errors of both methods are
less than $10^{-10}$ when $L\geq r=20$.

Figure 7:  (Left) shows the absolute error for estimated eigenvalues from rDMD and exact DMD when the dimension of projected space L=20. (Right) shows the absolute error for estimated 10th DMD mode by rDMD and exact DMD methods. In this case, both methods have very accurate results and error is less than $10^{-10}$.  Method | Projected by | Computational Time($\mathrm{s}$)
---|---|---
Exact DMD | SVD | 521.09
rDMD | Random Projection | 2.35
Table 1: Computational cost for the SVD based exact DMD and random projection
based rDMD method for the data simulated by Eq. (15). Computational cost of
SVD for high dimensional snapshot matrix is relatively larger than random
projection.

Further, we examine the case when the projected dimension $L=17<r=20$ and
compared the results of rDMD with the exact DMD. We can notice that both
methods demonstrate similar errors and rDMD is almost good as the SVD
projection based Exact DMD(See Fig. (8) and (9)). When the number of actual
modes($r$) is larger than the dimension of the projected space ($L$), the
projected DMD operator only estimates the $L$ number of modes, leads to both
truncation errors and error for eigenpair estimation based on projected DMD
operator. The $L<r$ case can be modeled as,

$\displaystyle
z(x,t)=\sum_{j=1}^{L}\hat{b}_{j}\hat{\phi}_{j}(x)e^{\hat{\omega}_{j}it}+E_{L+1}$
(16)

Where $E_{L+1}=\sum_{L=j+1}^{m}b_{j}\phi_{j}(x)e^{\omega_{j}it}$ is the
truncated error that also affects the estimation process of eigenpairs.
Therefore, if $L<r$ then there exists an error in eigenvalues and eigenvectors
calculated by any method based on the projected DMD. However, this example
demonstrates that rDMD can provide the results as good as the SVD projection
based method with very low computational cost(See Table. (1)).

Figure 8:  (Left) compares the eigenvalues $\omega_{j}=\ln\lambda_{j}$
calculated from rDMD (random projection(RP) with L=17) and exact DMD (SVD
projection with L=17) methods with the expected true values $\gamma_{j}=ji$.
Here $L=17<20$ is the dimension of projected space. (Right) shows the absolute
error for estimated eigenvalues from rDMD and exact DMD.

Figure 9:  (Left) compares the modes $b_{10}\phi_{10}(x)$ calculated from rDMD
(random projection(RP) with L=17) and exact DMD (SVD projection with L=17)
methods with the expected true values $\Phi_{10}(x)=10\sech(0.1x+10)$. (Right)
shows the absolute error for estimated values from rDMD and exact DMD.

###  4.3 Gulf of Mexico

In this example we will consider the data from HYbrid Coordinate Ocean Model
(HYCOM) hyc (2021) which simulates the ocean data around the Gulf of Mexico.
We used hourly surface velocity component ($u,v$) with $1/25^{0}$ spatial
resolution ($N=541\times 347$ grid points) data for 10 days ($240$
$\mathrm{h}$ and $M=239$). Understanding the dynamics from the oceanographic
data is an interesting application of DMD because those dynamics can be
decomposed by tidal constituents . Hence, we are expected to isolate the
dynamics associated with the tidal period ; in other words, the final DMD mode
selection is based on the period $P_{i}=2\pi/\imaginary(\ln(\lambda_{i}))$ of
the modes(see table (2)). We constructed the snapshot matrix

$\displaystyle X=\begin{bmatrix}u\\\ v\end{bmatrix}$ (17)

by stacking the snapshots of velocity components ($u,v$) in each column to
perform the DMD analysis.

Mode | Period($\mathrm{h}$) | Associated Feature
---|---|---
| DMD | rDMD |
1 | $\infty$ | $\infty$ | Gulf stream around the GOM.(see fig. (11))
2 | $12.47$ | $12.47$ | Semi-diurnal tidal constituents.(see fig. (12))
3 | $23.85$ | $24.56$ | Diurnal tidal constituents.(see fig. (13))
4 | $6.08$ | $6.07$ | 2nd harmonic to semi-diurnal tidal constituents. (see fig. (13))
5 | $4.16$ | $4.17$ | 3rd harmonic to semi-diurnal tidal constituents. (see fig. (13))
Table 2: DMD modes for Gulf of Mexico data set. Modes are selected based on
the association to the tidal periods.

(a)

(b)

Figure 10:  Eigenvalues $\lambda_{i}$ calculated from exact DMD and rDMD
methods. (a) Full spectrum of the two methods with projected space dimension
$L=239$ and (b) shows the first 5 modes. The mode selection is based on the
comparison of the tidal periods with period of the DMD modes.

Figure (10) shows that most of the eigenvalues calculated from SVD based exact
DMD and random projection based rDMD are agree. Furthermore, eigenvalues that
isolated the specific dynamics are almost equal. Additionally, figure
(11)-(13) shows the spacial profile of those modes from exact DMD and rDMD
methods. Also, each mode clearly isolated the interesting oceanographic
features(see table (2)) and both methods provide almost the same spacial
structures(see fig. (11)-(13)) as expected.

(a)

(b)

Figure 11: This figure compares the (a) DMD and (b) rDMD background mode
identified by data from the Gulf of Mexico (GOM). This background mode
captures the ocean current passing through the GOM.

(a)

(b)

Figure 12: This figure compere the (a) DMD and (b) rDMD mode associated with
M2 tidal frequency. This mode capture the "red tides".

Notice that the dimension of the snapshot matrix is $375454\times 239$ and SVD
calculation of this matrix is more costly for both computation and storage. On
the other hand, random projection performs only by matrix multiplication which
can be done at relatively low cost. Hence, we can achieve almost the same
results by using random projection method at a relatively much lower
computational and storage cost.

Figure 13: First row represent the exact DMD modes 3, 4 and 5 (left to right)
and second row shows the rDMD modes 3, 4 and 5 ((left to right)). Mode 3 is a
diurnal mode with period $23.85$ $\mathrm{h}\mathrm{o}\mathrm{u}\mathrm{r}$
for exact DMD case and $24.56$$\mathrm{h}\mathrm{o}\mathrm{u}\mathrm{r}$ for
rDMD case. Mode 4 and 5 are associated with 2nd and 3rd harmonic of a semi-
diurnal tidal constituents respectively.

##  5 Conclusion

We have demonstrated that our rDMD can achieve very accurate results with low-
dimensional data embedded in high-dimensional observable space. Recent
analytic technology from the concepts of high-dimensional geometry of data,
and concentration of measure have born out that perhaps surprising if not
initially intuitively that even random projection methods can be quite
powerful and capable. Here, in the setting of DMD methods approximating and
projecting the action of a Koopman operator, we show that randomized
projection can be developed and analyzed rigorously by the Johnson-
Lindenstrass theorem formalism, this showing a powerful and simple approach.
We provided a theoretical framework and experimental results to address those
issues raised from SVD based methods by introducing our new rDMD algorithm.
The theoretical framework is based on generalizing the SVD based concept as a
projection of high dimensional data into a low dimensional space. We proved
that eigenpairs of DMD in original space can be estimated by using any rank
$L$ projection matrix $P$. Being able to estimate eigenpairs allowed us to use
the powerful and simple Johnson-Lindenstrauss lemma and the random projection
method allowing us to project data with matrix multiplication. Therefore our
proposed random projection-based DMD(rDMD) can estimate eigenpairs of the DMD
operator with low storage and computational cost. Further, the error of the
estimation can be controlled by choosing the dimension of the projected space
and we demonstrated this error bound through the "logistic map" example.

DMD promises the separation of the spatial and time variables from data.
Hence, we experimentally demonstrated how well the rDMD algorithm performed
this task by a toy example. Notice that the number of those isolated modes
($m$) are relatively (to spatial and temporal resolution) low in practical
applications. If $m<<M$, then the rank of the data matrix is much lower, and
those eigenvalues and vectors of interest can be estimated accurately by
projecting data into the much lower dimensional space $L\geq m$. The SVD
projection-based exact DMD method still needs to calculate the SVD of a high
dimensional(roughly $10^{10}\times 10^{3}$ ) data matrix while rDMD only
requires to multiply the data matrix by a much lower-dimensional projection
matrix. Furthermore, we noticed that both exact and random DMD methods are
experiencing similar errors. However random projection is much faster and
needs less space for the calculations. We also demonstrate that practical
applications also provide similar results by using oceanographic data from the
Gulf of Mexico.

Since the size of the DMD matrix is enormous in those applications (this could
be roughly $10^{10}\times 10^{10}$), the eigenpairs of it must be estimated by
projecting data into low dimensional space. Estimating eigenvalues and
eigenvectors of a DMD operator using a high dimensional snapshot data matrix
(in applications this could be $10^{10}\times 10^{3}$) with existing SVD based
methods is expensive. The computational efficiency of the rDMD led to a new
path of current Koopman analysis. It allows using more observable variables in
the data matrix without need of much extra computational power. Hence, state
variables and more non-linear terms of them can be used in analysis with low
cost to improve the Koopman modes. JL theory can be adopted further into the
field of numerical methods of Koopman theory. As a next step, we can use the
random projection concept in the extended DMD and kernel DMD methods.

EB gratefully acknowledges funding from the Army Research Office
W911NF16-1-0081 (Dr Samuel Stanton) as well as from DARPA.

yes

References

yes

## References

  * Schmid (2010) Schmid, P.J.  Dynamic mode decomposition of numerical and experimental data.  Journal of Fluid Mechanics 2010, 656, 5–28.  doi:black10.1017/s0022112010001217.
  * Rowley et al. (2009) Rowley, C.W.; Mezić, I.; Bagheri, S.; Schlatter, P.; Henningson, D.S.  Spectral analysis of nonlinear flows.  Journal of Fluid Mechanics 2009, 641, 115–127.  doi:black10.1017/S0022112009992059.
  * Tu et al. (2014) Tu, J.H.; Rowley, C.W.; Luchtenburg, D.M.; Brunton, S.L.; Kutz, J.N.  On dynamic mode decomposition: Theory and applications.  Journal of Computational Dynamics 2014, 1, 391–421.
  * Dasgupta (2000) Dasgupta, S.  Experiments with Random Projection.  Proceedings of the 16th Conference on Uncertainty in Artificial Intelligence; Morgan Kaufmann Publishers Inc.: San Francisco, CA, USA, 2000; UAI ’00, p. 143–151.
  * (5) Johnson, W.B.; Lindenstrauss, J.; Schechtman, G.  Extensions of lipschitz maps into Banach spaces.
  * Bingham and Mannila (2001) Bingham, E.; Mannila, H.  Random Projection in Dimensionality Reduction: Applications to Image and Text Data.  Proceedings of the Seventh ACM SIGKDD International Conference on Knowledge Discovery and Data Mining; Association for Computing Machinery: New York, NY, USA, 2001; KDD ’01, p. 245–250.  doi:black10.1145/502512.502546.
  * Koopman (1931) Koopman, B.O.  Hamiltonian Systems and Transformation in Hilbert Space.  Proceedings of the National Academy of Sciences 1931, 17, 315–318, [https://www.pnas.org/content/17/5/315.full.pdf].  doi:black10.1073/pnas.17.5.315.
  * Bollt (2021) Bollt, E.M.  Geometric considerations of a good dictionary for Koopman analysis of dynamical systems: Cardinality, “primary eigenfunction,” and efficient representation.  Communications in Nonlinear Science and Numerical Simulation 2021, 100, 105833.  doi:blackhttps://doi.org/10.1016/j.cnsns.2021.105833.
  * Chen et al. (2012) Chen, K.K.; Tu, J.H.; Rowley, C.W.  Variants of Dynamic Mode Decomposition: Boundary Condition, Koopman, and Fourier Analyses, 2012.
  * Achlioptas (2003) Achlioptas, D.  Database-friendly random projections: Johnson-Lindenstrauss with binary coins.  Journal of Computer and System Sciences 2003, 66, 671–687.  Special Issue on PODS 2001, doi:blackhttps://doi.org/10.1016/S0022-0000(03)00025-4.
  * Kutz et al. (2016) Kutz, N.J.; Brunton, S.L.; Brunton, B.W.; Proctor, J.L.  Dynamic mode decomposition: data-driven modeling of complex systems; Society for industrial and applied mathematics., 2016.
  * hyc (2021) HYCOM + NCODA Gulf of Mexico 1/25° Analysis, 2021, [https://www.hycom.org/data/gomu0pt04/expt-90pt1m000].