Preconditioned dynamic mode decomposition and mode selection algorithms for large datasets using incremental proper orthogonal decomposition

This note proposes a simple and general framework of dynamic mode decomposition (DMD) and a mode selection for large datasets. The proposed framework explicitly introduces a preconditioning step using an incremental proper orthogonal decomposition to DMD and mode selection algorithms. By performing the preconditioning step, the DMD and the mode selection can be performed with low memory consumption and small computational complexity and can be applied to large datasets. In addition, a simple mode selection algorithm based on a greedy method is proposed. The proposed framework is applied to the analysis of a three-dimensional flows around a circular cylinder.

the mode selection of DMD modes. In the proposed framework, preprocessing steps are introduced before performing DMD and mode selection. Specifically, a small number of principal components of the input datasets are extracted using an incremental-type POD algorithm, then the input datasets are projected onto the principal component bases (namely, POD bases) to reduce the dimension. Subsequently, the DMD and mode selection steps are performed using the low-dimensionalized datasets. Mode selection is performed using a simple algorithm based on the greedy approach. The proposed framework is applied to a three-dimensional flow field around a circular cylinder to demonstrate its effectiveness.
First, as the preconditioning step, the low-dimensionalization of the input datasets is performed using POD. We use incremental-type POD algorithms because standard POD algorithms also need to store all the input datasets in memory. There are many incremental-type POD algorithms in the literature (e.g., Refs [12,13]). In this letter, the incremental POD proposed by Arora et al. [12] is used. It has been shown that the incremental POD is an effective tool for extracting dominant structures from fluid flow datasets. [6] We denote the input datasets as X = [x 1 x 2 · · · x N ], where a column vector x n ∈ R d represents the input data at time t = n∆t. The incremental POD updates the POD bases using the following algorithm each time new data x n is obtained. Suppose that C n−1 ∈ R l×l is a rankl approximation of a covariance matrix constructed using the datasets [x 1 x 2 · · · x n−1 ], and its eigendecomposition is C n−1 = U n−1 D n−1 U T n−1 , where the orthogonal matrix U n−1 ∈ R d×l and diagonal matrix D n−1 ∈ R l×l represent the POD bases and corresponding eigenvalues, respectively. The POD bases are obtained using the eigendecomposition of the approximate covariance matrix C n−1 . The update rule for C n−1 is as follows: First, an average µ n−1 can be updated by and then we definex n = U T n−1x n andx ⊥ n =x n − U n−1 U T n−1x n wherex n = x n − µ n .x n ,x n , andx ⊥ n represent the fluctuating component of x n , and the parallel and orthogonal components ofx n with respect to U n−1 , respectively. Using these variables, the covariance matrix C n−1 can be updated as follows [12]: where Therefore, the updated POD bases U n and corresponding eigenvalues D n are obtained using the eigendecomposition Q n = U S U T and If the rank of C n becomes greater than a user setting parameter r p , we delete the column and row of U n and D n that correspond to the smallest eigenvalues so that the rank becomes r p . The aforementioned update rule is repeated until all the input datasets are used. The memory consumption of the incremental POD is O(r p d), whereas that of the standard POD (and DMD) algorithms is O(N d). Therefore, we can apply the incremental POD to large datasets by setting the number of POD bases r p sufficiently small.
Next, we low-dimensionalize the input datasets. Before that, we perform the orthogonalization of the r d (= r p + 1) bases composed of the r p columns of POD bases U N and average µ N using the modified Gram-Schmidt (MGS) orthogonalization algorithm. We denote the obtained orthogonal bases as P . Using P , the preconditioning (i.e., low-dimensionalization of the input datasets) is achieved usingX = P T X, whereX ∈ R r d ×N is a small matrix. Note that this preprocessing step causes a loss of information in the input datasets. Generally, small spatial structures tend to disappear when low-dimensionalization using POD bases is applied to fluid flow data. [14,15] This loss of information decreases as rank r p increases. By contrast, the memory consumption and computational complexity of the incremental POD algorithm are O(r p d) and O(r 2 p d), respectively. Therefore, r p should be determined by considering this trade-off relation.
Most DMD algorithms can be applied to large datasets using the low-dimensionalized input datasets obtained using the aforementioned preconditioning step, for example, standard DMD [2], total least squares (TLS) DMD (tlsDMD) [16,17], and DMD with augmented input data. [18,19] In this letter, we use tlsDMD proposed by Hemati et al. [17] and Dawson et al. [16] Empirically, tlsDMD has good performance in terms of reconstructing input datasets using a small number of DMD modes because tlsDMD can accurately compute the growth rate and frequency of corresponding DMD modes.
DMD algorithms compute the approximate eigenvalues and eigenvectors of the linear operator A that satisfies tlsDMD computes the linear operator A using the TLS method, whereas standard DMD [2] uses the least squares method. As a first step of tlsDMD computation, set r d , that is, the number of columns of P , so that r d < (N − 1)/2 (this is not a strict restriction) and definẽ . Then, perform the reduced singular value decomposition ofZ and partition U d into four r d × r d sub-matrices as U d = matrix because Σ has 2r d non-zero singular values. Using U 11 and U 21 , the low-dimensional representation of the linear operator A is obtained as Finally, solve the eigenvalue problem ofÃφ =λφ, then the approximate eigenvalues of A and the corresponding DMD modes are obtained as respectively. Another DMD method for large datasets is a streaming DMD proposed by Hemati et al. [5] The main difference between the proposed method and streaming DMD is that in the latter, the updating of POD bases and projection of input datasets onto the bases are performed simultaneously, whereas in the present method, POD is performed as preprocessing for DMD (and mode selection). Therefore, the streaming DMD is suitable for online processing of streaming data. The advantage of the proposed framework is that the POD, DMD, and mode selection methods can be considered and performed separately. Therefore, various POD, DMD, and mode selection methods can be easily applied to this framework, and it is useful for analysis by trial and error.
Finally, algorithms for mode selection and the reconstruction of input datasets using the selected modes are shown. Selecting physically important DMD modes is important for the understanding of phenomena and constructing reduced order models. One promising approach to achieve this is to use compressed sensing, which was first introduced to DMD analysis by Jovanović et al. [9] The proposed method in this letter uses the compressed sensing approach. To adopt compressed sensing, the low-dimensionalized input datasetsX = P T X created by the preconditioning step are used instead of the raw large input datasets. Using the eigenvalues λ and eigenvectorsφ,x n can be written as the following expression: Therefore, we can approximateX using the following matrix form: In Eq. (12), the diagonal matrix D α and Vandermonde matrix V and represent the initial amplitudes and temporal variations of the corresponding DMD modes, respectively. Note thatΦ, D α , and V and are complex-valued matrices. Based on the compressed sensing approach, we select physically important DMD modes as the solution of the following optimization problem: where J(α) = X −ΦD α V and 2 .
α 0 is the number of non-zero elements in α = [α 1 α 2 · · · α r d ], and is a small positive number. Exact solutions of this optimization problem cannot be easily obtained because this problem is a combinatorial optimization problem. To obtain approximate solutions, Jovanović et al. [9] used the L1 regularization α 1 instead of the L0 regularization α 0 . This approach is the so-called least absolute shrinkage and selection operator [20] (LASSO). Algorithms based on the greedy approach are also often used to solve Eq. (13). It is well known that greedy approaches provide good solutions despite their quite simple algorithms. [21] This letter proposes the method based on the greedy approach. The proposed method selects the DMD mode (a column ofΦ) that minimizes a residual J S , and adds the corresponding index of the column to a support set S at each iteration step. We define the residual as J S = J(α sp ), where α sp is a solution of the following optimization problem: minimize α J(α) subj.to supp{α} = S, where supp{α} is a set of indices of α whose elements have non-zero values. This optimization problem was introduced by Jovanović et al. [9] to determine the optimized amplitude of α with a fixed sparsity structure, and they found that α sp is obtained using the following computation: where E is a matrix composed of unit column vectors for which the positions of the nonzero elements correspond to the zero components of α, for example, E = 0 1 0 0 0 1 T for α = α 1 0 0 T . F and g are F = (Φ * Φ ) • (V and V * and ) and g = diag(V andX TΦ ), respectively. An asterisk denotes the conjugate transpose, a overline denotes the complex conjugate, • denotes elementwise multiplication, and diag(·) denotes a vector whose elements are the diagonal elements of the matrix. Note that the computational cost to calculate J S is not expensive because the dimension ofX is typically small (e.g., r d = O(10 1 ) and N = O(10 2 )) as a result of the preconditioning step. If a stopping criterion set by a user is satisfied, the iteration is terminated. In this letter, α 0 = K is used as the stopping criterion, where K is a user setting parameter and the number of DMD modes to be selected. The proposed algorithm is summarized in Table  1. Initialize: initial support S = ∅ Repeat until stopping criterion is met: 1)Compute J S∪{j} for all the column indices j / ∈ S where, J S∪{j} = min α J(α) subj. to supp{α} = S ∪ {j} 2)Add j 0 = arg min j J S∪{j} to the support S S ← S ∪ {j 0 } Finally, the reconstructed input datasets X R using the selected modes can be computed as The proposed framework for DMD analysis is summarized in Table 2. To demonstrate the effectiveness of the proposed framework, we applied the framework to the analysis of three-dimensional laminar flow around a circular cylinder. The Reynolds and Mach numbers based on the diameter D of the cylinder and freestream velocity U ref were Re = 350 and M = 0.2, respectively. The numerical simulation was performed using in-house code that has been verified for several analyses. [22,6] The simulation was performed with a sixth-order compact finite difference scheme [23,24] with tenth-order filtering [25] for spatial discretization, and a third-order TVD Runge-Kutta method [26] for time marching. Periodic boundary conditions were applied in the spanwise direction, with the length of the computational domain L z = 3D. The computational grid was composed of 9.5×10 6 grid points. The number of snapshots was N = 800, and the time interval for each snapshot was ∆t = 0.25. Each snapshot x n was composed of three components of the velocity in the x, y, and z directions of each grid point.  Although the Reynolds number is relatively low, it is not so easy to extract dominant dynamics from this flow field because the flow includes chaotic behavior.
Memory consumption for the present analysis was approximately 24 GB, with r d = 51. This memory requirement is sufficiently small for performing the analysis on recent workstations. Note that if we analyze the present input datasets using standard DMD, over 10 times more memory is required.   Figure 2 shows the distribution of eigenvalues obtained by the present DMD. The eigenvalues selected by the mode selection algorithm are also shown. The DMD mode of (σ, St) = (0, 0) was first selected, where σ = Real{log(λ)}/∆t and St = Imag{log(λ)}/(2π∆t). This mode corresponds to the mean flow field. The first selected oscillation modes had a frequency of St = 0.20. This mode corresponds to the well-known two-dimensional vortex shedding, as shown in Figure 3a. The second oscillation modes had a relatively low frequency of St = 0.07. Interestingly, according to Figure 3b, this mode seems to represent oblique streamwise vortex phenomena. [27,28] Though not shown here, the third oscillation modes (St = 0.13) also represent oblique streamwise vortices. K-selected modes obtained by the present mode selection algorithm and LASSO [9] are shown. Figure 4 illustrates the dependence of the absolute values of the optimized initial amplitudes α sp on the growth rate σ. This figure shows that the modes with small growth rates selected by the mode selection algorithms had large initial amplitudes. This means that, even if a mode has a small growth rate, the mode can behave as a dominant phenomenon if its initial amplitude is sufficiently large. Therefore, both growth rates and initial amplitudes are important in selecting dominant modes that represent the input datasets. Comparing the proposed mode selection algorithm with the previous algorithm (LASSO) [9], it can be seen that the proposed algorithm tends to select modes that have large initial amplitudes rather than large growth rates.  Figure 5 shows the effect of the number of selected modes K used to reconstruct input datasets on the variance of the reconstruction error ε, where ε is defined by the following equation: whereX is the perturbation component of X, that is, Figure 5 clearly shows that the error ε monotonically decreases as K increases. In particular, the first oscillation modes made a large contribution to the reconstruction of the input datasets. Additionally, the second and third oscillation modes also made a relatively large contribution. This means that the proposed mode selection algorithm correctly selected the dominant DMD modes. Furthermore, in this fluid dataset case, the proposed algorithm had fewer reconstruction error than LASSO.[9] Figure 5 indicates there are still certain errors, even if all the modes are used for the reconstruction because of chaotic behavior included in the present flow. Generally, most DMD algorithms cannot manage the dynamics that cannot be approximated by a local linear operator. Additionally, the preconditioning step, that is, the dimensionality reduction using POD bases, causes a loss of information, mainly about small spatial structures. It is necessary to note that the proposed method is suitable for extracting the dynamics of large spatial structures where local linear approximation is valid. The temporal history of the original and reconstructed flows at the point (x, y, z) = (1.51D, 0.51D, 0D) is shown in Figure 6. We can confirm that large parts of the fluctuations of the velocity are well reproduced using seven modes, although the velocity in the z direction has a relatively large reconstruction error because of its chaotic behavior.
This letter proposed a simple and efficient framework of DMD and its mode selection for large datasets. The analysis of three-dimensional flow around a circular cylinder showed that the proposed framework can perform DMD and mode selection with low memory consumption,