OOPTIMAL LOW-RANK DYNAMIC MODE DECOMPOSITION
Patrick H´eas and C´edric Herzet
INRIA Centre Rennes - Bretagne Atlantique, Campus universitaire de Beaulieu, 35000 Rennes, France
ABSTRACT
Dynamic Mode Decomposition (DMD) has emerged as apowerful tool for analyzing the dynamics of non-linear sys-tems from experimental datasets. Recently, several attemptshave extended DMD to the context of low-rank approxima-tions. This extension is of particular interest for reduced-order modeling in various applicative domains, e.g. , forclimate prediction, to study molecular dynamics or micro-electromechanical devices. This low-rank extension takes theform of a non-convex optimization problem. To the best ofour knowledge, only sub-optimal algorithms have been pro-posed in the literature to compute the solution of this problem.In this paper, we prove that there exists a closed-form optimalsolution to this problem and design an effective algorithm tocompute it based on Singular Value Decomposition (SVD).A toy-example illustrates the gain in performance of theproposed algorithm compared to state-of-the-art techniques.
Index Terms — Low-Rank Approximations, Reduced-Order Models, Dynamical Mode Decomposition, SVD
1. INTRODUCTION
In many fields of Sciences, one is interested in studying thespatio-temporal evolution of a state variable characterized bya partial differential equation. Numerical discretization inspace and time leads to a high dimensional system of equa-tions of the form: (cid:40) x t = f t ( x t − ) ,x = θ, (1)where each element of the sequence of state variables { x t } t belongs to R n , f t : R n → R n with the initial condition θ ∈ R n . Because (1) may correspond to a very high-dimensionalsystem in some applications, computing a trajectory { x t } t given an initial condition θ may lead to a heavy computa-tional load, which may prohibit the direct use of the originalhigh-dimensional system.The context of uncertainty quantification provides an ap-pealing example. Assume we are interested in characterizingthe distribution of random trajectories generated by (1) withrespect to the distribution of the initial condition. A straight-forward approach would be to sample the initial condition and run the high-dimensional system. However, in many applica-tive contexts, it is impossible to generate enough trajecto-ries to make accurate approximations with Monte-Carlo tech-niques.As a response to this computational bottleneck, reduced-order models aim to approximate the trajectories of the sys-tem for a range of regimes determined by a set of initial con-ditions [1]. A common approach is to assume that the tra-jectories of interest are well approximated in a sub-space of R n . In this spirit, many tractable low-rank approximationsof high-dimensional systems have been proposed in the lit-erature, the most familiar being proper orthogonal decom-position (POD) [2], balanced truncation [3], Taylor expan-sions [4] or reduced-basis techniques [5]. Other popular sub-space methods, such as linear inverse modeling (LIM) [6],principal oscillating patterns (POP) [7], or more recently, dy-namic mode decomposition (DMD) [8, 9, 10, 11, 12], areknown as Koopman operator approximations.In this paper, we consider the setting where system (1)is a black-box. In other words, we assume that we do notknow the exact form of f t in (1) and we only have accessto a set of representative trajectories { x it } t,i , i = 1 , ..., N , t = 1 , ..., T so-called snapshots , obtained by running thehigh-dimensional system for N different initial conditions.Moreover, we focus on the low-rank DMD approximationproblem studied in [9, 10]. In a nutshell, these studies pro-vide a procedure for determining a matrix ˆ A k ∈ R n × n of rank k (cid:28) n , which substitutes for function f t in (1) as (cid:40) ˜ x t = ˆ A k ˜ x t − , ˜ x = θ, (2)and generates the approximations ˜ x t ∈ R n with a low com-putational effort. Alternatively, given ˆ A k , and its k non-zeroeigenvalues λ i ∈ C , i = 1 · · · k and associated eigenvectors φ i ∈ C n , i = 1 · · · k , trajectories of (2) can be computed byusing the reduced-order model ˜ x t = k (cid:88) i =1 ν i,t φ i , ν i,t = λ t − i φ ∗ i θ, (3)as long as matrix ˆ A k is symmetric. We will assume it is al-ways the case for simplification issues. In what follows, wewill refer to the parameters φ i and ν i,t ∈ C as the i -th low-rank DMD mode and amplitude at time t . a r X i v : . [ s t a t . M L ] M a y atrix ˆ A k targets the solution of the following non-convex optimization problem, which we will refer to as the low-rank DMD approximation problem A (cid:63)k ∈ arg min A : rank ( A ) ≤ k (cid:88) t,i (cid:107) x it − Ax it − (cid:107) , (4)where (cid:107) · (cid:107) refers to the (cid:96) norm. In order to compute asolution A (cid:63)k , the authors in [9, 10] propose to rely on the as-sumption of linear dependence of recent snapshots on previ-ous ones. This assumption may not be reasonable, especiallyin the case of non-linear systems.Beyond the reduced modeling context discussed above,there has been a resurgence of interest for low-rank solutionsof linear matrix equations [13]. This class of problems is verylarge and includes in particular problem (4). Problems in thisclass are generally nonconvex and do not admit explicit so-lutions. Howewer, important results have arisen at the theo-retical and algorithmic level, enabling the characterization ofthe solution for this class of problems by convex relaxation[14]. Applications concern scenarios such as low-rank ma-trix completion, image compression or minimum order linearsystem realization, see [13]. Nevertheless, there exists certaininstances with a very special structure, which admit closed-form solutions [15, 16]. This occurs typically when the solu-tion can be deduced from the well-known Eckart-Young the-orem [17].The contribution of this paper is to show that the specialstructure of problem (4) enables the characterization of an ex-act closed-form solution and an easily implementable solverbased on singular value decomposition (SVD). In the case k ≥ N ( T − , the proposed algorithm computes the solutionof [12]. More interestingly, for k < N ( T − , i.e. , in the con-strained case, our approach enables to solve exactly the low-rank DMD approximation problem without 1) any assump-tion of linear dependence, 2) the use of an iterative solver, onthe contrary to the approaches proposed in [9, 10, 14].The paper is organized as follows. In section 2, we pro-vide a brief review of state-of-the-art techniques to computelow-rank DMD of experimental data. Section 3 details our an-alytical solution and the algorithm solving (4). Given this op-timal solution, it then presents the reduced-order model solv-ing (2). Finally, a numerical evaluation of the method is pre-sented in Section 4 and concluding remarks are given in a lastsection.
2. STATE-OF-THE-ART OVERVIEW
In what follows, we assume that we have at our disposal N trajectories of T snapshots. We will need in the followingsome matrix notations. The symbol (cid:107) · (cid:107) F and the upperscript · ∗ will respectively refer to the Frobenius norm and thetranspose operator. I k will denote the k -dimensional identitymatrix. Let consecutive elements of the i -th snapshot trajec-tory between time t and t be gathered in a matrix X it : t = ( x it , · · · , x it ) , and let two large matrices X , Y ∈ R n × m with m = ( T − N be defined as X = ( X T − , ..., X N T − ) , Y = ( X T , ..., X N T ) . (5)Without loss of generality, this work will assume that m ≤ n and that rank ( X ) = rank ( Y ) = m . We introduce theSVD decomposition of a matrix M ∈ R p × q with p ≥ q : M = W M Σ M V ∗ M with W M ∈ R p × q , V M ∈ R q × q and Σ M ∈ R q × q so that W ∗ M W M = V ∗ M V M = I q and Σ M isdiagonal. The Moore-Penrose pseudo-inverse of a matrix M will be defined as M † = V M Σ − M W ∗ M .With these notations, problem (4) can be rewritten as A (cid:63)k ∈ arg min A : rank ( A ) ≤ k (cid:107) Y − A X (cid:107) F . (6)In what follows, we begin by presenting two state-of-the-art methods which enable to compute an approximation of thesolution of problem (6). As detailed herafter, the original DMD approach first pro-posed in [8], so-called projected
DMD in [12], assumes thatcolumns of A X are in the span of X . The assumption is writ-ten by the authors in [8, 10] as the existence of A c ∈ R m × m ,the so-called companion matrix of A parametrized by m co-efficients, such that A X = X A c . (7)We remark that this assumption is in particular valid when the i -th snapshot x iT can be expressed as a linear combination ofthe columns of X i T − and when f t is linear. Using the SVDdecomposition X = W X Σ X V ∗ X and noticing X is full rank,we obtain from (7) a projected representation of A in the basisspanned by the columns of W X , W ∗ X AW X = ˜ A c , (8)where ˜ A c = Σ X V ∗ X A c V X Σ − X ∈ R m × m . Therefore, the low-rank formulation in [10] proposes to approach the solutionof (6) by determining the m coefficients of matrix A c whichminimize the Frobenius norm of the residual Y − A X . Thisyields after some algebraic manipulations to solve the prob-lem arg min ˜ A c : rank ( ˜ A c Σ X ) ≤ k (cid:107) W ∗ X Y V X − ˜ A c Σ X (cid:107) F . (9)The Eckart-Young theorem [17] then provides the optimal so-lution to this problem based on a rank- k SVD approximationof matrix B = W ∗ X Y V X given by W B Λ B V ∗ B where Λ B is adiagonal matrix containing only the k -largest singular valuesf Σ B and with zero entries otherwise. Exploiting the low-dimensional representation (8), a reduced-order model for tra-jectories can then be obtained by inserting in (2) the low-rankapproximation ˆ A k = W X W B Λ B V ∗ B Σ − X W ∗ X . (10)As an alternative, the authors propose a reduced-order modelfor trajectories relying on the so-called DMD modes and theiramplitudes. These modes are related to the eigenvectors ofthe solution of (9). The amplitudes are given by solving aconvex optimization problem with an iterative gradient-basedmethod, see details in [10]. If we remove the low-rank constraint, (6) becomes a least-squares problem whose solution is ˆ A m = YX † = Y V X Σ − X W ∗ X . (11)Based on the approximation ˆ A m , DMD modes and ampli-tudes serve to design a model to reconstruct trajectories of(2). We note that the DMD modes are simply given by theeigendecomposition of ˆ A m , which can be efficiently com-puted using SVD, as proposed in [12]. The associated DMDamplitudes can then easily be derived.It is important to remark that truncating to a rank- k the so-lution of the above unconstrained minimization problem willnot necessarily yield the solution of (6). This approach willgenerally be sub-optimal. However surprisingly, the solutionto problem (6) remain to our knowledge overlooked in the lit-erature, and no algorithms enabling non-projected low-rankDMD approximations have yet been proposed.
3. THE PROPOSED APPROACH3.1. Closed-form Solution to (6)Let the columns of matrix P ∈ R n × k be the real orthonormaleigenvectors associated to the k largest eigenvalues of matrix YY ∗ . Theorem 1
A solution of (6) is A (cid:63)k = P P ∗ YX † . This theorem states that (6) can be simply solved by comput-ing the orthogonal projection of the unconstrained problemsolution (11) onto the subspace spanned by the k first eigen-vectors of YY ∗ . A detailed proof is provided in the technicalreport associated to this paper [18]. The matrix YY ∗ is of size n × n . Since n is typically verylarge, this prohibits the direct computation of an eigenvaluedecomposition. The following well-know remark is useful toovercome this difficulty. Algorithm 1
Solver for (6) input : N -sample { X i T } Ni =1
1) Form matrix X and Y as defined in (5).2) Compute the SVD of X .3) Compute the columns of P using Remark 1. output : matrix A (cid:63)k = P P ∗ Y V X Σ − X W ∗ X Algorithm 2
Low-rank DMD modes and amplitudes input : matrices ( P, Q, θ ) , with Q = ( YX † ) ∗ P
1) Compute the SVD of matrix Q .2) Solve for i = 1 · · · k the eigen equation ˜ A k w i = λ i w i , where w i ∈ C m and λ i ∈ C denote eigenvectors and eigen-values of ˜ A k = W ∗ Q P V Q Σ Q ∈ R m × m . output : DMD modes φ i = W Q w i and amplitudes ν i,t = λ t − i φ ∗ i θ Remark 1
The eigenvectors associated to the m ≤ n non-zero eigenvalues of matrix YY ∗ ∈ R n × n with Y ∈ R n × m can be obtained from the eigenvectors V Y = ( v , ..., v m ) ∈ R m × m and eigenvalues of the smaller matrix Y ∗ Y ∈ R m × m . Indeed, the SVD of a matrix Y of rank m is Y = W Y Σ Y V ∗ Y , where the columns of matrix W Y ∈ R n × m are the eigenvectors of YY ∗ . Since V Y is unitary, we obtainthat the sought vectors are the first k columns of W Y , i.e., of Y V Y Σ − Y . In the light of this remark, it is straightforward to designAlgorithm 1, which will compute efficiently the solution of(6) based on SVDs.
We now discuss the resolution of the reduced-order model(2) given the solution A (cid:63)k of (6). Trajectories of (2) are fullydetermined by a k -dimensional recursion involving the pro-jected variable z t = P ∗ ˜ x t : (cid:40) z t = P ∗ YX † P z t − ,z = P ∗ YX † θ. (12)Then, by multiplying both sides by matrix P , we obtain thesought low-rank approximation ˜ x t = P z t .Alternatively, we can employ reduced-order model (3).The parameters of this model, i.e. , low-rank DMD modes andamplitudes, are efficiently computed without any minimiza-tion procedure, in contrast to what is proposed by the authorin [10]. Indeed, we rely on the following remark stating thatDMD modes and amplitudes can be obtained by means ofSVDs using Algorithm 2. The remark is proved in the techni-cal report [18]. Remark 2
Each pair ( φ i , λ i ) generated by Algorithm 2 isone of the k eigenvector/eigenvalue pair of A (cid:63)k . ank E rr o r no r m f o r s e tt i ng i ) method a)method b)method c) Rank E rr o r no r m f o r s e tt i ng ii ) method a)method b)method c)Rank E rr o r no r m f o r s e tt i ng iii ) method a)method b)method c) Fig. 1 . Evaluation of error norms (cid:107) Y − ˆ A k X (cid:107) F as a function ofrank k . Setting i ) (top) and ii ) (middle) imply both a linear modelbut the former satisfies the snapshots linear dependence assumption.Setting iii ) (bottom) implements a non-linear model. We evaluate3 algorithms: method a ) is the proposed optimal algorithm, method b ) provides the rank- k SVD approximation of the unconstrained so-lution given in [12] and method c ) is the low-rank projected DMDmethod proposed in [10]. See details in Section 4.
4. NUMERICAL EVALUATION
In what follows, we evaluate on a toy model the different ap-proaches for solving the low-rank DMD approximation prob-lem. We consider a high-dimensional space of n = 50 dimen-sions, a low-dimensional subspace of r = 30 dimensions and m = 40 snapshots. Let G be a matrix of rank r generated ran-domly according to G = (cid:80) ri =1 ξ i ξ ∗ i , where entries of ξ i ’s are n independent samples of the standard normal distribution.Let the initial condition θ be randomly chosen according tothe same distribution. The snapshots, gathered in matrices X and Y , are generated using (1) for three configurations of f t : i ) f t ( x t − ) = Gx t − , s.t. ∃ A c satisfying G X = X A c , ii ) f t ( x t − ) = Gx t − , iii ) f t ( x t − ) = Gx t − + G diag ( x t − ) diag ( x t − ) x t − .Setting i ) corresponds to a linear system satisfying the as-sumption (7), as made in the projected DMD approaches[8, 10]. Setting ii ) and iii ) , do not make this assumptionand simulate respectively linear and non-linear dynamicalsystems. We assess three different methods for computing ˆ A k : a ) optimal rank- k approximation given by Algorithm 1, b ) k th-order SVD approximation of (11), i.e. , k -th orderapproximation of the rank- m non-projected DMD so-lution [12], c ) rank- k approximation by (10), corresponding to theprojected DMD approach [10] (or [8] for k ≥ m ).The performance is measured in terms of the error norm (cid:107) Y − ˆ A k X (cid:107) F with respect to the rank k . Results for the three set-tings are displayed in Figure 1.As a first remark, we notice that the solution provided byAlgorithm 1 (method a ) yields the best results, in agreementwith Theorem 1.Second, in setting i ) , the experiments confirm that whenthe linearity assumption is valid, the low-rank projected DMD(method c ) achieves the same performance as the optimal so-lution (method a ). Moreover, truncating the rank- m DMDsolution (method b ) induces as expected an increase of theerror norm. This deterioration is however moderate in ourexperiments.Then, in settings ii ) and iii ) we remark that the behaviorof the error norms are analogous (up to an order of magni-tude). The performance of the projected approach (method c )differs notably from the optimal solution. A significant de-terioration is visible for k > . This is the consequence ofthe non-validity of the assumption made in method c . Never-theless, we notice that method c accomplishes a slight gain inperformance compared to method b up to a moderate rank( k < ). Besides, we also notice that the error norm ofmethod b in the case k < is not optimal.Finally, as expected, all methods succeed in properly char-acterizing the low-dimensional subspace as soon as k ≥ r .
5. CONCLUSION
Following recent attempts to characterize an optimal low-rankapproximation based on DMD, this paper provides a closed-form solution to this non-convex optimization problem. Tothe best of our knowledge, state-of-the-art methods are allsub-optimal. The paper further proposes effective algorithmsbased on SVD to solve this problem and run reduced-ordermodels. Our numerical experiments attest that the proposedalgorithm is more accurate than state-of-the-art methods. Inparticular, we illustrate the fact that simply truncating the full-rank DMD solution, or exploiting too restrictive assumptionsfor the approximation subspace is insufficient. cknowledgements
This work was supported by the “Agence Nationale de laRecherche” through the GERONIMO project (ANR-13-JS03-0002).
A. REFERENCES [1] A. Cohen and R. Devore, “Approximation of high-dimensional parametric PDEs,”
ArXiv e-prints , Feb.2015.[2] P. Holmes, J. L. Lumley, and G. Berkooz,
Turbulence,coherent structures, dynamical systems and symmetry ,Cambridge University Press, 1996, Cambridge BooksOnline.[3] A. C. Antoulas, “An overview of approximation meth-ods for large-scale dynamical systems,”
Annual Reviewsin Control , vol. 29, no. 2, pp. 181–190, Jan. 2005.[4] J. P. Fink and W. C. Rheinboldt, “On the error behav-ior of the reduced basis technique for nonlinear finiteelement approximations,”
ZAMM - Journal of AppliedMathematics and Mechanics , vol. 63, no. 1, pp. 21–28,1983.[5] A. Quarteroni, G. Rozza, and A. Manzoni, “Certified re-duced basis approximation for parametrized partial dif-ferential equations and applications,”
Journal of Math-ematics in Industry , vol. 1, no. 1, pp. 1–49, Dec. 2011.[6] C. Penland and T. Magorian, “Prediction of nino 3sea surface temperatures using linear inverse modeling,”
Journal of Climate , vol. 6, no. 6, pp. 1067–1076, 1993.[7] K. Hasselmann, “PIPs and POPs: The reduction of com-plex dynamical systems using principal interaction andoscillation patterns,”
Journal of Geophysical Research:Atmospheres , vol. 93, no. D9, pp. 11015–11021, 1988.[8] P. J. Schmid, “Dynamic mode decomposition of numeri-cal and experimental data,”
Journal of Fluid Mechanics ,vol. 656, pp. 5–28, 2010.[9] K. K. Chen, J. H. Tu, and C. W. Rowley, “Variantsof dynamic mode decomposition: boundary condition,koopman, and fourier analyses,”
Journal of NonlinearScience , vol. 22, no. 6, pp. 887–915, 2012.[10] MR Jovanovic, PJ Schmid, and JW Nichols, “Low-rankand sparse dynamic mode decomposition,”
Center forTurbulence Research Annual Research Briefs , pp. 139–152, 2012.[11] M. O. Williams, I.G Kevrekidis, and C.W. Rowley, “Adata–driven approximation of the koopman operator:extending dynamic mode decomposition,”
Journal ofNonlinear Science , vol. 25, no. 6, pp. 1307–1346, 2015. [12] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brun-ton, and J. N. Kutz, “On dynamic mode decomposition:theory and applications,”
Journal of Computational Dy-namics , vol. 1, no. 2, pp. 391–421, 2014.[13] B. Recht, M. Fazel, and P. A. Parrilo, “Guaranteedminimum-rank solutions of linear matrix equations vianuclear norm minimization,”
SIAM review , vol. 52, no.3, pp. 471–501, 2010.[14] M. Fazel,
Matrix rank minimization with applications,Stanford University , Ph.D. thesis, 2002.[15] P. A. Parrilo and S. Khatri, “On cone-invariant linearmatrix inequalities,”
IEEE Transactions on AutomaticControl , vol. 45, no. 8, pp. 1558–1563, 2000.[16] M. Mesbahi and G. P Papavassilopoulos, “On the rankminimization problem over a positive semidefinite lin-ear matrix inequality,”
IEEE Transactions on AutomaticControl , vol. 42, no. 2, pp. 239–243, 1997.[17] C. Eckart and G. Young, “The approximation of onematrix by another of lower rank,”
Psychometrika , vol.1, no. 3, pp. 211–218, 1936.[18] P. H´eas and C. Herzet, “Low-rank Approximation andDynamic Mode Decomposition,”