Randomized Projection Learning Method forDynamic Mode Decomposition

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 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. Keyword: Koopman Operator, Dynamic Mode Decomposition(DMD), Johnson-Lindenstrauss Lemma, Random Projection, Data-driven method.


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) [1][2][3] 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 projection [4] 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 [5] 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×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 O(MLN) while SVD has computational complexity O(NM 2 ) when M < N [6]. 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.

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.

Koopman Operator
Consider a discrete-time dynamical system where S : M → M and 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 M as phase space. The associated Koopman operator is described as the evaluation of observable functions ( Fig. 1) ψ : M → R in function space F . Instead of analysing the individual trajectories in phase space, the Koopman operator operates on the observations [2,3,7,8].
Definition 1 (Koopman operator [7]). The Koopman operator K for a map S is defined as following composition, on the function space F .
It is straight forward to prove [7], Figure 1. This figure shows the behavior of Koopman operator K in observable space F associated with a dynamical system S. The Koopman operator evaluate the observable ψ at downstream or future for ψ 1 , ψ 2 ∈ F and a, b ∈ C, and therefore the Koopman operator is linear on 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, 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 [8]. 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 K is given by then vector-valued observables g g g : M → R N (or C N ) can be represented by where φ i φ i φ i ∈ R N (or C N ) are the vector coefficients of the expansion and called "Koopman modes"(here we assumed that components of g g g lie within the span of the eigenfunctions of K). Note that the observable value at time n + 1 is given by 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.

Dynamic Mode Decomposition
The dynamic mode decomposition is a data-driven method to estimate the Koopman modes from numerical or experimental data [2]. Suppose dynamics are govern by eq. (1) for any state x x x and vector valued measurements are given by observable g(x) g(x) g(x) ∈ R N . For a given set of data , the Koopman modes and eigenvalues of the Koopman operator can be estimated through solving the least-squares problem and K = YX † (here X † is the pseudo-inverse of X) is defined as the "Exact DMD" operator [3]. The eigenvalue (λ) of K is an approximation of an eigenvalue (λ) of K; the corresponding right eigenvector(φ) is called the DMD mode and approximates the Koopman mode (φ). Then the observable value g(x(t)) g(x(t)) g(x(t)) at time t can be modeled as where r is the number of selected DMD modes and demonstrates the finite dimensional approximation for vector-valued observable g g g under the Koopman operator. Based on this decomposition, data matrices can be expressed as for i = 1, 2, . . . r, j = 1, 2, . . . , M and Λ = diag{λ 1 , λ 2 , . . . , λ r }. Note that with the above decomposition K = YX † = ΦΛTT † Φ † . We will suppose K has distinct eigenvalues λ i , columns of X are linearly independent and r ≤ 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 R N (N >> 1 and roughly O(10 15 ) to O(10 10 )), but the number of snapshots or time steps (M) are small and often it is O(10 3 ) to O(10 1 ) [9]. Hence, computing the spectrum of matrix 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 ≤ L ≤ M << N, therefore need to estimate the spectrum of 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.

Random projection
The random projection method is based on the Johnson-Lindenstrauss lemma which is applied by many data analysis methods.
Theorem 1 (Johnson-Lindenstrauss Lemma [5]). For any 0 < < 1 and any integer M > 1, let L be a positive integer such that L ≥ L 0 with L 0 = C ln M 2 , where C is a suitable constant (C ≈ 8 in practice,C = 2 is good enough). Then for any set X of M data points in R N , there exists a map f : R N → R L such that for all x 1 , x 2 ∈ X, Theorem 2 (Random Projection [4]). For any 0 < , δ < 1 2 and positive integer N, there exists a random matrix of B of size L × N such that for L ≥ L 0 with L 0 = C ln(1/δ) 2 . and for any unit-length vector x ∈ R N 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 × N with the L ≥ C ln N 2 which provides − isometry to R N . (See Fig. 2 for details).

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

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 ∈ R L×N be any rank L projection matrix, then dimension of data matrices X, Y ∈ R N×M can be reduced to L × M by the projection X L = PX, Y L = PY. The DMD operator on projected space (see Fig. (2)) is given by, where K = YX † is the DMD operator on original space.
Proof. Let (λ L , φ L ) be an eigenpair ofK. ThenKφ L = λ L φ L and by eq. (11), Kφ − λ L φ = 0 is a solution to the above equation, λ L is an eigenvalue and the corresponding eigenvector In other words, we can lift up the dimension of eigenvectors in projected space by P † 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φ = Y(PX) † φ L is an eigenvector of K for corresponding non-zero eigenvalues.
Also notice, Pφ = PY(PX) † φ L =Kφ L and thereforφ 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 Next, the focus move to the spatial-temporal decomposition of projected data matrices by spectrum of the DMD operator. Note that observable value g(x(t)) g(x(t)) g(x(t)) at time t can be modeled as g( i and similar to the eq. (9), data can be decomposed as This decomposition leads to K = YX † = P †Φ ΛTT †Φ † P and if r ≤ L all the non-zero eigenvalues and corresponding eigenvectors of 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 (Projection by SVD). Commonly used projection matrix is based on SVD of the input matrix X = UΣ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 asK = U * YVΣ −1 .
Remark (Standard DMD and Exact DMD). Let eigenpair of a SVD basedK = U * YVΣ −1 be given by (λ, φ L ). Remark. QR decomposition based projection methods on both input and output data [X Y] can be used [3].
In this paper we are proposing a simple random projection based method to estimate the spectrum of the DMD operator.

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×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 byK = RY(RX) † , and if an eigenpair ofK is given by (λ, φ L ), then the eigenpair of K is given (λ, Y(RX) † φ 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 [10].
One time step forecasting error for any given snapshot by using rDMD algorithm can be bounded by using the JL theory.

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.

Logistic Map
We will first consider a dataset of 300 snapshots from a logistic map, 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 ≥ C ln(300) 2 ≈ 34.22 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 ( fig. (5) (a)) and dimension of the projected space (( fig. (5) (b))). Since the rank of the snapshot matrix is 256, any L ≥ 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 ( or the projected dimension) from the JL theory.

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 [11]) , where γ j 's are constants, and let Φ j (x) = j sech(0.1x + j) and T j (t) = e iγ j 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 γ j = j (see fig. (6)). As discussed (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 and (b) shows the error with dimenstion of the projected space that will guarantee the bound for this example.  z(x, t)) plot for all x, t values at 20000 × 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 Φ j (x) = j sech(0.1x + j) and T j (t) = e ijt from the given data constructed by z(x, t).
in the previous section, if L ≥ 20 then those expected modes can be isolated and there exist eigenvalues λ j of rDMD operator such that ω j = ln λ j = ji for j = 1, 2, . . . , 20 (see fig. (7)). Furthermore, we are expecting corresponding rDMD modes equal to spatial variables of the model such that As expected, we noticed that calculated modes have negligible error when the dimension of projected space L ≥ 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 ≥ r = 20.
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  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.
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, Where E L+1 = ∑ m L=j+1 b j φ j (x)e ω 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)).

Gulf of Mexico
In this example we will consider the data from HYbrid Coordinate Ocean Model (HYCOM) [12] 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 × 347 grid points) data for 10 days (240 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π/ Im(ln(λ i )) of the modes(see table (2)). We constructed the snapshot matrix by stacking the snapshots of velocity components (u, v) in each column to perform the DMD analysis. 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.
Notice that the dimension of the snapshot matrix is 375454 × 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.

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 ≥ m. The SVD projection-based exact DMD method still needs to calculate the SVD of a high dimensional(roughly 10 10 × 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 × 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 × 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.
Funding: EB gratefully acknowledges funding from the Army Research Office W911NF16-1-0081 (Dr Samuel Stanton) as well as from DARPA.