Reduced-order modeling using Dynamic Mode Decomposition and Least Angle Regression

TL;DR

This paper introduces LARS4DMD, a novel method combining Dynamic Mode Decomposition with Least Angle Regression to efficiently select modes for reduced-order modeling without needing a regularization parameter.

Contribution

The paper presents LARS4DMD, a new algorithm that integrates LARS with DMD for automatic mode selection, simplifying the process and maintaining performance.

Findings

LARS4DMD achieves comparable accuracy to DMDSP.

LARS4DMD does not require a regularization parameter.

Numerical tests on Poiseuille flow validate the method.

Abstract

Dynamic Mode Decomposition (DMD) yields a linear, approximate model of a system's dynamics that is built from data. We seek to reduce the order of this model by identifying a reduced set of modes that best fit the output. We adopt a model selection algorithm from statistics and machine learning known as Least Angle Regression (LARS). We modify LARS to be complex-valued and utilize LARS to select DMD modes. We refer to the resulting algorithm as Least Angle Regression for Dynamic Mode Decomposition (LARS4DMD). Sparsity-Promoting Dynamic Mode Decomposition (DMDSP), a popular mode-selection algorithm, serves as a benchmark for comparison. Numerical results from a Poiseuille flow test problem show that LARS4DMD yields reduced-order models that have comparable performance to DMDSP. LARS4DMD has the added benefit that the regularization weighting parameter required for DMDSP is not needed.

Equations

$$\underset{A}{\min }||{\Psi }_1-A{\Psi }_0|{|}_{\mathcal{F}}^2,$$

$$A\approx U_{r}F{U}_r^H.$$

$$F_{\text{DMD}}=U_r^H{\Psi }_1{V}_r{\Sigma }_r^{-1}.$$

$${\mathbf{\varphi }}_j=U_r{\mathbf{v}}_j.$$

$${\mathbf{\psi }}_k\approx \sum _{j=1}^r{\mathbf{\varphi }}_j{\lambda }_j^k{\alpha }_j,$$

$$\displaystyle\Phi=\left[\begin{array}[]{cccc}\bm{\phi}_{1}&\bm{\phi}_{2}&\dots&\bm{\phi}_{r}\end{array}\right],$$

$$\displaystyle\Xi=\left[\begin{array}[]{cccc}\lambda_{1}^{0}&\lambda_{1}^{1}&\dots&\lambda_{1}^{N-1}\\ \lambda_{2}^{0}&\lambda_{2}^{1}&\dots&\lambda_{2}^{N-1}\\ \vdots&\vdots&\ddots&\vdots\\ \lambda_{r}^{0}&\lambda_{r}^{1}&\dots&\lambda_{r}^{N-1}\\ \end{array}\right]$$

$${\Psi }_0^{\text{rec}}=\Phi D_{\alpha }\Xi .$$

$$\underset{\mathbf{\alpha }}{\min }J(\mathbf{\alpha })=||{\Psi }_0-\Phi D_{\alpha }\Xi |{|}_{\mathcal{F}}^2.$$

$$P_{\text{loss}}=100\times \frac{||{\Psi }_0-{\Psi }_0^{\text{rec}}|{|}_{\mathcal{F}}}{||{\Psi }_0|{|}_{\mathcal{F}}}.$$

$$\mathbf{\mu }=\sum _{j=1}^r{\mathbf{x}}_j{\alpha }_j.$$

$$\hat{\mathbf{c}}=X^T(\mathbf{y}-{\hat{\mathbf{\mu }}}_S).$$

$$L_S$$

$${\mathbf{w}}_S$$

$${\mathbf{u}}_S$$

$$\mathbf{g}=X^T{\mathbf{u}}_S.$$

$$\hat{\gamma}=\left\{\begin{array}[]{ll}\min^{+}_{j\in S^{c}}\quad\left\{\frac{\hat{C}-\hat{c}_{j}}{L_{S}-g_{j}},\frac{\hat{C}+\hat{c}_{j}}{L_{S}+g_{j}}\right\}&\text{if $S^{c}$ is nonempty,}\\ \frac{\tilde{C}}{L_{S}}&\text{if $S^{c}$ is empty,}\end{array}\right.$$

$${\hat{\mathbf{\mu }}}_{S,k}={\hat{\mathbf{\mu }}}_{S,k-1}+\hat{\gamma }{\mathbf{u}}_S.$$

$$\alpha_{k,j}=\left\{\begin{array}[]{ll}\alpha_{k-1,j}+\hat{\gamma}s_{j}w_{S,j}&\text{if }j\in S,\\ 0&\text{if }j\notin S\end{array}\right.$$

$$X_S^T{\mathbf{u}}_S=L_S\text{mathbbm}1,$$

$$L_S={\left({\text{mathbbm}1}^T{\left(X_S^T{X}_S\right)}^{-1}\text{mathbbm}1\right)}^{-1/2},$$

$${\mathbf{w}}_S={\left(X_S^T{X}_S\right)}^{-1}{L}_S\text{mathbbm}1,$$

$${\mathbf{u}}_S=X_S{\mathbf{w}}_S.$$

$$⟨\cdot ,\cdot ⟩:\text{mathbbm}{C}^n\times \text{mathbbm}{C}^n⟶\text{mathbbm}C,$$

$$⟨\mathbf{x},\mathbf{y}⟩={\mathbf{x}}^H\mathbf{y},$$

$$\hat{\mathbf{c}}=X^H(\mathbf{y}-{\hat{\mathbf{\mu }}}_S).$$

$$L_S$$

$${\mathbf{w}}_S$$

$${\mathbf{u}}_S$$

$$\mathbf{g}=X^H{\mathbf{u}}_S.$$

$$\hat{\gamma}=\left\{\begin{array}[]{ll}\min^{+}_{j\in S^{c}}\quad\frac{\left(\Re\left\{\left<g_{j},\hat{c}_{j}\right>\right\}-\hat{C}L_{S}\right)\pm\sqrt{\left(\Re\left\{\left<g_{j},\hat{c}_{j}\right>\right\}-\hat{C}L_{S}\right)^{2}-\left(|g_{j}|^{2}-L_{S}^{2}\right)\left(|\hat{c}_{j}|^{2}-\hat{C}^{2}\right)}}{|g_{j}|^{2}-L_{S}^{2}}&\text{if $S^{c}$ is nonempty,}\\ \frac{\tilde{C}}{L_{S}}&\text{if $S^{c}$ is empty,}\end{array}\right.$$

$${\hat{\mathbf{\mu }}}_{S,k}={\hat{\mathbf{\mu }}}_{S,k-1}+\hat{\gamma }{\mathbf{u}}_S.$$

Full text

1 Introduction

Data-driven system analysis has become an increasingly popular technique for studying features of dynamical systems. Dynamic Mode Decomposition (DMD) is a data analysis algorithm that can identify dynamical features that appear in the output data of a system [1]. DMD works by finding the best-fit linear operator $A$ that marches the measurements forward in time. The best-fit linear operator can be analyzed through modal analysis and its eigenvectors are known as the DMD modes of the system with corresponding DMD eigenvalues. These modes can be utilized to give insight into features of dynamical systems such as coherent structures in fluid flows [2] and patterns in neural recordings [3]. It has also been shown that there is a strong connection between DMD and Koopman operator theory [4].

DMD first appeared in the fluids literature [1] and is particularly useful for data-driven analysis of fluid systems [5, 6, 7]. For fluid systems, DMD is frequently used to analyze Particle Image Velocimetry (PIV) data, which is a measurement technique that creates vector-field images of a flowfield in time. Measurement techniques such as PIV often require very fine spatial resolution to resolve the fluid dynamics [8]. When the number of spatial data points is much larger than the total number of temporal observations or snapshots, the number of snapshots often determines the number of DMD modes. Including more snapshots can help to capture the system’s dynamics in the DMD modes. However, since the number of modes increases with the number of snapshots, often the DMD model is large. To deal with a large number of DMD modes, one typically reduces the size of the model by retaining only relatively important modes. Determining the DMD modes that are most important is a reduced-order modeling problem and is the focus of this paper.

A reduced-order modeling algorithm that has been effective in creating sparse DMD models is Sparsity-Promoting Dynamic Mode Decomposition (DMDSP) [9, 10, 11]. DMDSP has also been extended to systems with inputs [12]. It has been used to cancel DMD modes that are associated with noise [13] and for model selection in filtering applications [14]. The authors of [9] observe that the reduced-order modeling problem amounts to appropriately weighting DMD modes in the model, because a mode’s contribution can be neglected by scaling its weight to zero. One process for weighting modes involves optimizing the reconstruction of the original data using the model. The reconstruction problem seeks to minimize the reconstruction error, which is the difference between the full data set and the data predicted by the model. DMDSP creates a reduced-order model through the multi-objective optimization of minimizing the reconstruction error and the $l_{1}$ norm of the vector of mode amplitudes. The $l_{1}$ regularization encourages sparsity in the mode amplitude vector by scaling the mode amplitudes to zero for modes that have minimal influence on the reconstruction error. The final step in DMDSP removes the regularization term and re-solves for the mode amplitudes for the selected set of modes.

The number of modes that DMDSP deactivates depends on the relative weighting between the reconstruction error and the $l_{1}$ regularization term in the optimization. The relative weighting is mediated by a user-specified weighting coefficient. Unfortunately, the user does not know values of this parameter a priori that will produce useful model reconstructions. Further, the user does not know how the weighting parameter will influence the induced sparsity and the associated tradeoff in performance. The user therefore searches for a range of regularization values for which DMDSP provides model reduction. The user then examines the resulting models and selects the one that produces an acceptable tradeoff between performance and model size.

In this work, we propose a new method for determining the DMD modes that are most relevant to the system dynamics through adaption of the Least Angle Regression (LARS) algorithm [15]. LARS is a regression algorithm that sequentially selects vectors from a set of candidate vectors and appropriately scales them to fit a vector of data. These user-specified candidate vectors for model construction are called covariates. As an example, if the measurement data is a disease diagnosis in a set of patients, the covariates may consist of other relevant data collected about the patients such as their age, height, weight, and so forth. LARS determines which of these covariates are important and scales the selected covariates appropriately to yield an estimate of the diagnosis.

LARS is a member of a class of algorithms known as forward selection algorithms that build predictive models in a stepwise manner. The "S" in the LARS acronym references the popular Stagewise and Lasso algorithms that are shown to be variants of LARS [15]. These algorithms sequentially build up a model by traveling in the space of covariates along directions that are determined by the covariates that are most correlated with the residual (i.e., the difference between the data and the current model’s prediction of the data) [15]. The direction of travel that LARS selects is the equiangular direction between the most correlated covariates [15]. We choose to apply a forward selection algorithm for DMD mode selection to create a principled approach for examining the tradeoff between performance and model size during reduced-order model construction. Using a forward selection algorithm, a user can more easily track the influence of an individual mode on the model.

LARS functions by sequentially adding covariates to an active set based on their correlation with the residual. After each selection step, LARS generates an interim estimate of the data that is formed by traveling along an equiangular direction with all covariates in the active set. The contribution of each covariate to travel along the equiangular direction determines its regression coefficient, which weights the covariate in the model.

This paper contributes a modified version of the LARS algorithm, known as Complex LARS, that handles complex-valued data and complex-valued covariates, and the LARS4DMD algorithm that applies Complex LARS for reduced-order DMD model construction. The LARS algorithm was originally developed for data fitting in statistics using real-valued data and real-valued covariates. We adapt LARS to operate using complex data and complex covariates by replacing the dot product with a complex-valued inner product over the vector space $\mathbbm{C}^{n}$. These contributions are significant because they enhance automation in constructing reduced-order DMD models by eliminating the need for a problem-dependent regularization weighting parameter. The performance of the LARS4DMD algorithm is demonstrated on synthesized Poiseuille flow data that is available from the original DMDSP paper [9, 16].

The outline of this paper is as follows. Section 2 presents necessary background on the DMD, DMDSP, and LARS algorithms. Section 3 derives Complex LARS. Section 4 applies Complex LARS for use in DMD mode selection. Section 5 shows the effectiveness of LARS4DMD in generating reduced-order DMD models using DMDSP as a benchmark. Section 6 concludes the paper and discusses ongoing work.

2 Data-driven, reduced-order modeling

This section presents the techniques necessary for development of LARS4DMD: Section 2.1 reviews DMD; Section 2.2 describes DMDSP, a state-of-the-art method for DMD mode selection; and Section 2.3 introduces the original LARS algorithm.

2.1 Dynamic Mode Decomposition (DMD)

The DMD data analysis begins with the collection and proper arrangement of measurements for processing. Although generalized definitions of DMD exist (e.g., see [8, 17]), we focus on the case of sequential, constant-interval measurements of a process evolving in time, similar to [9]. Let $\bm{\psi}_{k}$ be a measurement vector (or snapshot) of the system for time steps $k=0,\dots,N$. DMD seeks the best-fit linear operator $A$ that advances each snapshot one time step such that $\bm{\psi}_{k+1}\approx A\bm{\psi}_{k}$. By constructing two data matrices, \Psi_{0}=[\begin{array}[]{cccc}\bm{\psi}_{0}&\bm{\psi}_{1}&\dots&\bm{\psi}_{(N-1)}\end{array}] and \Psi_{1}=[\begin{array}[]{cccc}\bm{\psi}_{1}&\bm{\psi}_{2}&\dots&\bm{\psi}_{N}\end{array}], that have column entries offset by one timestep, this condition can be expressed as $\Psi_{1}\approx A\Psi_{0}$ [9]. The DMD optimization problem [9]

[TABLE]

where ${\cal F}$ is the Frobenius norm, provides the best-fit $A$ matrix $A^{\ast}=\Psi_{1}\Psi_{0}^{\dagger}$, where $(\cdot)^{\dagger}$ is the Moore-Penrose pseudo-inverse [9].

The eigenvectors of $A$ are the DMD modes with associated DMD eigenvalues. In practice, the $A$ matrix can often be too large to form [1]. When this occurs, the DMD modes and eigenvalues are still accessible by first solving for a projected version of $A$. Let $\Psi_{0}=U\Sigma V^{H}$ be an economy Singular Value Decomposition (SVD) of $\Psi_{0}$, where $(\cdot)^{H}$ denotes the Hermitian or conjugate-transpose operation, and let $r$ be the rank of $\Psi_{0}$. To address the possibility of linearly dependant snapshots in $\Psi_{0}$, truncate $U$, $\Sigma$, and $V$ according to $r$. Let $U_{r}$ be the first $r$ columns of $U$, let $\Sigma_{r}$ be an $r\times r$ matrix extracted from the upper left corner of $\Sigma$, and let $V_{r}$ be the first $r$ columns of $V$.

Consider a version $F$ of the $A$ matrix that is transformed into the basis formed by the columns of $U_{r}$ such that [9]

[TABLE]

Inserting (2) and the SVD of $\Psi_{0}$ into (1) and optimization yields the best-fit projected version of the A matrix, given by [9]

[TABLE]

Performing modal analysis on $F_{\text{DMD}}$ can provide insight on the DMD modes of $A$. Let $\bm{v}_{j}$ be and eigenvector of $F_{\text{DMD}}$ such that $F_{\text{DMD}}\bm{v}_{j}=\lambda_{j}\bm{v}_{j}$ for eigenvalue $\lambda_{j}$. A DMD mode $\bm{\phi}_{j}$ of $A$ can be recovered from the corresponding eigenvector $\bm{v}_{j}$ in the subspace spanned by the columns of $U_{r}$ by [9]

[TABLE]

Using the DMD modes, it is possible to approximate the dynamics of the measurements or snapshots. A measurement vector $\bm{\psi}_{k}$ at time $k$ can be expressed as a linear combination of the DMD modes. The DMD modes evolve in time by repeated multiplication with their corresponding DMD eigenvalues, yielding the snapshot dynamics [9]

[TABLE]

where $\alpha_{j}$ is, in general, a complex scalar that corresponds to the contribution of the mode $\bm{\phi}_{j}$ to the initial snapshot matrix $\Psi_{0}$. The amplitudes $\alpha_{j}$ have also been shown to be equivalent to the values of the Koopman eigenfunctions calculated at the initial condition [4].

Using the snapshot dynamics (5), it is possible to reconstruct the data matrix $\Psi_{0}$ from an initial snapshot vector by letting $\bm{\psi}_{k}^{\text{rec}}=\bm{\psi}_{k}$ for each timestep $k$. The DMD mode amplitudes are computed using $\alpha_{j}=\bm{v}_{j}^{H}U_{r}^{H}\bm{\psi}_{0}$. Solving for the DMD mode amplitudes to best-fit DMD modes to a data set is referred to as the reconstruction problem. In the context of DMD, reduced-order modeling seeks to identify a subset of the DMD modes that perform well in data reconstruction for a dataset or a variety of datasets.

2.2 Sparsity-Promoting Dynamic Mode Decomposition (DMDSP)

DMDSP is a reduced-order modeling technique that selects modes by increasing the sparsity in the vector of DMD amplitudes $\bm{\alpha}$ during a reconstruction optimization problem. To state the reconstruction problem mathematically, put the DMD modes $\bm{\phi}_{j}$ for $j=1,\dots,r$ into matrix form

[TABLE]

Using the DMD eigenvalues $\lambda_{j}$, construct a Vandermonde matrix

[TABLE]

which expresses the temporal evolution of the eigenvalue portion of each mode’s coefficient in the snapshot dynamics (5). Form a diagonal matrix that is constructed from the DMD mode amplitudes such that D_{\alpha}=\text{diag}([\begin{array}[]{cccc}\alpha_{1}&\alpha_{2}&\dots&\alpha_{r}\end{array}]). Using the diagonal amplitude matrix $D_{\alpha}$, the DMD mode matrix (7), and the Vandermonde matrix (12), we can reconstruct $\Psi_{0}$ using the snapshot dynamics [9]

[TABLE]

The reconstruction problem can be expressed as

[TABLE]

DMDSP is a two-step procedure that minimizes $J(\bm{\alpha})+\beta||\bm{\alpha}||_{1}$, where $||\bm{\alpha}||_{1}$ is an $l_{1}$ regularization penalty that serves as a proxy for penalizing the number of nonzero entries in $\bm{\alpha}$, and $\beta$ is a user-defined, regularization weighting term. The second step in DMDSP is a polishing step in which the algorithm re-solves for the mode amplitudes $\bm{\alpha}$ that best solve the reconstruction problem, but the regularization term is not present and the desired sparsity structure from the first step is strictly enforced. Typically, a user of DMDSP considers an array of regularization weights $\beta$. For each $\beta$ value, DMDSP solves for a mode amplitude vector $\bm{\alpha}$ that contains zeros for modes deactivated by the method and nonzero values for the amplitudes of selected modes. Each $\beta$ value therefore corresponds to a separate reduced-order model. However, it should be noted that it is possible, and common in practice, multiple $\beta$ values to yield the same reduced-order model.

The authors of [9] define the percent performance loss to be a measure of the error in reconstruction, normalized by size of the original data set such that

[TABLE]

The user can calculate the performance loss for each model, examine the performance and model size tradeoff, and select a model appropriate for the application.

2.3 Least Angle Regression (LARS)

Let $X=[\bm{x}_{1}\;\bm{x}_{2}\;\dots\;\bm{x}_{r}]$ be a set of zero-mean, unit-variance, and linearly independent covariates $\bm{x}_{j}$ for $j=1,\;\dots,\;r$, and let $\bm{y}$ represent a zero-mean data vector [15]. The LARS algorithm searches for regression coefficients $\alpha_{j}$ for $j=1,\dots,r$, to construct a linear estimate

[TABLE]

Note that we use $\alpha_{j}$ notation to denote the DMD mode amplitudes and the LARS regression coefficients; these quantities correspond when we apply LARS4DMD for mode selection in Section 4.

Algorithm 1 presents the sequential procedure for selection of covariates by LARS [15]. In each iteration of the algorithm, LARS produces an estimate $\bm{\hat{\mu}}_{S}$ of the data based on currently selected covariates in a set of active covariates. The hat $\hat{(\cdot)}$ notation denotes the current iteration, and the subscript $(\cdot)_{S}$ indicates that the quantity is based on selected covariates in the active set. The current estimate initializes with $\bm{\hat{\mu}}_{S}=0$. At each iteration, the difference between the data $\bm{y}$ and the current estimate $\bm{\hat{\mu}}_{S}$ is the residual $(\bm{y}-\bm{\hat{\mu}}_{S})$. LARS calculates the current correlation between each covariate $\bm{x}_{j}$ and the residual $(\bm{y}-\bm{\hat{\mu}}_{S})$ in Step 1.1. After finding the maximum absolute current correlation $\hat{C}=\max_{j}|\hat{c}_{j}|$ in Step 1.2, LARS selects covariates for which $\hat{C}=|\hat{c}_{j}|$ (or covariates for which this equality holds within a small tolerance). When adding covariate $\bm{x}_{j}$ to the active set, LARS multiplies each selected covariate by the sign of its correlation with the residual $\text{sgn}(\hat{c}_{j})$ so that all covariates in the active set have positive correlations with the residual.

Algorithm 1 (LARS)

Inputs: Zero-mean, real data vector $\bm{y}$, and a set of real, zero-mean, unit-variance covariates $X$.

1. 1.1)

Obtain a vector of correlations with the current residual $(\bm{y}-\bm{\hat{\mu}}_{S})$,

[TABLE] 2. 1.2)

Find the current maximum absolute correlation, $\hat{C}=\max_{j}\left|\hat{c}_{j}\right|$. 3. 1.3)

Form the active set $S=\{j\in 1,...,r\,\,\lvert\,\,|\hat{c}_{j}|=\hat{C}\}$, and collect the signum-aligned covariates $s_{j}\bm{x}_{j}$ of the active set, where $s_{j}=\text{sgn}(\hat{c}_{j})$, within the columns of the signum-aligned covariate matrix X_{S}=[\begin{array}[]{ccc}\dots&s_{j}\bm{x}_{j}&\dots\end{array}]. 4. 1.4)

Create an equiangular direction of travel $\bm{u}_{S}$,

[TABLE] 5. 1.5)

Find correlations with the direction of travel for all covariates

[TABLE] 6. 1.6)

Find the length $\hat{\gamma}$ to travel along $\bm{u}_{S}$

[TABLE]

where $S^{c}$ is the complement of $S$ and $\min^{+}$ indicates the minimum taken over positive values only. 7. 1.7)

Update the estimate of the data,

[TABLE]

(Note: Initialize with $\bm{\hat{\mu}}_{S,0}=0$.) 8. 1.8)

Update the regression coefficients for $j=1,\dots,r$,

[TABLE]

(Note: Initialize with $\bm{\alpha}_{0}=0$.) 9. 1.9)

Repeat Steps 1.1-1.8 until all covariates have zero correlation with the residual or until there are no covariates remaining in $S^{c}$.

Output: Vector of regression coefficients $\bm{\alpha}$.

Using the covariates in the active set, LARS determines an equiangular direction of travel $\bm{u}_{S}$ in which LARS can step to reduce the current residual. To calculate the step direction, the equiangular condition [15]

[TABLE]

where $\bm{\mathbbm{1}}$ is a vector of ones, ensures that the dot products between each aligned covariate in the active set $s_{j}\bm{x}_{j}$ and the direction $\bm{u}_{S}$ have an equal value $L_{S}$. Requiring $\bm{u}_{S}$ to be a unit vector, i.e. $\bm{u}_{S}^{T}\bm{u}_{S}=1$, one can utilize (15) to derive the value [15]

[TABLE]

in the equi-angle condition (15). $L_{S}$ used to calculate a vector of weighting coefficients $\bm{w}_{S}$ for the active covariates and the associated equi-angle direction such that [15]

[TABLE]

and

[TABLE]

Traveling along the equiangular direction reduces the current correlation equally among all covariates in the active set [15]. LARS selects a step size to travel in the equi-angle direction $\bm{u}_{S}$ that is as small as possible until another covariate enters the active set. Step 1.6 chooses the distance $\hat{\gamma}$ to travel along the equiangular direction $\bm{u}_{S}$. Subsequently, Step 1.7 uses $\hat{\gamma}\bm{u}_{S}$ to update the current estimate, and Step 1.8 provides the new regression coefficients. The LARS algorithm repeats until all covariates have zero correlation with the residual or until all candidate covariates have joined the active set. The algorithm returns the regression coefficients needed for model construction.

3 Complex Least Angle Regression

As formulated, LARS selects real covariates to fit real data. Often complex covariates can arise in applications, such as in the DMD mode-selection problem that this paper studies. Complex data can also occur, such as if the components of a planar vector field are stored together in complex format. The example in Section 5 uses complex data.

This section modifies LARS to allow for complex covariates and complex data. The adaption of LARS for complex covariates and complex data requires replacement of the inner product space over which the algorithm evolves. An inner product can be used to describe the angular relationship between two vectors, so its selection is important for the LARS algorithm, which seeks an equiangular direction of travel during covariate selection. Let covariates $\bm{x}_{j}$, for $j=1,\dots,r$, and data vector $\bm{y}$ reside in the vector space $\mathbbm{C}^{n}$. Similar to LARS, let $\bm{y}$ be a zero-mean data vector and let $\bm{x}_{j}$ for $j=1,\dots,r$ be zero-mean, unit-variance555The variance for a vector of complex entries is $\text{var}(\bm{q})=1/n\sum_{j=1}^{n}|q_{j}-\text{mean}(\bm{q})|^{2}$. covariates. Consider an inner product of the form [18]

[TABLE]

with the following properties for $\bm{x},\bm{y},\bm{z}\in\mathbbm{C}^{n}$ and $a\in\mathbbm{C}$ [18]:

1. (i)

positive definiteness

          $\left<\bm{x},\bm{x}\right>\geq 0,\text{ with }\left<\bm{x},\bm{x}\right>=0\text{ if and only if }\bm{x}=0,$ 2. (ii)

conjugate symmetry

          $\left<\bm{x},\bm{y}\right>=\overline{\left<\bm{y},\bm{x}\right>},$ 3. (iii)

right linearity

          $\left<\bm{x},a\bm{y}\right>=a\left<\bm{x},\bm{y}\right>,$

          $\left<\bm{x},\bm{y}+\bm{z}\right>=\left<\bm{x},\bm{y}\right>+\left<\bm{x},\bm{z}\right>,$ 4. (iv)

left-conjugate linearity

          $\left<a\bm{x},\bm{y}\right>=\overline{a}\left<\bm{x},\bm{y}\right>,$

          $\left<\bm{x}+\bm{y},\bm{z}\right>=\left<\bm{x},\bm{z}\right>+\left<\bm{y},\bm{z}\right>.$

Note that Property (iv) is a consequence of Properties (ii) and (iii). Left-conjugate linearity represents a choice of inner product convention that is often called the quantum mechanical definition; right-conjugate linearity is another common choice [18]. We select the Euclidean inner product over $\mathbbm{C}^{n}$ defined by [19]

[TABLE]

to satisfy left-conjugate linearity and Properties (i) through (iii).

Algorithm 2 presents Complex LARS. Each step of Complex LARS closely reflects the corresponding step of the LARS Algorithm. Step 1.1 in the LARS algorithm calculates a vector of correlations $\bm{\hat{c}}$ based on the dot product $\bm{x}_{j}^{T}\left(\bm{y}-\bm{\hat{\mu}}_{S}\right)$ between each covariate $\bm{x}_{j}$, for $j=1,\dots,r$, and the current residual $(\bm{y}-\bm{\hat{\mu}}_{S})$. Substitution of the inner product (18) leaves the first step unchanged except for replacing the transpose operation $(\cdot)^{T}$ with the Hermitian operation $(\cdot)^{H}$. Note that the choice of inner product causes minimal change to the equations in Steps 2.1 and 2.4. Although the inner product (18) yields a complex value, we still refer to the elements of $\bm{\hat{c}}$ as correlations to retain the terminology from [15].

**Algorithm 2 (Complex LARS)

**

Inputs: Zero-mean data vector $\bm{y}$, and a set of zero-mean, unit-variance covariates $X$.

1. 2.1)

Obtain a vector of correlations with the current residual $(\bm{y}-\bm{\hat{\mu}}_{S})$,

[TABLE] 2. 2.2)

Find the current maximum absolute correlation $\hat{C}=\max_{j}\left(|\hat{c}_{j}|\right)$. 3. 2.3)

Form the active set $S=\{j\in 1,...,r\,\,\lvert\,\,|\hat{c}_{j}|=\hat{C}\}$, and collect the signum-aligned covariates $s_{j}\bm{x}_{j}$ of the active set, where $s_{j}=\text{sgn}_{\mathbbm{C}}(\hat{c}_{j})$, within the columns of the signum-aligned covariate matrix X_{S}=[\begin{array}[]{ccc}\dots&s_{j}\bm{x}_{j}&\dots\end{array}]. 4. 2.4)

Create an equiangular direction of travel

[TABLE] 5. 2.5)

Find correlations with the direction of travel for all covariates

[TABLE] 6. 2.6)

Find the length $\hat{\gamma}$ to travel along $\bm{u}_{S}$

[TABLE]

where $S^{c}$ is the complement of $S$ and $\min^{+}$ indicates the minimum taken over positive values only. 7. 2.7)

Update the estimate of the data,

[TABLE]

(Note: Initialize with $\bm{\hat{\mu}}_{S,0}=0$.) 8. 2.8)

Update the regression coefficients for $j=1,\dots,r$,

[TABLE]

(Note: Initialize with $\bm{\alpha}_{0}=0$.) 9. 2.9)

Repeat Steps 1.1-1.8 until all covariates have zero correlation with the residual or until there are no covariates remaining in $S^{c}$.

Output: Vector of regression coefficients $\bm{\alpha}$.

Step 2.2 finds the maximum absolute correlation $\hat{C}$ by examining $|\hat{c}_{j}|$ for $j=1,...,r$. In the original LARS algorithm, the covariates corresponding to the maximum absolute correlation are multiplied by the signs of their correlations with the residual to provide sign-aligned versions of the covariates that are added to the active set. To adapt this step for complex LARS, consider the complex signum function .

[TABLE]

which is the complex extension of the real-valued sign function $\text{sgn}(x)=x/|x|$. The $\text{sgn}_{\mathbbm{C}}(\cdot)$ function returns a complex number that lies on the unit circle. Note that if $z$ is strictly real, $\text{sgn}_{\mathbbm{C}}(z)$ agrees with the real-valued sign function $\text{sgn}(z)$. The following proposition shows that multiplying a covariate by the signum of its correlation with the residual yields a positive correlation with the residual.

Proposition 1

*Given covariate $\bm{x}_{j}$ and current residual $(\bm{y}-\bm{\hat{\mu}}_{S})$, the inner product $\left<s_{j}\bm{x}_{j},\bm{y}-\bm{\hat{\mu}}_{S}\right>$, where $s_{j}=\text{sgn}_{\mathbbm{C}}\left(\left<\bm{x}_{j},\bm{y}-\bm{\hat{\mu}}_{S}\right>\right)$ is a positive, real value.

Proof:*

[TABLE]

$\square$**

Proposition 1 shows that the signum of the complex correlation provides a method of aligning a covariate so that it is positively correlated with the residual. Using the signum-aligned covariates in the active set, the next step of LARS is the construction of an equiangular search direction. The next proposition provides the new direction of travel.

Proposition 2

Let an equi-inner product condition for covariates in the active set be given by

[TABLE]

The resulting inner product value $L_{S}$, active set covariate weights $\bm{w}_{S}$, and the equi-inner product direction $\bm{u}_{S}$ are, respectively,

[TABLE]

Proof:* Let $X_{S}\in\mathbbm{C}^{n\times q}$ be a nonempty active set where n is the covariates dimension and $q$ is the number of covariates. Since $X_{S}$ (by construction), we have that $\text{rank}(X_{S})=q$. Therefore, the associated $q\times q$ Gram matrix $X_{S}^{H}X_{S}$ has [20]*

[TABLE]

The matrix $X_{S}^{H}X_{S}$ is therefore invertible. Since we construct $\bm{u}_{S}$ to reside in the column space of $X_{S}$, then the projection of $\bm{u}_{S}$ onto the column space of $X_{S}$ using the projection operator $X_{S}(X_{S}^{H}X_{S})^{-1}X_{S}^{H}$ returns $\bm{u}_{S}$, providing the identity

[TABLE]

Taking the inner product of both sides with $\bm{u}_{S}$ and invoking the constraint $\bm{u}_{S}^{H}\bm{u}_{S}=1$, yields

[TABLE]

Inserting the equi-inner product condition (20) and solving leads to

[TABLE]

Multiplying the equi-inner product condition (20) on the left by $X_{S}(X_{S}^{H}X_{S})^{-1}$ and invoking the identity (24) provides the travel direction

[TABLE]

$\square$**

Note that the change in inner product results in minimal alteration of the calculations for the travel direction. The third proposition of this section gives the equation for length of travel along $\bm{u}_{S}$.

Proposition 3

If the active set $S$ and its complement $S^{c}$ are nonempty, then traveling along the equi-inner product direction $\bm{u}_{S}$ by the length

[TABLE]

*causes a new covariate to enter the active set.

Proof: After a step $\gamma\bm{u}_{S}$ in the estimate, the absolute current correlation between each covariate $\bm{x}_{j}$ and the residual becomes,*

[TABLE]

where $g_{j}=\left<\bm{x}_{j},\bm{u}_{S}\right>$. For a member of the active set $S$, the current absolute correlation becomes $\left\lvert\hat{C}-\gamma L_{S}\right\rvert$, since $\hat{c}_{j}=\hat{C}$, and $g_{j}=L_{S}$, for each $j\in S$. Equating the absolute current correlations provides the condition

[TABLE]

When condition (26) holds, a new covariate achieves an equal absolute correlation and enters the active set $S$.

Note that since $X^{H}X$ is positive-definite, $\left(X^{H}X\right)^{-1}$ is positive definite. Further, note that $L_{S}$ given in (21) contains a quadratic form of a positive-definite matrix, and $L_{S}$ itself turns out to be a strictly positive, real number. Since $\hat{C}>0$ and $L_{S}>0$, the condition (26) leads to the quadratic equation

[TABLE]

The solution of this equation is given in (25) with the additional restriction that we seek $\hat{\gamma}$ to be the minimum value with $\gamma>0$ such that condition (26) holds. $\square$

Propositions 1 through 3 provide the necessary tools for adapting LARS for complex covariates and complex data. If the covariates exist in complex conjugate pairs, as often occurs for real-valued data, note that LARS will add both covariates in a conjugate pair to the active set during the same iteration, since the current absolute correlations are equal when $y$ and $\bm{\hat{\mu}}_{S}$ are real, i.e., $\left\lvert\left<\bm{x}_{j},\bm{y}-\bm{\hat{\mu}}_{S}\right>\right\rvert=\left\lvert\left<\overline{\bm{x}}_{j},\bm{y}-\bm{\hat{\mu}}_{S}\right>\right\rvert$.

4 Least Angle Regression for Dynamic Mode Decomposition

DMD often provides complex DMD modes and DMD eigenvalues. Section III provides the necessary modification of the LARS algorithm for complex quantities, and this section applies Complex LARS to DMD mode selection. To apply Complex LARS to the problem of DMD mode selection, we seek a mapping between the DMD modes and the covariates in Complex LARS so that covariate selection equates to mode selection.

Observe that Complex LARS takes a zero-mean data vector $\bm{y}$ as input. Let $\text{vec}(\cdot)$ denote the operation of vectorizing a two-dimensional array into a column vector by stacking the columns of the array. Vectorization of the data matrix $\Psi_{0}$ given in (1) and subtraction of the mean provides a zero-mean data vector

[TABLE]

containing the complete temporal evolution of the measurements. Similarly, to capture the temporal evolution of the DMD modes in vector form, consider the Vandermonde matrix $\Xi$ in (12), which captures the temporal evolution of the eigenvalue coefficients of the DMD modes. Define the $j$th row of $\Xi$ to be $\bm{\xi}_{j}\in\mathbb{C}^{1\times N}$ such that

[TABLE]

Using the Vandermonde row vectors, each covariate $\bm{x}_{j}$ for $j=1,\dots,r$, can be constructed by

[TABLE]

where the variance $\text{var}(\bm{q})=1/n\sum_{j=1}^{n}|q_{j}-\text{mean}(\bm{q})|^{2}$ for an $n$-dimensional vector $\bm{q}$. Each resulting $\bm{x}_{j}$ is zero-mean, unit-variance, and lies in $\mathbb{C}^{mN\times 1}$, where $m$ is the length of each snapshot and $N$ is the number of snapshots. The covariates $\bm{x}_{j}$ for $j=1,\dots,r,$ consist of the time resolved DMD modes but are not scaled by the DMD mode amplitudes. Appropriate scaling of the DMD mode amplitudes is the reduced-order modeling problem. Algorithm 3 presents LARS4DMD, the reduced-order modeling algorithm proposed by this paper.

**Algorithm 3 (LARS4DMD)

**

Inputs: Data matrices $\Psi_{0}$ and $\Psi_{1}$.

1. 3.1)

Calculate the DMD modes $\bm{\phi}_{j}$ and the DMD eigenvalues $\lambda_{j}$ for $j=1,\dots,r$, from $\Psi_{0}$ and $\Psi_{1}$ using the SVD of $\Psi_{0}$, the eigendecomposition of (3), and equation (4). 2. 3.2)

Form the data vector $\bm{y}$

[TABLE] 3. 3.3)

Form zero-mean, unit-variance covariates

[TABLE]

for $j=1,\dots,r$. 4. 3.4)

Call Complex LARS, forming the LARS4DMD estimate

[TABLE]

within Step 2.6. Within Step 2.7, also calculate the LARS4DMD regression coefficients

[TABLE]

Output: Vector of regression coefficients $\alpha^{\text{LD}}_{k}$ for each iteration $k$.

Once the zero-mean, unit-variance covariates $\bm{x}_{j}$ and the zero-mean data $\bm{y}$ are formed in Steps 3.1 through 3.3, the user can perform Complex LARS to obtain regression coefficients. However, the standard regression coefficients from LARS assume that the algorithm is allowed to run to completion (cf. Steps 1.9 Algorithm 1 and Step 2.9 of Algorithm 2). Generating a model using fewer modes than the total number available can be thought of as pausing the Complex LARS algorithm and temporarily discarding the remaining modes that are not in the active set (i.e., emptying $S^{c}$). The LARS algorithm (Algorithm 1) provides a separate calculation for the travel distance along the equiangular direction when $S^{c}$ is empty, given by $\hat{\gamma}=\hat{C}/L_{S}$. The authors of [15] note that this choice results in regression coefficients that correspond with the least-squares solution assuming a model based on the covariates in the active set only. Step 3.4 in Algorithm 3 implements this choice for travel distance when calculating the LARS4DMD estimate $\bm{\mu}_{\text{LD},k}$ and the LARS4DMD regression coefficients $\bm{\alpha}^{\text{LD}}_{k}$ at each iteration $k$. We note that this choice of travel distance $\hat{\gamma}$ corresponds to variant of LARS known as the LARS Ordinary Least Squares (OLS) hybrid method in [15] that performs re-fitting of the regression coefficients using least squares on the selected covariates.

Each $j$th regression coefficient $\alpha^{\text{LD}}_{k,j}$ corresponds to the $j$th covariate $\bm{x}_{j}$, which is a manipulated version of the corresponding $j$th DMD mode $\bm{\phi}_{j}$ and Vandermonde row $\bm{\xi}_{j}$. To generate a reconstruction of original dataset, calculate

[TABLE]

where the addition of of the scalar value $\text{mean}(\text{vec}(\Psi_{0}))$ occurs elementwise to add back the data mean that was previously subtracted when generating $y$ in (27). Reshaping the $\text{vec}\left(\Psi^{\text{LD,rec}}_{0}\right)$ vector provides the reconstruction of $\Psi_{0}$. The next section performs a numerical comparison of the performance loss in the reconstuction of $\Psi_{0}$ by the LARS4DMD and DMDSP algorithms.

5 Numerical experiment: Poiseuille flow

This section compares the performance and mode selection of the DMDSP and LARS4DMD algorithms on a Poiseuille flow test case. Pressure-driven flow in a channel between two parallel walls (of infinite depth perpendicular to the flow) is known as plane Poiseuille flow [21]. The Poiseuille flow solution can be obtained from the two-dimensional, linearized Navier-Stokes equations. These equations are solved numerically at Reynolds Number $\text{Re}=10,000$ using a pseudo-spectral scheme in [9]. The first author of [9] generously provides the data from the solution on the author’s website at [16]. Please refer to [9] for details on the numerical methods used in generation of the Poiseuille flow dataset.

Using the numerically generated data for Poiseuille flow from [9, 16], we implemented the DMDSP and LARS4DMD algorithms to obtain models of various sizes. For DMDSP, we chose the regularization parameter $\beta$ to range logarithmically from $5\times 10^{-6}$ to $160$, providing $2,500$ values for $\beta$. Although DMDSP treats each $\beta$ value as a separate optimization problem, we examine the sequence of models from DMDSP that result for decreasing $\beta$ values to compare with sequential model construction procedure of LARS4DMD.

Figure 1 compares the performance loss of the dataset reconstruction for these two methods. Note that LARS4DMD and DMDSP provide the same performance loss for a one-mode model. As $\beta$ decreases, DMDSP often produces the same model several times in a row. After crossing a certain $\beta$ threshold, DMDSP changes model size by adding one or more modes to the model. The performance loss for DMDSP and LARS4DMD remain close for increasing model size. For a two-mode model, LARS4DMD provides better performance. However, for model sizes 3 through 8, DMDSP provides better performance. LARS4DMD lies below the DMDSP performance loss for model sizes 9 through 20. For larger model sizes, the two methods yield comparable values. Across all model sizes, DMDSP and LARS4DMD yield similar performance loss, with DMDSP only notably outperforming LARS4DMD in the intermediate range of model sizes 3 through 8.

Figure 1 contains a black dot for each $\beta$ value in DMDSP, and each marker denotes a corresponding system model generated by DMDSP. Note that there are ranges of $\beta$ values that yield models of the same size. In total, there were only $26$ unique models (i.e., models with unique $\alpha$ vectors) out of the $2,500$ models generated by DMDSP, and there were only $24$ unique model sizes. DMDSP skips model size 17. DMDSP generated two models containing $21$ DMD modes and two models containing $22$ DMD modes. In contrast, LARS4DMD sequentially constructed $26$ unique models by stopping the LARS procedure at each iteration to generate a new model.

Figures 2(a) and 2(b) illustrate the DMD modes selected in the unique models constructed by LARS4DMD and DMDSP, respectively. The shading of each square represents the absolute value of the DMD mode amplitude $|\alpha_{j}|$ for the $j$th mode. For some model sizes, LARS4DMD and DMDSP produce very similar models. For example, the one-mode model and the three largest models match very closely between the two techniques. Mode 1 is selected first and mode 2 is selected last by both algorithms. Generally, both methods weight the most recently added mode the least, however exceptions occur (e.g. in Fig. 2(a), the addition of mode 21 at model size 16 outweighs mode 17 which was added at model size 10). However, there are also notable differences for the models produced by these methods. For model size 2, LARS4DMD selects mode 7, but DMDSP selects mode 26. However, DMDSP eventually scales mode 26 to zero in model 21(b). Note that DMDSP has the ability to remove modes from the model, because DMDSP solves a separate optimization problem for each $\beta$ value.

DMDSP produces two models each for model sizes 21 and 22. In order of increasing regularization parameter $\beta$, DMDSP increases from model size 21 to 22, then decreases back to 21 before increasing again to 22. This sequence can be seen in inset plot of Figure 1. This behavior is also depicted in Figure 2(b) where regularization parameter $\beta$ decreases from left to right.

6 Conclusion

In this paper, we adopt a useful algorithm known as Least Angle Regression (LARS) from the statistics and machine learning literature and adapt it for construction of reduced-order DMD models. Utilizing a complex inner product, we create a complex version of LARS that can be applied to select DMD modes, which are often complex. We refer to the DMD-based reduced-order modeling algorithm as LARS4DMD. LARS4DMD constructs a reduced-order model in a sequential manner by selecting the DMD modes that are highly correlated with the residual between the measured data and the current estimate. The results from LARS4DMD are comparable to results from Sparsity-Promoting Dynamic Mode Decomposition (DMDSP), a popular DMD-based modeling algorithm. In ongoing work, we are testing the performance of LARS4DMD on data collected experimentally for which measurement noise is present.