---
abstract: 'Dynamic Mode Decomposition (DMD) has emerged as a powerful tool for analyzing the dynamics of non-linear systems from experimental datasets. Recently, several attempts have extended DMD to the context of low-rank approximations. This extension is of particular interest for reduced-order modeling in various applicative domains, [*e.g.*, ]{}for climate prediction, to study molecular dynamics or micro-electromechanical devices. This low-rank extension takes the form of a non-convex optimization problem. To the best of our knowledge, only sub-optimal algorithms have been proposed in the literature to compute the solution of this problem. In this paper, we prove that there exists a closed-form optimal solution to this problem and design an effective algorithm to compute it based on Singular Value Decomposition (SVD). A toy-example illustrates the gain in performance of the proposed algorithm compared to state-of-the-art techniques.'
address: 'INRIA Centre Rennes - Bretagne Atlantique, Campus universitaire de Beaulieu, 35000 Rennes, France'
bibliography:
- './bibtex.bib'
title: 'Optimal Low-Rank Dynamic Mode Decomposition'
---

Low-Rank Approximations, Reduced-Order Models, Dynamical Mode Decomposition, SVD

Introduction
============

In many fields of Sciences, one is interested in studying the [spatio-temporal]{} evolution of a state variable characterized by a partial differential equation. Numerical discretization in space and time leads to a high dimensional system of equations [of the form:]{} $$\begin{aligned}
\label{eq:model_init} 
 \left\{\begin{aligned}
& x_{t}= f_t(x_{t-1}) , \\&x_1={\theta},
\end{aligned}\right. \vspace{-0.cm}$$ [where]{} each element of the sequence of state variables $\{x_t\}_t$ belongs to ${\mathds{R}}^n$, $f_t:{\mathds{R}}^n \to {\mathds{R}}^n$ with the initial condition $\theta \in {\mathds{R}}^n$. [Because may correspond to a very high-dimensional system in some applications, computing a trajectory $\{x_t\}_t$ given an initial condition ${\theta}$]{} may lead to a heavy computational load, which may prohibit the direct use of the original high-dimensional system.

The context of uncertainty quantification provides an appealing example. Assume we are interested in characterizing the distribution of random trajectories generated by with respect to the distribution of the initial condition. A straightforward approach would be to sample the initial condition and run the high-dimensional system. However, in many applicative contexts, it is impossible to generate enough trajectories to make accurate approximations with Monte-Carlo techniques.

As a response to this computational bottleneck, reduced-order models aim to approximate the trajectories of the system for a range of regimes determined by a set of initial conditions [@2015arXiv150206797C]. A common approach is to assume that the trajectories of interest are well approximated in a sub-space of ${\mathds{R}}^n$. In this spirit, many tractable low-rank approximations of high-dimensional systems have been proposed in the literature, the most familiar being proper orthogonal decomposition (POD) [@9780511622700], balanced truncation [@Antoulas2005Overview], Taylor expansions [@ZAMM:ZAMM19830630105] or reduced-basis techniques [@Quarteroni2011Certified]. Other popular sub-space methods, such as linear inverse modeling (LIM) [@penland1993prediction], principal oscillating patterns (POP) [@Hasselmann88], or more recently, dynamic mode decomposition (DMD) [@Schmid10; @Chen12; @Jovanovic12; @williams2015data; @Tu2014391], are known as Koopman operator approximations.

In this paper, we consider the setting where system is a black-box. In other words, we assume that we do not know the exact form of $f_t$ in and we only have access to a set of representative trajectories $\{x^i_t\}_{t,i}$, $i=1,...,N$, $t=1,...,T$ so-called *snapshots*, obtained by running the high-dimensional system for $N$ different initial conditions. Moreover, we focus on the low-rank DMD approximation problem studied in [@Chen12; @Jovanovic12]. In a nutshell, these studies provide a procedure for determining a matrix $ \hat A_k{\in{\mathds{R}}^{n\times n}}$ of rank $k \ll n$, which substitutes [for]{} function $f_t$ in $$\begin{aligned}
\label{eq:model_koopman_approx} 
 \left\{\begin{aligned}
& \tilde x_{t}=  \hat A_k \tilde  x_{t-1} , \\&\tilde x_1={\theta},
\end{aligned}\right. \vspace{-0.cm}$$ and generates the approximations $\tilde x_{t}\in {\mathds{R}}^n$ with a low computational effort. Alternatively, given $\hat A_k$, and its $k$ non-zero eigenvalues $\lambda_i \in {\mathds{C}},\, i=1 \cdots k$ and associated eigenvectors $\phi_i \in {\mathds{C}}^n,\, i=1 \cdots k$, trajectories of can be computed by using the reduced-order model $$\begin{aligned}
\tilde x_{t}=\sum_{i=1}^k \nu_{i,t} \phi_i,\label{eq:koopman1}\quad \nu_{i,t} =  \lambda_i^{t-1} \phi_i^*\theta,\end{aligned}$$ as long as matrix $\hat A_k$ is symmetric. We will assume it is always the case for simplification issues. In what follows, we will refer to the parameters $\phi_i$ and $\nu_{i,t}\in {\mathds{C}}$ as the *$i$-th low-rank DMD mode and amplitude* at time $t$.

Matrix $\hat A_k$ targets the solution of the following non-convex optimization problem, which we will refer to as the [*low-rank DMD approximation problem*]{} $$\begin{aligned}
\label{eq:prob} 
A_k^\star \in &\operatorname*{arg\,min}_{A:\textrm{rank}(A)\le k}  \sum_{t,i} \| x^i_{t} - A   x^i_{t-1}\|^2_2, $$ where $\|\cdot\|_2$ refers to the $\ell_2$ norm. In order to compute [a]{} solution $A_k^\star$, the authors in [@Chen12; @Jovanovic12] propose to rely on the assumption of linear dependence of recent snapshots on previous ones. This assumption may not be reasonable, especially in the case of non-linear systems.

Beyond the reduced modeling context discussed above, there has been a resurgence of interest for low-rank solutions of linear matrix equations [@recht2010guaranteed]. This class of problems is very large and includes in particular problem . Problems in this class are generally nonconvex and do not admit explicit solutions. Howewer, important results have arisen at the theoretical and algorithmic level, enabling the characterization of the solution for this class of problems by convex relaxation [@fazel2002matrix]. Applications concern scenarios such as low-rank matrix completion, image compression or minimum order linear system realization, see [@recht2010guaranteed]. Nevertheless, there exists certain instances with a very special structure, which admit closed-form solutions [@parrilo2000cone; @mesbahi1997rank]. This occurs typically when the solution can be deduced from the well-known Eckart-Young theorem [@eckart1936approximation].

The contribution of this paper is to show that the special structure of problem enables the characterization of an exact closed-form solution and an easily implementable solver based on singular value decomposition (SVD). In the case $k \ge N(T-1)$, the proposed algorithm computes the solution of [@Tu2014391]. More interestingly, for $k < N(T-1)$, [*i.e.*, ]{}in the constrained case, our approach enables to solve exactly the low-rank DMD approximation problem without 1) any assumption of linear dependence, 2) the use of an iterative solver, on the contrary to the approaches proposed in [@Chen12; @Jovanovic12; @fazel2002matrix]. The paper is organized as follows. In section \[sec:stateArt\], we provide a brief review of state-of-the-art techniques to compute low-rank DMD of experimental data. Section \[sec:contrib\] details our analytical solution and the algorithm solving . Given this optimal solution, it then presents the reduced-order model solving . Finally, a numerical evaluation of the method is presented in [Section]{} \[sec:numEval\] and concluding remarks are [given]{} in a last section.

State-Of-The-Art Overview {#sec:stateArt}
=========================

In what follows, we assume that we have at our disposal $N$ trajectories of $T$ snapshots. We will need in the following some matrix notations. The symbol $\|\cdot\|_F$ and the upper script $\cdot^*$ will respectively refer to the Frobenius norm and the transpose operator. $I_k$ will denote the $k$-dimensional identity matrix. Let consecutive elements of the $i$-th snapshot trajectory between time $t_1$ and $t_2$ be gathered in a matrix ${{X}}_{t_1:t_2}^i = (x^i_{t_1},\cdots,x^i_{t_2}),$ and let two large matrices $ {\mathbf{X}}, {\mathbf{Y}}\in {\mathds{R}}^{n \times  m} $ with $m=(T-1)N$ be defined as $$\begin{aligned}
\label{eq:matrixAB}
 {\mathbf{X}}= ({{X}}^1_{1:T-1}, ..., {{X}}^N_{1:T-1}), \quad 
 {\mathbf{Y}}= ( {{X}}^1_{2:T}, ...,  {{X}}^N_{2:T}) .
 \end{aligned}$$ Without loss of generality, this work will assume that $m\leq n$ and that $\textrm{rank}({\mathbf{X}})=\textrm{rank}({\mathbf{Y}})=m$. We introduce the SVD decomposition of a matrix $M\in {\mathds{R}}^{p \times q }$ with $p\ge q$: $M=W_M\Sigma_M V_M^*$ with $W_M\in {\mathds{R}}^{p \times q }$, $V_M\in {\mathds{R}}^{ q \times  q}$ and $\Sigma_M\in {\mathds{R}}^{q  \times q }$ so that $W_M^*W_M=V_M^*V_M=I_q$ and $\Sigma_M$ is diagonal. The Moore-Penrose pseudo-inverse of a matrix $M$ will be defined as $M^\dagger=V_M\Sigma_M^{-1} W_M^*$. With these notations, problem can be rewritten as $$\begin{aligned}
\label{eq:prob1} 
        A_k^\star \in &\operatorname*{arg\,min}_{A:\textrm{rank}(A)\le k} \|{\mathbf{Y}}-A {\mathbf{X}}\|_F^2.
        \end{aligned}$$ In what follows, we begin by presenting two state-of-the-art methods which enable to compute an approximation of the solution of problem .

Projected DMD and Low-Rank Formulation
--------------------------------------

As detailed herafter, the original DMD approach first proposed in [@Schmid10], so-called *projected* DMD in [@Tu2014391], assumes that columns of $A{\mathbf{X}}$ are in the span of ${\mathbf{X}}$. The assumption is written by the authors in [@Schmid10; @Jovanovic12] as the existence of $A^c\in {\mathds{R}}^{ m \times m}$, the so-called *companion matrix* of $A$ parametrized by $m$ coefficients, such that $$\begin{aligned}
A {\mathbf{X}}={\mathbf{X}}A^c.\label{eq:companion}
  \end{aligned}$$ We remark that this assumption is in particular valid when the $i$-th snapshot $x^i_T$ can be expressed as a linear combination of the columns of ${{X}}_{1:T-1}^i $ and when $f_t$ is linear. Using the SVD decomposition ${\mathbf{X}}=W_{\mathbf{X}}\Sigma_{\mathbf{X}}V_{\mathbf{X}}^*$ and noticing ${\mathbf{X}}$ is full rank, we obtain from a projected representation of $A$ in the basis spanned by the columns of $W_{\mathbf{X}}$, $$\begin{aligned}
\label{eq:DMDassumption}
 W_{\mathbf{X}}^*AW_{\mathbf{X}}=\tilde A^c,\end{aligned}$$ where $ \tilde A^c=\Sigma_{\mathbf{X}}V_{\mathbf{X}}^*A^cV_{\mathbf{X}}\Sigma_{\mathbf{X}}^{-1}\in {\mathds{R}}^{ m \times  m}.$ Therefore, the low-rank formulation in [@Jovanovic12] proposes to approach the solution of by determining the $m$ coefficients of matrix $A^c$ which minimize the Frobenius norm of the residual $ {\mathbf{Y}}-A {\mathbf{X}}$. This yields after some algebraic manipulations to solve the problem $$\begin{aligned}
\label{eq:DMDSVD}
\operatorname*{arg\,min}_{ \tilde A^c: \textrm{rank}(\tilde A^c\Sigma_{\mathbf{X}})\le k} \|W_{\mathbf{X}}^*{\mathbf{Y}}V_{\mathbf{X}}- \tilde A^c \Sigma_{\mathbf{X}}\|^2_F.\end{aligned}$$ The Eckart-Young theorem [@eckart1936approximation] then provides the optimal solution to this problem based on a rank-$k$ SVD approximation of matrix $B=W_{\mathbf{X}}^*{\mathbf{Y}}V_{\mathbf{X}}$ given by $W_B\Lambda_BV_B^*$ where $\Lambda_B$ is a diagonal matrix containing only the $k$-largest singular values of $\Sigma_B$ and with zero entries [otherwise]{}. Exploiting the low-dimensional representation , a reduced-order model for trajectories can then be obtained by inserting in the low-rank approximation $$\begin{aligned}
\label{eq:projDMD}
 \hat A_k=W_{\mathbf{X}}W_B\Lambda_BV_B^*\Sigma_{\mathbf{X}}^{-1}W_{\mathbf{X}}^*.
 \end{aligned}$$ As an alternative, the authors propose a reduced-order model for trajectories relying on the so-called DMD modes and their amplitudes. These modes are related to the eigenvectors of the solution of . The amplitudes are given by solving a convex optimization problem with an iterative gradient-based method, see details in [@Jovanovic12].

Non-projected DMD
-----------------

If we remove the low-rank constraint, becomes a least-squares problem whose solution is $$\begin{aligned}
\label{eq:exactDMD}
\hat A_m={\mathbf{Y}}{\mathbf{X}}^{\dagger}={\mathbf{Y}}V_{\mathbf{X}}\Sigma_{\mathbf{X}}^{-1}W_{\mathbf{X}}^*.\end{aligned}$$ Based on the approximation $\hat A_m$, DMD modes and amplitudes serve to design a model to reconstruct trajectories of . We note that the DMD modes are simply given by the eigendecomposition of $\hat A_m$, which can be [efficiently]{} [computed]{} using SVD, as proposed in [@Tu2014391]. The [associated]{} DMD amplitudes can then easily be derived.

It is important to remark that truncating to a rank-$k$ the solution of the above unconstrained minimization problem will not necessarily yield the solution of . This approach will generally be sub-optimal. However surprisingly, the solution to problem remain to our knowledge overlooked in the literature, and no algorithms enabling non-projected low-rank DMD approximations have yet been proposed.

 The Proposed Approach {#sec:contrib}
======================

Closed-form Solution to 
------------------------

Let the columns of matrix ${P}\in {\mathds{R}}^{n \times k}$ be the real orthonormal eigenvectors associated to the $k$ largest eigenvalues of matrix ${\mathbf{Y}}{\mathbf{Y}}^*.$

\[prop22\] A solution of is $A_k^\star={P}{P}^*{\mathbf{Y}}{\mathbf{X}}^{\dagger}.$

This theorem states that can be simply solved by computing the orthogonal projection of the unconstrained problem solution onto the subspace spanned by the $k$ first eigenvectors of ${\mathbf{Y}}{\mathbf{Y}}^*$. A detailed proof is provided in [the technical report associated to this paper [@Heas16_DMD].]{}

[Efficient]{} Solver
--------------------

The matrix ${\mathbf{Y}}{\mathbf{Y}}^*$ is of size $n\times n$. [Since $n$ is typically very large, t]{}his prohibits the direct computation of an eigenvalue decomposition. The following well-know remark is useful to overcome this difficulty.\

\[rem:1\] The eigenvectors associated to the $m~\le~n$ non-zero eigenvalues of matrix ${\mathbf{Y}}{\mathbf{Y}}^*\in {\mathds{R}}^{n \times n}$ with $ {\mathbf{Y}}\in {\mathds{R}}^{n \times m}$ can be obtained [from]{} the eigenvectors $ V_{\mathbf{Y}}=(v_1,...,v_{m}) \in {\mathds{R}}^{m \times m}$ and eigenvalues of the smaller matrix ${\mathbf{Y}}^*{\mathbf{Y}}\in {\mathds{R}}^{m\times m}$. Indeed, the SVD of a matrix $ {\mathbf{Y}}$ of rank $m$ is ${\mathbf{Y}}=  W_{\mathbf{Y}}\Sigma_{\mathbf{Y}}V_{\mathbf{Y}}^*,$ where the columns of matrix $  W_{\mathbf{Y}}\in {\mathds{R}}^{n \times m} $ are the eigenvectors of ${\mathbf{Y}}{\mathbf{Y}}^*$. Since $  V_{\mathbf{Y}}$ is unitary, we obtain that the sought vectors are the [first $k$]{} columns of $ W_{\mathbf{Y}}$, [*i.e.*, ]{}of $ {\mathbf{Y}}V_{\mathbf{Y}}\,\Sigma_{\mathbf{Y}}^{-1}.$

In the light of this remark, it is straightforward to design Algorithm \[algo:1\], which will compute efficiently the solution of based on SVDs.

**input**: $N$-sample $\{{{X}}_{1:T}^i\}_{i=1}^N$1) Form matrix ${\mathbf{X}}$ and ${\mathbf{Y}}$ as defined in . 2) Compute the SVD of ${\mathbf{X}}$. 3) Compute the columns of ${P}$ using Remark \[rem:1\].**output**: matrix $A_k^\star ={P}{P}^*{\mathbf{Y}}V_{\mathbf{X}}\Sigma_{\mathbf{X}}^{-1}W_{\mathbf{X}}^*$

**input**: matrices $({P}, {Q},\theta)$, with ${Q}= ({\mathbf{Y}}{\mathbf{X}}^{\dagger})^* {P}$ 1) Compute the SVD of matrix ${Q}$.2) Solve for $ i=1 \cdots k$ the eigen equation $ \tilde A_k w_i = \lambda_i w_i, $ where $w_i \in {\mathds{C}}^m$ and $\lambda_i \in {\mathds{C}}$ denote eigenvectors and eigenvalues of $\tilde  A_k=W_{{Q}}^*{P}V_{{Q}}\Sigma_{{Q}} \in {\mathds{R}}^{m\times m}.$

**output**: DMD modes $  \phi_i=W_{{Q}} w_i$ and amplitudes $\nu_{i,t} =  \lambda_i^{t-1} \phi_i^*\theta$

Reduced-Order Models {#sec:ROM}
--------------------

We now discuss the resolution of the reduced-order model given the solution $A_k^\star$ of . Trajectories of are fully determined by a $k$-dimensional recursion involving the projected variable $z_t= {P}^* \tilde x_t$: $$\label{eq:GaussLinSystReduced}
 \left\{\begin{aligned}
& z_t={P}^*{\mathbf{Y}}{\mathbf{X}}^{\dagger}{P}z_{t-1} ,\\
&z_2={P}^*{\mathbf{Y}}{\mathbf{X}}^{\dagger}\theta.  \\
\end{aligned}\right.$$ Then, by multiplying both sides by matrix ${P}$, we obtain the sought low-rank approximation $\tilde x_{t} = {P}z_t$. Alternatively, we can employ reduced-order model . The parameters of this model, [*i.e.*, ]{}low-rank DMD modes and amplitudes, are efficiently computed without any minimization procedure, in contrast to what is proposed by the author in [[@Jovanovic12]]{}. Indeed, we rely on the following remark stating that DMD modes and amplitudes can be obtained by means of SVDs using Algorithm \[algo:2\]. The remark is proved in the [technical report [@Heas16_DMD]]{}.

\[rem:2\] Each pair $(\phi_i,\lambda_i)$ generated by Algorithm \[algo:2\] is one of the $k$ eigenvector/eigenvalue pair of $A_k^\star$.

Numerical Evaluation {#sec:numEval}
====================

  ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
    ![Evaluation of error norms $\|{\mathbf{Y}}-\hat A_k {\mathbf{X}}\|_F$ as a function of rank $k$. Setting $i)$ [([top]{})]{} and $ii)$ (middle) imply both a linear model but the former satisfies the snapshots linear dependence assumption. Setting $iii)$ [([bottom]{})]{} implements a non-linear model. We evaluate 3 algorithms: [method $a)$]{} is the proposed optimal algorithm, [method $b)$]{} provides the rank-$k$ SVD approximation of the unconstrained solution given in [@Tu2014391] and [method $c)$]{} is the low-rank projected DMD method proposed in [@Jovanovic12]. See details in Section \[sec:numEval\]. \[fig:1\]](res2.pdf "fig:"){width="0.9\columnwidth"}
    ![Evaluation of error norms $\|{\mathbf{Y}}-\hat A_k {\mathbf{X}}\|_F$ as a function of rank $k$. Setting $i)$ [([top]{})]{} and $ii)$ (middle) imply both a linear model but the former satisfies the snapshots linear dependence assumption. Setting $iii)$ [([bottom]{})]{} implements a non-linear model. We evaluate 3 algorithms: [method $a)$]{} is the proposed optimal algorithm, [method $b)$]{} provides the rank-$k$ SVD approximation of the unconstrained solution given in [@Tu2014391] and [method $c)$]{} is the low-rank projected DMD method proposed in [@Jovanovic12]. See details in Section \[sec:numEval\]. \[fig:1\]](res1.pdf "fig:"){width="0.9\columnwidth"}
   ![Evaluation of error norms $\|{\mathbf{Y}}-\hat A_k {\mathbf{X}}\|_F$ as a function of rank $k$. Setting $i)$ [([top]{})]{} and $ii)$ (middle) imply both a linear model but the former satisfies the snapshots linear dependence assumption. Setting $iii)$ [([bottom]{})]{} implements a non-linear model. We evaluate 3 algorithms: [method $a)$]{} is the proposed optimal algorithm, [method $b)$]{} provides the rank-$k$ SVD approximation of the unconstrained solution given in [@Tu2014391] and [method $c)$]{} is the low-rank projected DMD method proposed in [@Jovanovic12]. See details in Section \[sec:numEval\]. \[fig:1\]](res3_NL.pdf "fig:"){width="0.9\columnwidth"}
  ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

In what follows, we evaluate on a toy model the different approaches for solving the low-rank DMD approximation problem. We consider a high-dimensional space of $n=50$ dimensions, a low-dimensional subspace of $r=30$ dimensions and $m=40$ snapshots. Let $G$ be a matrix of rank $r$ generated randomly according to $G=\sum_{i=1}^{r} \xi_i \xi_i^*$, where entries of $\xi_i$’s are $n$ independent samples of the standard normal distribution. Let the initial condition $\theta$ be randomly chosen according to the same distribution. The snapshots, gathered in matrices ${\mathbf{X}}$ and ${\mathbf{Y}}$, are generated using for three configurations of $f_t$:

-   $f_t(x_{t-1})=G x_{t-1}$, s.t. $\exists A^c$ satisfying $G {\mathbf{X}}={\mathbf{X}}A^c$,

-   $f_t(x_{t-1})=G x_{t-1}$,

-   $f_t(x_{t-1})=G x_{t-1}+G \textrm{diag}(x_{t-1})\textrm{diag}(x_{t-1})x_{t-1}$.

Setting $i)$ corresponds to a linear system satisfying the assumption , as made in the projected DMD approaches [@Schmid10; @Jovanovic12]. Setting $ii)$ and $iii)$, do not make this assumption and simulate respectively linear and non-linear dynamical systems. We assess three different methods for computing $\hat A_k$:

-   optimal rank-$k$ approximation given by Algorithm \[algo:1\],

-   $k$th-order SVD approximation of , [*i.e.*, ]{}$k$-th order approximation of the rank-$m$ non-projected DMD solution [@Tu2014391],

-   rank-$k$ approximation by , corresponding to the projected DMD approach [@Jovanovic12] (or [@Schmid10] for $k\ge m$).

The performance is measured in terms of the error norm $\|{\mathbf{Y}}-\hat A_k {\mathbf{X}}\|_F$ with respect to the rank $k$. Results for the three settings are displayed in Figure \[fig:1\]. As a first remark, we notice that the solution provided by Algorithm \[algo:1\] (method $a$) yields the best results, in agreement with [T]{}heorem \[prop22\].

Second, in setting $i)$, the experiments confirm that when the linearity assumption is valid, the low-rank projected DMD (method $c$) achieves the same performance as the optimal solution (method $a$). Moreover, truncating the rank-$m$ DMD solution (method $b$) induces as expected an increase of the error norm. This deterioration is however moderate in our experiments.

Then, in settings $ii)$ and $iii)$ we remark that the behavior of the error norms are analogous (up to an order of magnitude). The [performance]{} of the projected approach (method $c$) differs notably from the optimal solution. A significant deterioration is visible for $k>10$. This is the consequence of the non-validity of the assumption made in method $c$. Nevertheless, we notice that method $c$ accomplishes a slight gain in performance compared to method $b$ up to a moderate rank ($k<5$). Besides, we also notice that the error norm of method $b$ in the case $k<30$ is not optimal.

Finally, as expected, all methods succeed in properly characterizing the low-dimensional subspace as soon as $k \ge r$.

Conclusion
==========

Following recent attempts to characterize an optimal low-rank approximation based on DMD, this paper provides a closed-form solution to this non-convex optimization problem. To the best of our knowledge, state-of-the-art methods are all sub-optimal. The paper further proposes effective algorithms based on SVD to solve this problem and run reduced-order models. Our numerical experiments attest that the proposed algorithm is more accurate than state-of-the-art methods. In particular, we illustrate the fact that simply truncating the full-rank DMD solution, or exploiting too restrictive assumptions for the approximation subspace is insufficient.

Acknowledgements {#acknowledgements .unnumbered}
================

This work was supported by the “Agence Nationale de la Recherche" through the GERONIMO project (ANR-13-JS03-0002).
