Bayesian Dynamic Mode Decomposition

1 Introduction

Analyzing nonlinear dynamical systems is fundamental for the understanding of complex phenomena in a variety of scientific and industrial fields. For example, the analysis of the Navier-Stokes equation has been one of the fundamental problems for understanding fluid flows. One of popular approaches for this purpose is decomposition of the dynamics into multiple components based on some criteria; the individual aspects of complex phenomena can be investigated by examining each decomposed component. For example, proper orthogonal decomposition (POD) (see e.g. [Holmes et al., 2012]) decomposes the dynamics into orthogonal modes that optimally capture the energy of the dynamics, and it has been extensively applied in physics [Bonnet et al., 1994; Noack et al., 2003]. In machine learning and pattern recognition, a method equivalent to data-driven POD is known as principal component analysis (PCA) and has been applied for modal decomposition and dimensionality reduction of a wide variety of numerical datasets [Jolliffe, 2002].

Dynamic mode decomposition (DMD) [Rowley et al., 2009; Schmid, 2010; Kutz et al., 2016a] has been attracting attention in various fields of science and engineering as an approach for the above purpose that works without explicit knowledge on the governing equations (see Section 5

for some examples). Although DMD is a data-driven decomposition technique like PCA, it generates modes that are directly related to the dynamics behind the data; thus, these modes are a useful tool for the diagnostics of complex dynamic phenomena. Insofar, several algorithmic variants of DMD have been utilized according to given datasets or purposes. However, most of these variants are deterministic (i.e., lack corresponding probabilistic models), and thus it could be difficult to incorporate uncertainty in data into the analysis. Building a probabilistic model for DMD enables us to treat the data statistically and consider observation noises explicitly, as well as to enrich the DMD techniques systematically by modifying the involved probabilistic distributions.

In this paper, we propose Bayesian DMD, which provides a principled way to transfer the advantages of the Bayesian formulation into DMD. To this end, we first develop a probabilistic model corresponding to DMD, whose maximumlikelihood estimator coincides with the solution of DMD. Then, we provide the Gibbs sampler for the posterior inference in Bayesian DMD. Due to the Bayesian treatment, we can infer posteriors of DMD-related quantities, such as dynamic modes and eigenvalues. Moreover, we can consider variants of DMD within the unified Bayesian framework by modifying the probabilistic model. In particular, we discuss the case of using a sparsity-promoting prior for dynamic modes, which allows us to automatically determine the number of modes in the light of data.

The remainder of this paper is organized as follows. We briefly review the underlying theory of DMD and its numerical procedure in Section 2. The probabilistic model for DMD is described in Section 3, and based on that model, Bayesian DMD is introduced in Section 4. In Section 5, we review the related work. In Section 6, we show the experimental results with synthetic and real-world datasets. This paper is concluded in Section 7.

2 Background

We briefly review the underlying theory of DMD, the spectral decomposition of nonlinear dynamical systems based on the Koopman operator. We recommend readers to consult papers such as [Mezi´c, 2005; Budi˘si´c et al., 2012; Mezi´c, 2013] for more details on the theory of the Koopman operator and decomposition based on it.

2.1 Koopman Spectral Analysis

Our interest lies in analyzing a (possibly, nonlinear) discretetime nonlinear dynamical system

whereM is the state space and t ∈ {0}∪N is the time index. Let g : M→ C be a measurement function (observable) in a function space G. Koopman operator K : G → G is defined as an infinite-dimensional linear operator such that

Defining K, we can lift nonlinear dynamics f to a linear (but infinite-dimensional) regime. However, it is difficult to numerically compute K from data because of its infinite dimensionality. Nonetheless, the infinite-dimensional system can further be lifted to a finite-dimensional one as follows.

Suppose that there exists an invariant subspace of K, i.e., a subspace G ⊂ G such that

and that a set of observables {g1, . . . , gn} (n < ∞) spans G. Moreover, consider the restriction of K to G and denote it by K : G → G. NowK is a finite-dimensional linear operator. Note that K has a matrix-form representation K with respect to {g1, . . . , gn}, i.e.,

wherein g = [g1 · · · gn]T . In the sequel, we assume that such invariant subspace G and the set of observables that spans G exist.

Since K is linear, modal decomposition based on it can be derived as follows. Let ϕ : M→ C be an eigenfunction of K with eigenvalue λ ∈ C, i.e.,

This eigenfunction with respect to {g1, . . . , gn} is obtained by ϕ(x) = zHg(x), wherein z is the left-eigenvector of K corresponding to eigenvalue λ. Moreover, let vi and zi respectively be the right- and the left- eigenvectors of K corresponding to eigenvalue λi, and suppose that they are normalized so that vH i zi′ = δii′ , without loss of generality. If all the eigenvalues of K are distinct (i.e., their multiplicity is one), any values of g are expressed as

Applying K to both sides of Eq. (1), we have

In this way, modal decomposition of observables via the Koopman operator is given by

wherein the values of g are described as a sum of Koopman modes w, whose temporal frequency and decay rate are given by ∠λ and |λ|, respectively. From the above, we could see that, unlike the classical modal decomposition of linear time invariant systems, this theory is applicable even to nonlinear dynamical systems.

2.2 Dynamic Mode Decomposition

DMD [Rowley et al., 2009; Schmid, 2010; Kutz et al., 2016a] is a numerical decomposition method, and it coincides with modal decomposition via the Koopman operator under some conditions. Suppose we have data matrices:

where g is again a concatenation of n observables, and m+1 is the number of snapshots in the dataset. The popular implementation of DMD [Tu et al., 2014] is shown in Algorithm 1, with which we finally compute the eigenvectors of a matrix A = Y1Y † 0 corresponding to its nonzero eigenvalues, wherein Y † 0 denotes the Moore-Penrose pseudoinverse of Y0. If the components of g span a subspace invariant to K and all modes are sufficiently excited in the data (i.e., rank(Y0) = dim(G)), then the dynamic modes computed by Algorithm 1 converge to the Koopman modes in modal decomposition (3) in the large sample limit.

Algorithm 1 (DMD [Tu et al., 2014]).

1. Compute the compact SVD of Y0 = UrSrV H r .
2. Define a matrix ˜A= U H r Y1VrS-1 r .
3. Calculate eigendecomposition of ˜A, i.e., compute ˜w and λ such that ˜A ˜w = λ ˜w.
4. Return dynamic modes w = λ-1Y1VrS-1 r ˜w and corresponding eigenvalues λ.

Note that, for DMD to obtain the theoretical rationale, the components of g need to span an (approximately) invariant subspace. There are several studies that address this issue (e.g., the use of nonlinear basis functions [Williams et al., 2015], reproducing kernels [Kawahara, 2016], and delay coordinates [Susuki and Mezi´c, 2015; Arbabi and Mezi´c, 2016]). In this paper, however, we simply assume data are generated with observables that intrinsically span an (approximately) invariant subspace, as in the previous studies on DMD.

3 Probabilistic DMD

We develop a probabilistic model associated with modal decomposition via the Koopman operator (Eqs. (1) and (2)). The maximum-likelihood estimator (MLE) of this model coincides with the solution of DMD in the no-noise limit. As will be described in the next section, this probabilistic model forms the foundation for Bayesian DMD.

3.1 Generative Model

Let yℓ,j ∈ Cn be the j-th column of Yℓ in Eq. (4) plus observation noise, for ℓ = 0, 1. Following the relations in Eqs. (1) and (2), the probabilistic DMD model for such data can be given by

where we assume that the observation noise is Gaussian, and CN (µ, σ2I) is the complex Gaussian distribution [Goodman, 1963] whose density is defined as

.

Here, w1:k, λ1:k, and σ2 are the parameters to be estimated (λ1:k denotes a set {λ1, . . . , λk}), and k is the tunable hyperparameter that determines the number of modes (usually k ≤ n). In addition, we treat ϕi,j as a latent variable with the standard Gaussian prior

3.2 Maximum-likelihood Estimator

To derive the MLE of probabilistic DMD, let us rewrite likelihood (5) in a matrix form, i.e.,

where we use notations as follows:

Marginalizing out ϕ with prior (6), we have

In the following, we describe the relationship between probabilistic model (7), total-least-squares DMD [Dawson et al., 2016], which is a "noise-aware" variant of DMD, and standard DMD (Algorithm 1). In short, their estimation results coincide in the no-noise limit.

Proposition 1. Suppose we have a dataset that is possibly contaminated by observation noises E:

and let ˆ Y = [ ˆ Y T 0 ˆ Y T 1 ]T and Σy = m-1 ˆ Y ˆ Y H . In addition, let (σ2 )⋆ , W ⋆ and Λ⋆ be the MLEs of Eq. (7) given ˆ Y . If k = n, then the columns of W ⋆ and the elements of diag(Λ⋆ ) coincide with the dynamic modes and eigenvalues obtained by total-least-squares DMD, respectively.

Proof. Following [Tipping and Bishop, 1999], the MLEs for probabilistic model (7) are given as

with Uk = [u1 . . . uk] and Mk = diag(µ1, . . . , µk), where µi is the i-th largest eigenvalue of Σy with corresponding eigenvector ui, and R is an arbitrary unitary matrix. If k = n, we haveW ⋆Λ⋆ (W ⋆ )-1 = U1,nU -1 0,n , where U0,n comprises the first n rows and U1,n comprises the last n rows of Un. Hence the columns of W ⋆ and the elements of diag(Λ⋆ ) are obtained by the eigendecomposition of U1,nU -1 0,n , which is exactly the same procedure with the one in total-least-squares DMD [Dawson et al., 2016].

Proposition 2. If Y0 and Y1 are linearly consistent [Tu et al., 2014] and there is no observation noise (i.e., E = 0), then the estimation results of total-least-squares DMD coincides with those of standard DMD (Algorithm 1).

Proof. From the definition of the linear consistency [Tu et al., 2014], when there is no observation noise, rank(Σy) = n. Hence, m- 1 2 ˆ Y = UnM 1 2 n V H n (Vn is comprising first n right singular vectors of m- 12 ˆ Y ). Consequently,

which shows the equivalence of the outputs of total-leastsquares DMD and standard DMD.

4 Bayesian DMD

For the Bayesian treatment of DMD, we consider the following priors on the parameters in probabilistic model (5). First, we put a Gaussian prior on w1:k:

with an inverse gamma hyperprior on v2 i,d (d = 1, . . . , n):

whose shape parameter is αv and rate parameter is βv. Moreover, we consider priors on λ1:k and σ2 as

A graphical model of Bayesian DMD is shown in Figure 1.

4.1 Posterior Inference by Gibbs Sampling

The conditional probabilistic distribution on each latent variable in the above model becomes a complex Gaussian or an inverse gamma distribution and thus is easy to sample. Consequently, we can develop a Gibbs sampler for inferring the latent variables. In the sequel, we use notations

The sampling procedures of the Gibbs sampler are summarized in Algorithm 2. In the following, ¯λ denotes the complex conjugate of λ, and [wi]d denotes the d-th element of wi. In our implementation, the hyperparameters α and β were set to 10-3 .

Algorithm 2 (Gibbs sampling for Bayesian DMD).

1. Sample w1:k from CN(mwi ,P -1 wi ) with

1. Sample v2 1:k,1:n from InvGamma (av2 i,d , bv2 i,d ) with

1. Sample λ1:k from CN(mλi , p-1 λi ) with

1. Sample ϕ1:m from CN(mϕj ,P -1 ϕj ) with

1. Sample σ2 from InvGamma(aσ2, bσ2) with

1. Repeat (1)-(5) for a sufficient number of iterations.

4.2 Sparsity-promoting Prior

One of the difficulties when applying DMD to noisy data in practice is how to determine the effective number of dynamic modes. [Jovanovi´c et al., 2014] proposed sparsity-promoting DMD, in which the number of dynamic modes are determined by a lasso-like post-processing. In this work, we develop a Bayesian approach for automatic determination of the number of dynamic modes using a sparsity-promoting prior. This approach works without manual tuning of hyperparameters through the empirical Bayes technique.

Following [Park and Casella, 2008], we incorporate the two-level Laplacian prior on w1:k, replacing prior (8) and hyperprior (9) by

with new hyperparameters γ1:k. They change the parameters of the conditional distributions for w1:k and σ2 (at Steps 1

and 5 in Algorithm 2) as follows:

Further, the distribution for v2 1:k,1:n (at Step 2) becomes the generalized inverse Gaussian distribution (see e.g. [Devroye, 2014]) with the following parameters:

To draw a sample from the generalized inverse Gaussian distribution, we used an efficient sampler of [Devroye, 2014].

Empirical Bayes for hyperparameter The set of hyperparameters γ1:k needs to be chosen appropriately for successful model selection. We determine it by maximizing the marginal likelihood, since we empirically found that this was more stable than using gamma distribution as a hyperprior for γ1:k. We use a Monte Carlo EM algorithm [Casella, 2001], which comprises iterations between the Gibbs sampling with the modified parameters (E-step) and the maximization of the marginal likelihood (M-step) by

where γ (Q) i denotes the hyperparameter at the Q-th iteration of the EM, and Eγ (Q-1) i [·] denotes the expectation under the hyperparameter at the previous iteration.

5 RelatedWork

DMD was originally proposed as a tool for diagnostics of fluid flows [Rowley et al., 2009; Schmid, 2010], and it has been utilized in various fields of science and engineering, including fluid mechanics [Schmid et al., 2011], analyses of power systems [Susuki and Mezi´c, 2014], epidemiology [Proctor and Eckhoff, 2015], robotic control [Berger et al., 2015], neuroscience [Brunton et al., 2016a], chaotic systems [Brunton et al., 2016b], image processing [Kutz et al., 2016b], and nonlinear system identification [Mauroy and Goncalves, 2016]. Moreover, there are several algorithmic variants such as the use of nonlinear basis functions [Williams et al., 2015], formulation in a reproducing kernel Hilbert space [Kawahara, 2016], and consideration for controlled systems [Proctor et al., 2016].

While no previous work incorporates the probabilistic and Bayesian point of view to DMD, several studies elaborated on the effects of the observation noise; [Duke et al., 2012] and [Pan et al., 2015] conducted error analyses on the outputs of DMD, and there is a line of research on low-rank approximation of DMD [Chen et al., 2012; Wynn et al., 2013; Jovanovi´c et al., 2014; Dicle et al., 2016; H´eas and Herzet, 2017], with which we can mitigate the noise by ignoring insignificant components of data. In addition, [Dawson et al., 2016] proposed total-least-squares DMD, which explicitly considered the presence of observation noise in datasets by formulating DMD as a total least-squares problem. Note that, in Proposition 1, we have shown that the MLE of probabilistic DMD coincides with the solution of total-least-squares DMD.

6 Numerical Examples

We conducted experiments to demonstrate the performance of Bayesian DMD (termed BDMD in this section) regarding the tolerance to noise, the posterior inference, and the automatic determination of the number of modes. In addition, we examined the applications of BDMD to dimensionality reduction and time-series denoising tasks.

6.1 Estimation with Noisy Observations

We validated the performance of BDMD on two types of noisy datasets: one was obtained from a limit cycle, and the other was generated from a system with damping modes.

Limit cycle We generated data from the discrete-time Stuart-Landau equation in polar-coordinates:

and the noisy observable (i is the imaginary unit here):

where each element of et was sampled independently from zero-mean Gaussian with variance 10-4 . The Stuart-Landau equation contains a limit cycle at r = √µ. We set the parameters by µ = 1, γ = 1, β = 0, ∆t = 0.01, r0 = √µ, and θ0 = 0, generated 10,000 snapshots, and fed them into standard DMD (Algorithm 1), total-least-squares DMD (TLSDMD) [Dawson et al., 2016], and BDMD (with k = 5). The estimated eigenvalues are plotted in Figure 2 wherein the ellipses denote the 95% credible interval of the samples generated from the Gibbs sampler of BDMD, for each eigenvalue. While there is the bias on the estimation by standard DMD due to the observation noise, the estimations by TLS-DMD and BDMD coincide, which agrees with Proposition 1. Note that one of the advantages of BDMD is that it returns the posterior distribution of the parameters, instead of the point estimation like TLS-DMD.

Damping modes We also investigated the performance for identifying damping modes, i.e., modes that decay rapidly over time. Generally, it is more difficult to identify damping modes than to identify modes in a limit cycle. The dataset was generated by

where et was zero-mean Gaussian noise with different variances σ2 (σ = 0, 0.05, 0.1, 0.15, 0.2, 0.25), and we set λ1 = 0.9 and λ2 = 0.8 as the eigenvalues. We compared the performances of standard DMD, TLS-DMD, and BDMD (with k = 2). A typical instance of the results is depicted in Figure 3 wherein the box plots show the statistics of the samples generated from the Gibbs sampler of BDMD. The sample medians of BDMD and the estimations by TLS-DMD lie near, and both are more accurate than the estimations by standard DMD. In addition, we ran 100 trials on the same type of datasets generated with different random seeds. In Table 1, the averages of the absolute errors of estimated eigenvalues are listed. We can observe that the point-estimate performance of BDMD is comparable to that of TLS-DMD.

6.2 Automatic Relevance Determination

We conducted an experiment to investigate how well BDMD can determine the number of modes automatically. We gen-

Figure 5: Averages of the RMSEs of estimated dynamic modes over 100 trials for each number of snapshots.

6

erated a dataset by

where et was zero-mean Gaussian noise with variance 10-4 . For determining the number of modes given noisy datasets, standard DMD (Algorithm 1) may utilize the truncation of small singular values at SVD step (Step 1), but the truncation threshold is not trivial in practice. Sparsity-promoting DMD (SP-DMD) [Jovanovi´c et al., 2014] uses a lasso-like post-processing for automatic determination of the number of modes, but it still requires to tune the regularization parameter. However, BDMD with the sparsity-promoting prior (termed BDMD-sp hereafter) can automatically determine the number of modes and the hyperparameter in the light of data inherently, without any need for manual tuning.

We applied standard DMD, SP-DMD (with γ = 10 tuned to give the best results), and BDMD-sp (with k = 4) to the above-mentioned data. As for BDMD-sp, we adopted the medians of the samples generated from the Gibbs sampler as the point estimation values. A typical instance of the results is depicted in Figure 4 wherein the structures of the true modes (matrix A) and the estimated modes are shown. In this case, SP-DMD and BDMD-sp successfully recover the structure of dynamic modes. Furthermore, we ran 100 trials with the same type of datasets generated with different random seeds, varying the number of snapshots fed into the algorithms from m = 4 to 9. We investigated root-mean-square errors (RMSEs) between the estimated and the true modes, which were calculated after normalizing the maximum absolute values and sorting the order of the modes. The results are summarized in Figure 5 wherein the averages (and the standard deviations) of the RMSEs are plotted. We can see that BDMD-sp achieves smaller errors than SP-DMD does.

6.3 Applications

We show two examples of BDMD applications: the dimensionality reduction and the time-series denoising.

Dimensionality reduction

BDMD-sp provides a way for dimensionality reduction of time-series data, since it can concentrate their information on a small number of dynamic modes. To demonstrate the performance, we address the task of data visualization using BDMD-sp on the motion capture data of human activities. 1 We chose locomotion data of three subjects (Subjects #2, #16 and #35), for which both "walk" and "run/jog" motions were recorded. We concatenated the recordings of "walk" and "run/jog" of the three subjects and subsampled them by 1/4, finally obtaining 62-dimensional 421 measurements.

The results of the dimensionality reduction by PCA, t-SNE [van der Maaten and Hinton, 2008], and BDMD-sp (with k = 32) are plotted in Figure 6. As for PCA, we plot only the first and the second principal scores, since the characteristics of the first eight principal scores were all similar. As for BDMD-sp, we focus on latent variable ϕ corresponding to the dynamic mode of the largest magnitude and plot the medians of its samples generated from the Gibbs sampler. Now let us elaborate on the features of the results in Figure 6. The distinction between "walk" and "run/jog" is clearly observed as the different distributions of the trajectories in every plot of Figure 6. The distinction between Subjects #2, #16 and #35 is less obvious, while the distribution of the trajectories implies the difference of the locomotive behavior of Subject #2 from those of the other two. On this point, BDMD-sp (Figure 6c) shows the most consistent result wherein the trajectories of Subject #2 are consistently distributed in the upper part of the plot.

1 Downloaded from http://mocap.cs.cmu.edu/.

Univariate time-series denoising

We prepared a time-series dataset by extracting single series {x} from the Lorenz attractor [Lorenz, 1963] (with ρ = 28, σ = 10 and β = 8/3) and contaminated them with zeromean Gaussian noise of variance 16. The task was recovering the original series from the noisy series. We applied BDMD (with k = 1) on the noisy series and reconstructed them using samples generated by the Gibbs sampler.

The original and the reconstructed series are plotted in Figure 7. The RMSE decreased from 3.2 to 2.3 by the denoising. A simple moving average as a baseline achieved RMSE 2.5 at the best, but note that we cannot necessarily obtain such performance by moving average since it needs to tune the window size.

7 Conclusions

We have introduced the probabilistic model corresponding to DMD and based on that model, proposed Bayesian DMD to conduct posterior inference on the DMD parameters and to enrich the DMD techniques systematically in the unified Bayesian framework. We have shown that the MLE of the proposed probabilistic model coincides with the solution of the standard DMD algorithm in the no-noise limit. Moreover, we have provided the Gibbs sampler for the posterior inference in Bayesian DMD. We have also discussed the case of using the sparsity-promoting prior for automatic determination of the effective number of dynamic modes. Finally, we have presented the results of the experiments with the synthetic and the real-world datasets, which show the effectiveness of Bayesian DMD.

Based on the Bayesian framework proposed in this study, there would be various possible extensions of DMD. One of the promising extensions would be the use of structured priors on dynamic modes. For example, the dynamic modes modeled with Markov random fields fit for images, and applications in natural language processing are possible with discrete probability distributions as prior. Then a challenge would be an efficient inference; we relied on the simple Gibbs sampler in this study, but developing more fast and efficient ways is of great importance.

Acknowledgements

This work was supported by JSPS KAKENHI Grant Numbers JP15J09172, JP16H01548, JP26280086, JP15K12639, and JP26289320.