Dynamic mode decomposition using a Kalman filter for parameter estimation

A novel dynamic mode decomposition (DMD) method based on a Kalman filter is proposed. This paper explains the fast algorithm of the proposed Kalman filter DMD (KFDMD) in combination with truncated proper orthogonal decomposition for many-degree-of-freedom problems. Numerical experiments reveal that KFDMD can estimate eigenmodes more precisely compared with standard DMD or total least-squares DMD (tlsDMD) methods for the severe noise condition if the nature of the observation noise is known, though tlsDMD works better than KFDMD in the low and medium noise level. Moreover, KFDMD can track the eigenmodes precisely even when the system matrix varies with time similar to online DMD, and this extension is naturally conducted owing to the characteristics of the Kalman filter. In summary, the KFDMD is a promising tool with strong antinoise characteristics for analyzing sequential datasets.


I. INTRODUCTION
Recently, fluid analysis has been conducted with high-resolution numerical simulations and experimental measurements. Such simulations and experiments provide large-scale data, so precise extraction is necessary in order to understand and model essential phenomena from the data provided. Mode decomposition 1 is a useful method for conducting such processes.
Proper orthogonal decomposition (POD) proposed by Lumrey has been applied to fluid analyses, especially turbulent analyses. 2,3 Modes obtained by POD are known to be orthogonal to each other. Furthermore, the original flow can be reconstructed with a limited number of modes. Proper orthogonal decomposition is equivalent to principal component analysis (PCA) and Karuhunen-Loéve expansion. Note that fluid phenomena can be approximated and modeled by several methods using POD modes, e.g., the Galerkin projection method.
Proper orthogonal decomposition is optimum from the viewpoint of energy reconstruction with fewer modes, although the POD modes are not derived from the original fluid equations.
Moreover, global linear stability analysis (GLSA) 4-6 is a major method that can extract the eigenmodes of perturbations using governing equations (e.g., the Navier-Stokes equations) linearized around a nonlinear steady state. If GLSA is applied to the Navier-Stokes equations linearized on a steady state solution, the most unstable eigenmodes are extracted and are used to judge whether the steady state solution is stable. The modes obtained by GLSA satisfy the original linearized equation(s), although GLSA is more complicated than POD.
Note that these modes are generally not orthogonal.
Dynamic mode decomposition (DMD) 7 is developed as an intermediate method of POD and GLSA and is applied to numerous applications. [8][9][10] In DMD, a sequential dataset of an unsteady flow solution, which is generally nonlinear, is given as input. This dataset is considered to be explained by a linear system (x k+1 = Ax k ), and the eigenvalues and corresponding eigenvectors of A are calculated for mode decomposition. Although these modes are generally not orthogonal, they represent the single-frequency response with amplification or damping. This means that physical phenomena of the DMD mode can be understood more simply than those of POD. The original DMD 7 uses singular value decomposition (SVD) to compute a low-rank approximation of matrix A. In exact DMD (EDMD), a Moore-Penrose pseudoinverse matrix is applied instead. 11 Various studies have recently been conducted in this field. Total least-squares DMD (tlsDMD) 12 improves the 2 estimation accuracy of the eigenvalues and eigenmodes using truncated SVD (POD) for the combined data of successive two snapshots that filter out the critical noise for estimation of dynamics. Sparsity promoting DMD (spDMD) 13 chooses modes with which the original flows can be effectively reconstructed in the framework of DMD by introducing sparse modeling and compressing sensing ideas. Recently, for identification of a time-variant system, a short-time DMD method is proposed. 14 As mentioned above, DMD is a more promising method for extracting the modes that can directly describe the system dynamics, as compared to POD, and further development is expected.
Following previous studies, DMD is reconsidered for parameter estimation of matrix A or system identification in the present study. In other words, if matrix A of DMD is considered to be a kind of filter, then the DMD problem is regarded as coefficient identification of the filter. Conventional approaches to solve this kind of problem are a recursive least-squares (RLS) method and a Kalman filter method. 15 In the present study, we propose a novel method by which to use a Kalman filter to identify matrix A.
The following advantages are expected when adopting a Kalman filter for the estimation of DMD modes.
• More arbitrary treatment for denoising when the noise characteristics are known, • System identification of the transient system, and • Most likelihood estimation of the system. The first two of the above advantages are demonstrated using the Kalman filter in the present paper. With regard to the first advantage, the system is considered to be estimated more precisely by the Kalman filter than by standard DMD methods if the observation noise covariance is known in advance. Data for space science, astronomy, and meteorology are contaminated by severe time-dependent noise, the characteristics of which are known, and system identification based on such observations appears to be useful. This type of problem, in which the noise level is known in advance, was solved by the group of astrophysics 16 .
Moreover, the proposed method will help in conducting DMD for extremely severe measurements at low-signal-noise ratios, such as the measurement of compressible turbulence. In addition, with regard to the second advantage, the Kalman filter can be adopted inherently for a time-variant system. Matrix A is expected to be naturally identified, even if the matrix is time dependent, as shown in a previous study. 14

3
In Section II of the present paper, we introduce the Kalman filter for DMD and the fast algorithm in combination with POD to improve the poor computational efficiency of the straightforward implementation. In Section III, test problems are solved and two of the above-described advantages are demonstrated.

A. Proposed Algorithm
In the present study, a discretized system in the temporal direction is considered, as is usual in standard DMD methods. The subscript k represents a kth quantity in discretized time, k∆t, where ∆t is the time interval of snapshots. A linear temporal evolution system is considered: in vector form or in tensor form with the Einstein summation convention. Here, A = (a ij ) ∈ R n×n is a system matrix, x = x i are fluid variables, and n is the dimensionality of the fluid variables. Only the snapshots of the system, i.e., a dataset assumed to be generated by A, can be observed.
The snapshot data x can also be expressed as follows: where we define X = Z 1:m−1 and Y = Z 2:m , the elements of which are expressed as x k and y k (k = 1, · · · , m − 1).
Then, we consider a system identification problem. Each element of matrix A is estimated 4 in the present study, and the parameter vector θ is introduced as follows: where θ is considered to be a constant or slowly and randomly varying parameter vector according to the system noise, and the time evolution of θ can be given as follows: where F = I is an identity matrix for a constant or slowly and randomly varying parameter vector, and v is the system noise.
An observation equation for the next input y k = x k+1 = Ax k is given as follows: whereas 5 Since H k varies with the time step, the system is a linear time-variant system, and the resulting algorithm of the Kalman filter is standard for a linear time-variant system and not a special implementation. Based on these equations, a standard linear Kalman filter can be used with θ.
Following the theory of a Kalman filter, a covariance matrix regarding a priori estimation P k|k−1 can be obtained using the covariance matrix of one step earlier, i.e., P k−1|k−1 , where Q is a covariance matrix regarding system noise, and a system matrix F k becomes an identity matrix from Eq. (6). In a priori estimation, θ does not change because of the The state variables are updated by the Kalman gain when observation takes place. A noise covariance matrix after observation, S k , is given by, where R k is a covariance matrix of observation noise and is generally time dependent. The Kalman gain is then directly computed by The amount of modification of state variables θ can be computed as follows: where A is generated from θ k . Note that during this process, snapshot y k is newly observed.
The quantity y k observed here is used as x k in the next time step to construct the observation matrix.
The covariance matrix after the observation is updated by, We update the state variable as follows: The disadvantage of this formulation is that inversion of the large matrix S k with a dimension of n 2 is required. In the next section, a novel algorithm with extremely low computational cost is introduced. 6

B. Fast Algorithm
The periodicity and sparsity of the matrices appearing in the previous algorithm are used, and the problem is further simplified.
The following assumptions are introduced for simplicity: 1. The initial covariance matrix P is assumed to be a block diagonal matrix, and all of the diagonal matrices are identical.
2. The covariance matrices of observation and system noises, R and Q, are assumed to be block diagonal matrices and all of the diagonal matrices are identical (where R k = r k I in this case) The above conditions are expressed as follows: where subscript (1, 1) represents the first block element in an original matrix.
In this case, S k = s k I is obtained, and its value s k is, 7 The Kalman gain becomes a vector and is expressed as follows: Note that the dimensionality of the Kalman gain derived here is n × 1. When the Kalman gain matrix is multiplied by Eq. (13), the element of the Kalman gain matrix is copied in the column direction, as follows: and its dimension is expanded to n × n as a result. Here, the subscript (1) indicates the first column in the original matrix.
The covariance matrix after observation can be updated as follows:

C. Combination with truncated POD
Although the use of the algorithm described in the previous subsection helps us to compute θ quickly, matrix P and state variables θ require memories of n 2 variables, and, for some fluid problems, it is impossible to store all of the matrix variables. Therefore, in the present study, truncated POD (truncated SVD) is used as a preconditioner and the number of degrees of freedom are reduced for applying the Kalman filter to the dataset of the fluid system. In the present study, 1) the batch POD is first applied, and 2) the proposed Kalman filter is then applied to the amplitude of each POD mode. Finally, the mode shape of the fluid system is recovered by multiplying the spatial POD modes.
More concretely, POD is applied to a data matrix: whereas U and V contain the spatial and temporal POD modes, respectively, as row vectors.
The r-rank approximation is obtained by applying the truncated POD, as follows: andX andỸ are treated similarly to X and Y in the Kalman filter DMD procedures. In addition, for a more flexible online procedure, we can use the following formulation: where the left singular vector is assumed to be fixed. After obtaining the right eigenvector of the reduced system using the Kalman filter DMD, it is necessary to recover the original dimension by multiplying matrix U.
In the present paper, Eqs. 26 and 27 are used for the truncated POD. This procedure is adopted for many-degree-of-freedom problems (n > 10, 000), and it is not used unless otherwise mentioned.

III. NUMERICAL EXPERIMENTS AND DISCUSSION
The fast Kalman filter algorithm described in Section II B is adopted in the numerical experiments below. White noise with N (0, σ 2 ) is added to sequential snapshot data after the transformation described above. The variance (σ 2 ) is time dependent, and σ 2 = σ 2 0 · (1.01 + cos (π∆tk)), where σ 2 0 = 0.05. For the initial adjustable parameters of the Kalman filter, the diagonal parts of the covariance matrix are set to 10 3 , and all of the elements of Q are set to zero. The same computation is conducted 100 times using different time-dependent seeds for the random number generation.

A. Static System Identification with Noise
First, the case in which σ 2 0 is set in the diagonal part of R is investigated. Figures 1  and 2 show the results with 100 (m = 100) and 500 (m = 500) snapshots, respectively.
The dashed line in the figures represents a unit circle. These figures show that the results obtained using the Kalman filter agree well with those obtained using the standard DMD.

B. Static System Identification with Noise and its Known Characteristics
Next, the performance of KFDMD for the problem with time-dependent R with the known characteristics of σ 2 = σ 2 0 · (1.01 + cos (π∆tk)) is investigated, whereas the other problem settings are the same as in the previous problem. Figures 3 and 4  Both methods produce the dynamic mode of Kármán vortex shedding, and these results show that KFDMD can estimate dynamic modes similar to DMD, when noise is absent.

D. Static Fluid System Identification with Noise
In this subsection, the same data as in the previous problem are used, but with known noise characteristics. The noise strength is determined using time-dependent R with σ 2 =

E. Dynamical System Identification with Noise
Next, dynamical system identification using KFDMD is conducted. The system has only two eigenmodes, the real part of which is zero, and the time-dependent imaginary part is given by Im(λ) = 2πf , where the frequency f is given by f = 1 + k∆t. The system is then projected to a new system of dimension r = 20 using the same method stated in the previous section. The time-dependent white noise, the magnitude of which is the same as the static system identification, is added. Then, σ 2 0 is fixed to 0.05, i.e., the element of the R matrix is not modified during this identification procedure. Figure 9 shows the results of the numerical experiment. In the figure, the horizontal axis represents the number of observed snapshots. Until the snapshots are sufficiently observed, the results of the proposed method are not converged. Therefore, it is difficult to identify which modes are appropriate to plot.
For this reason, the estimated frequency is plotted only after the convergence. Remarkably, the proposed method can naturally track the true frequency.

IV. CONCLUSIONS
A novel dynamic mode decomposition method based on the Kalman filter was proposed.
Our numerical experiments revealed that the proposed method can estimate matrix A more precisely than the standard dynamic mode decomposition method. This tendency is especially strengthened if the characteristics of the observation noise that contaminates the system to be identified are known when identification is performed. Furthermore, the proposed method can identify time-dependent systems in which the matrix is transited with time, and this expansion is naturally conducted owing to the characteristics of the Kalman filter. Note that all of these properties are preferred in data analysis. The dynamic mode decomposition method based on the Kalman filter is a promising tool for analyzing a noisy dataset.