Sparsity-Promoting Dynamic Mode Decomposition

dmdsp-visual 

Mihailo R. Jovanovic, Peter J. Schmid, Joseph W. Nichols

September 2013

Matlab Files
Download all Matlab files

Presentation
2013 SIAM Talk

Paper
2014 Physics of Fluids Paper

Purpose

This website provides a Matlab implementation of the Sparsity-Promoting Dynamic Mode Decomposition (DMDSP) algorithm. Dynamic Mode Decomposition (DMD) is an effective means for capturing the essential features of numerically or experimentally generated snapshots, and its sparsity-promoting variant DMDSP achieves a desirable tradeoff between the quality of approximation (in the least-squares sense) and the number of modes that are used to approximate available data. Sparsity is induced by augmenting the least-squares deviation between the matrix of snapshots and the linear combination of DMD modes with an additional term that penalizes the ell_1-norm of the vector of DMD amplitudes. We employ alternating direction method of multipliers (ADMM) to solve the resulting convex optimization problem and to efficiently compute the globally optimal solution.

Problem formulation

Dynamic Mode Decomposition

Starting with a sequence of snapshots

 left{   psi_0, psi_1, ldots, psi_N right}, 	~~ 	psi_i in {bf C}^M

DMD provides a low-dimensional representation of a discrete-time linear time-invariant system

 psi_{t + 1} 	; = ; 	A , psi_t, 	~~ 	t 	, = , 	left{   0, ldots, N - 1 	right}

on a subspace spanned by the basis resulting from the Proper Orthogonal Decomposion (POD) of the data sequence. In particular, given two data matrices

 begin{array}{rcl} Psi_0 & = & left[   begin{array}{cccc}   psi_0   &   psi_1   &   cdots   &   psi_{N-1}   end{array} right] , in ; {bf C}^{M times N} [0.25cm] Psi_1 & = & left[   begin{array}{cccc}   psi_1   &   psi_2   &   cdots   &   psi_{N}   end{array} right] , in ; {bf C}^{M times N} end{array}

the DMD algorithm provides optimal representation F in {bf C}^{r times r} of the matrix A in the basis spanned by the POD modes of Psi_0,

     A     ; = ;     U , F , U^*.

Here, r is the rank of the matrix of snapshots Psi_0, and U^* is the complex-conjugate-transpose of the matrix of POD modes U which is obtained from an economy-size singular value decomposition (SVD) of Psi_0,

 	Psi_0 ; = ; U , Sigma , V^*

where Sigma is an r times r diagonal matrix determined by the non-zero singular values { sigma_1, ldots, sigma_r } of Psi_0, and

     begin{array}{rcl}     U , in , {bf C}^{M times r}     &     mbox{with}     &     U^* , U     ; = ;     I     [0.25cm]     V , in , {bf C}^{r times N}     &     mbox{with}     &     V^* , V     ; = ;     I.     end{array}

The matrix F is determined from the matrices of snapshots Psi_0 and Psi_1 by minimizing the Frobenius norm of the difference between Psi_1 and A , Psi_0, thereby resulting into

 	F_{mathrm{dmd}}     ; = ;     U^* , Psi_1 , V , Sigma^{-1}

Optimal amplitudes of DMD modes

The matrix F_{mathrm{dmd}} in {bf C}^{r times r} determines an optimal low-dimensional representation of the inter-snapshot mapping A in {bf C}^{M times M} on the subspace spanned by the POD modes of Psi_0. The dynamics on this r-dimensional subspace are governed by

     x_{t + 1}     ; = ;     F_{mathrm{dmd}} , x_t

and the matrix of POD modes U can be used to map x_t into a higher dimensional space {bf C}^M,

   psi_t ; approx ; U , x_t.

If F_{mathrm{dmd}} has a full set of linearly independent eigenvectors { y_1, ldots, y_r }, with corresponding eigenvalues { mu_1, ldots, mu_r }, then it can be brought into a diagonal coordinate form and experimental or numerical snapshots can be approximated using a linear combination of DMD modes,

     psi_t     ; approx ;     sum_{i , = , 1}^{r}     ,     phi_i     ,     mu_i^t     ,     alpha_i,     ~~      t     ; in ,     left{     0, ldots, N - 1     right}

or, equivalently, in matrix form,

     overbrace{     left[     begin{array}{cccc}     psi_0     &     psi_1     &     cdots     &     psi_{N-1}     end{array}     right]     }^{Psi_0}     ~ approx ~     overbrace{     left[     begin{array}{cccc}     phi_1     &     phi_2     &     cdots     &     phi_{r}     end{array}     right]     }^{Phi}     ,     overbrace{     left[     begin{array}{cccc}     alpha_1     &     &     &          &     alpha_2     &     &          &     &     ddots     &          &     &     &     alpha_r     end{array}     right]     }^{D_{alpha} ; mathrel{mathop:}= ~ hbox{diag} , { alpha }}     ,     overbrace{     left[     begin{array}{cccc}     1     &     mu_1     &     cdots     &     mu_1^{N-1}     [0.15cm]     1     &     mu_2     &     cdots     &     mu_2^{N-1}          vdots     &     vdots     &     ddots     &     vdots     [0.15cm]     1     &     mu_r     &     cdots     &     mu_r^{N-1}     end{array}     right]     }^{V_{mathrm{and}}}.

Here, the matrix of POD modes U and the matrix of eigenvectors of F_{mathrm{dmd}},  Y  mathrel{mathop:}=  left[ begin{array}{ccc} y_1 & cdots & y_r end{array} right], are used to determine the matrix of DMD modes

     Phi     ; = ,     left[     begin{array}{ccc}     phi_1 & cdots & phi_r     end{array}     right]     , = ;     U , Y     , in , {bf C}^{M times r}.

Furthermore, the amplitude alpha_i quantifies the ith modal contribution of the initial condition x_0 on the subspace spanned by the POD modes of Psi_0, and the temporal evolution of the dynamic modes is governed by the Vandermonde matrix V_{mathrm{and}} in {bf C}^{r times N} which is determined by the r eigenvalues mu_i of F_{mathrm{dmd}}.

dmdsp-visual 

Dynamic mode decomposition can be used to represent experimentally or numerically generated snapshots as a linear combination of DMD modes, properly weighted by their amplitudes and advanced in time according to their temporal growth/decay rate.

Determination of the optimal vector of amplitudes  alpha mathrel{mathop:}=  left[  begin{array}{ccc}  alpha_1 & cdots & alpha_r  end{array}  right]^T then amounts to minimization of the Frobenius norm of the difference between Psi_0 and Phi , D_{alpha} , V_{mathrm{and}}. Using SVD of Psi_0, the definition of the matrix Phi, and a sequence of straightforward algebraic manipulations, this convex optimization problem can be brought into the following form

     begin{array}{rl}     {rm minimize}     &     J (alpha)     ; mathrel{mathop:}= ;     alpha^* P , alpha     ; - ;     q^* alpha     ; - ;     alpha^* q     ; + ;     s     end{array}

where alpha is the optimization variable and the problem data is determined by

     P     ; = ,     left(     Y^* , Y     right)     circ     left(     overline{V_{mathrm{and}} , V_{mathrm{and}}^*}     right),     ~~     q     ; = ;     overline{{rm diag} left( V_{mathrm{and}} , V , Sigma^* , Y right)},     ~~     s     ; = ;     {rm trace}     left(     Sigma^* Sigma     right).

Here, an asterisk denotes the complex-conjugate-transpose of a vector (matrix), an overline signifies the complex-conjugate of a vector (matrix), {rm diag} of a vector is a diagonal matrix with its main diagonal determined by the elements of a given vector, {rm diag} of a matrix is a vector determined by the main diagonal of a given matrix, and circ is the elementwise multiplication of two matrices.

The optimal vector of DMD amplitudes that minimizes the Frobenius norm of the difference between Psi_0 and Phi , D_{alpha} , V_{mathrm{and}} can thus be obtained by minimizing the quadratic function J (alpha) with respect to alpha,

     alpha_{mathrm{dmd}}     ; = ;     P^{-1} q     ; = ,     left(     left(     Y^* , Y     right)     circ     left(     overline{V_{mathrm{and}} , V_{mathrm{and}}^*}     right)     right)^{-1}     ,     overline{{rm diag} left( V_{mathrm{and}} , V , Sigma^* , Y right)}

Sparsity-promoting dynamic mode decomposition

We now direct our attention to the problem of selecting the subset of DMD modes that has the most profound influence on the quality of approximation of a given sequence of snapshots. In the first step, we seek a sparsity structure that achieves a user-defined tradeoff between the number of extracted modes and the approximation error (with respect to the full data sequence). In the second step, we fix the sparsity structure of the vector of amplitudes (identified in the first step) and determine the optimal values of the non-zero amplitudes.

dmdsp-visual 

The sparsity-promoting dynamic mode decomposition is aimed at identifying a subset of DMD modes that optimally approximate the entire data sequence.

We approach the problem of inducing sparsity by augmenting the objective function J (alpha) with an additional term that penalizes the the ell_1-norm of the vector of unknown amplitudes alpha,

     begin{array}{rl}     {rm minimize}     &     J (alpha)     ; + ;     gamma     ,     displaystyle{sum_{i , = , 1}^{r}} , | alpha_i |     end{array}     hspace{1.cm}     {rm (SP)}

In the modified optimization problem (SP), gamma is a positive regularization parameter that reflects our emphasis on sparsity of the vector alpha. Larger values of gamma place stronger emphasis on the number of non-zero elements in the vector alpha (relative to the quality of the least-squares approximation, J (alpha)), thereby encouraging sparser solutions.

J-vs-gamma 

Increased emphasis on sparsity encourages sparser solutions at the expense of compromising quality of least-squares approximation.

After a desired balance between the quality of approximation of experimental or numerical snapshots and the number of DMD modes is achieved, we fix the sparsity structure of the unknown vector of amplitudes and determine only the non-zero amplitudes as the solution to the following constrained convex optimization problem:

 	begin{array}{cc}     begin{array}{rl}     {rm minimize}     &     J (alpha)     [0.25cm]     {rm subject~to}     &     E^T , alpha     , = ,     0     end{array}     &     hspace{1.cm}     {rm (POL)}     end{array}

Here, the matrix E in {bf R}^{r times m} encodes information about the sparsity structure of the vector alpha. The columns of E are the unit vectors in {bf R}^{r} whose non-zero elements correspond to zero components of alpha. For example, for alpha in {bf C}^4 with

     alpha     ; = ,     left[     begin{array}{cccc}     alpha_1     &     0     &     alpha_3     &     0     end{array}     right]^T

the matrix E is given as

     E     ; = ,     left[     begin{array}{cc}     0 & 0          1 & 0          0 & 0          0 & 1     end{array}     right].

Alternating direction method of multipliers (ADMM)

ADMM is a simple but powerful algorithm well-suited to large optimization problems. In the sparsity-promoting DMD problem, the algorithm consists of four steps:

  • Step 1: introduce additional variable/constraint

     begin{array}{ll}     {rm minimize}     &     J(alpha)     ; + ;     gamma ,  g (beta)     [0.25cm]     {rm subject~to}     &     alpha ; - ; beta ; = ; 0     end{array}
  • Step 2: introduce the augmented Lagrangian

     {cal L}_rho (alpha,beta,lambda)     ~ = ~     J(alpha)     ~ + ~     gamma     ,     g(beta)     ~ + ~     displaystyle{ frac{1}{2} }     left(     lambda^*     left(     alpha     ; - ;     beta     right)     ~ + ~     left(     alpha     ; - ;     beta     right)^*     lambda     ~ + ~     rho     ,     |     alpha     ; - ;     beta     |_2^2     right)
  • Step 3: use ADMM for the augmented Lagrangian minimization

     begin{array}{rllll}     bf{alpha !! - !! minimization~problem:}         &&     alpha^{k+1}     & !! mathrel{mathop:}= !! &     {displaystyle {rm arg , min}_{alpha}}     ;     {cal L}_rho     (alpha,beta^k,lambda^k)     [0.35cm]     bf{beta !! - !! minimization~problem:}             &&     beta^{k+1}     & !! mathrel{mathop:}= !! &     {displaystyle {rm arg , min}_{beta}}     ;     {cal L}_rho     (alpha^{k+1},beta,lambda^k)     [0.35cm]     bf{lambda !! - !! update~step:}     &&         lambda^{k+1}     & !! mathrel{mathop:}= !! &     lambda^{k}     ;+;     rho left( alpha^{k+1} ; - ; beta^{k+1} right)     end{array}
  • Step 4: polishing - solve the structured quadratic programming problem for the identified sparsity pattern

     {bf structured~optimization~problem}!!:     ~~     left{     begin{array}{ll}     {rm minimize}     &     J(alpha)     [0.2cm]     {rm subject~to}     &     E^T , alpha     , = ,     0     end{array}     right.

Solution to (SP)

The respective structures of the functions J and g in the sparsity-promoting DMD problem can be exploited to show that the alpha-minimization step amounts to solving an unconstrained regularized quadratic program and that the beta-minimization step amounts to a use of the soft thresholding operator S_kappa (cdot):

     begin{array}{rllll}     bf{alpha !! - !! minimization~problem:}         &&     alpha^{k+1}     & !! = !! &     left(     P , + , (rho/2) , I     right)^{-1}     left(      q      , + ,      (rho/2)      left(       beta^k     ; - ;     ( 1/rho )     ,     lambda^k     right)     right)     [0.35cm]     bf{beta !! - !! minimization~problem:}             &&     beta_i^{k+1}     & !! = !! &     S_{kappa}     (       alpha_i^{k+1}     ; + ,     left( 1/rho right)     lambda_i^k ),     ~~     kappa     ; = ;     gamma/rho,     ~~     i      ; = ;      left{     1, ldots, r     right}     [0.35cm]     bf{lambda !! - !! update~step:}     &&         lambda^{k+1}     & !! = !! &     lambda^{k}     ;+;     rho left( alpha^{k+1} ; - ; beta^{k+1} right)     end{array}

where

     S_kappa (v_i^k)     ; mathrel{mathop:}= ;     left{     begin{array}{ll}     v_i^k ,-, kappa,     &     v_i^k     , > ,     kappa     [0.15cm]     0,     &     v_i^k     , in ,     [ -kappa, , kappa ]     [0.15cm]     v_i^k ,+, kappa,     &     v_i^k     , < ,     -kappa     end{array}     right.

Solution to (POL)

After the desired sparsity structure has been identified, the optimal vector of amplitudes with a fixed sparsity structure, alpha_{mathrm{pol}}, can be determined from:

    alpha_{mathrm{pol}}    ; = ;    left[    begin{array}{cc}    {I} & {0}    end{array}    right]    left[    begin{array}{cc}     {P} & {E}    [0.05cm]    {E^T} & {0}    end{array}    right]^{-1}    left[    begin{array}{c}    {q}    [0.1cm]    {0}    end{array}    right]

Acknowledgments

We gratefully acknowledge Prof. Parviz Moin for his encouragement to pursue this effort during the 2012 Center for Turbulence Research Summer Program at Stanford University.