Incremental DMD

[swsuh28]

Incremental DMD

Files added:

• lib/algo/IncrementalDMD.h and lib/algo/IncrementalDMD.cpp (Incremental DMD class)
• examples/dmd/heat_conduction_otf.cpp (Demo for heat equation solver w/ online DMD)

Files modified:

• lib/linalg/svd/IncrementalSVDBrand.h and lib/linalg/svd/IncrementalSVDBrand.cpp (save matrices to pass to IncrementalDMD)
• lib/linalg/BasisGenerator.h (add IncrementalDMD as friend)
• CMakeList.txt and lib/CMakeList.txt (include demo and add to libraries)

Background

IncrementalDMD constructs DMD in a fast, online fashion, leveraging the incremental nature of IncrementalSVD introduced by [Brand, 2006] (https://doi.org/10.1016/j.laa.2005.07.021). Starting from the original DMD formulation, let us first define snapshot matrices $$F^+ = [f^1 \vert f^2 \vert \cdots \vert f^m] \in \mathbb R^{N \times m} \quad \text{and} \quad F^- = [f^0 \vert f^1 \vert \cdots \vert f^{m-1}] \in \mathbb R^{N \times m}$$ where $f^i \in \mathbb R^N$ is a snapshot that has $N$ degrees of freedom at time $i$ for $i=0,1,\cdots,m$. Then, we can construct the reduced DMD matrix as $$\tilde A_r = U^* F^+ W S^{-1} \in \mathbb R^{r \times r},$$ where $F^- \approx USW^*$ is the truncated SVD of the snapshot matrix $F^-$ at rank $r$. In other words, $$U \in \mathbb R^{N \times r}, \quad S \in \mathbb R^{r \times r},\quad \text{and} \quad W \in \mathbb R^{m \times r}.$$

In case new snapshots stream into the DMD and $F^\pm$ increments every iteration, it is efficient to update the SVD incrementally rather than constructing all the matrices from scratch repetitively. For a fast rank-1 update of the SVD matrices, IncrementalSVD first decomposes the matrices as $$USW^* = \bar U \bar U' S \bar W'^* \bar W^*$$ where $\bar U \in \mathbb R^{N \times r}$, $\bar U' \in \mathbb R^{r \times r}$, $\bar W \in \mathbb R^{m \times r}$, and $\bar W' \in \mathbb R^{r \times r}$. The trick is only to append the large matrices $\bar U$ and $\bar W$, and all the matrix-matrix multiplications are done on small matrices $\bar U'$ and $\bar W'$. The updating rules are as follows.

Given the error $p = f_\mathrm{new} - UU^* f_\mathrm{new}$, first do the SVD

Q= \begin{bmatrix} 
S & U^* f_\mathrm{new} \\
0 & \Vert p \Vert
\end{bmatrix}=U_Q S_Q W_Q^*.

Then,

\bar U_\mathrm{new} \leftarrow \begin{bmatrix} \bar U_\mathrm{old} \\ p/\Vert p \Vert \end{bmatrix}, \quad \bar U'_\mathrm{new} \leftarrow \begin{bmatrix} \bar U'_\mathrm{old} & 0 \\ 0 & 1 \end{bmatrix} U_Q, \quad S_\mathrm{new} \leftarrow S_Q,

\bar W_\mathrm{new} \leftarrow \begin{bmatrix} \bar W_\mathrm{old} & 0 \\ 0 & 1 \end{bmatrix}, \quad \bar W'_\mathrm{new} \leftarrow \begin{bmatrix} \bar W'_\mathrm{old} & 0 \\ 0 & 1 \end{bmatrix} W_Q, \quad \bar W'_\mathrm{new} \leftarrow W_Q^* \begin{bmatrix} \bar W'_\mathrm{old} & 0 \\ 0 & 1 \end{bmatrix}

if $p$ is linearly independent from $U$ and

\bar U_\mathrm{new} \leftarrow \bar U_\mathrm{old}, \quad \bar U'_\mathrm{new} \leftarrow \bar U'_\mathrm{old} U_Q[:-1,:-1], \quad S_\mathrm{new} \leftarrow S_Q[:-1,:-1],

\bar W_\mathrm{new} \leftarrow \begin{bmatrix} \bar W_\mathrm{old} \\ W_Q[-1,:] \bar W'^\dagger_\mathrm{new} \end{bmatrix}, \quad \bar W'_\mathrm{new} \leftarrow \bar W'_\mathrm{old} W_Q[:-1,:], \quad \bar W'^\dagger_\mathrm{new} \leftarrow W_Q[:-1,:] \bar W'^\dagger_\mathrm{old}

otherwise. The pseudo-inverse of $\bar W'$ can also be updated efficiently using the Shermann-Woodbury-Morrison formula. See [Brand, 2006] for more details.

Plugging these into our first equation for $\tilde A_r$,

\begin{align*}
\tilde A_{r,\mathrm{new}} &= \bar U'^*_\mathrm{new} \bar U^*_\mathrm{new} F^+_\mathrm{new} \bar W_\mathrm{new} \bar W'_\mathrm{new} S^{-1}_\mathrm{new} \\
&=U_Q^* \begin{bmatrix} \bar U'^*_\mathrm{old} & 0 \\ 0 & 1\end{bmatrix} \begin{bmatrix} \bar U^*_\mathrm{old} \\ p^*/\Vert p \Vert\end{bmatrix} [F^+ \ f_\mathrm{new}] \begin{bmatrix} \bar W_\mathrm{old} & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \bar W'_\mathrm{old} & 0 \\ 0 & 1 \end{bmatrix} W_Q S_Q^{-1} \\
&= U_Q^* \begin{bmatrix} \tilde A_{r,\mathrm{old}} S_\mathrm{old} & \bar U'^*_\mathrm{old} \bar U^*_\mathrm{old} f_\mathrm{new} \\ p^* F^+ \bar W_\mathrm{old} \bar W'_\mathrm{old} / \Vert p \Vert & p^* f_\mathrm{new} / \Vert p \Vert \end{bmatrix} W_Q S_Q^{-1}
\end{align*}

if $p$ is linearly independent, and

\begin{align*}
\tilde A_{r,\mathrm{new}} &= \bar U'^*_\mathrm{new} \bar U^*_\mathrm{new} F^+_\mathrm{new} \bar W_\mathrm{new} \bar W'_\mathrm{new} S^{-1}_\mathrm{new} \\
&=U_Q^*[:-1,:-1] \bar U'^*_\mathrm{old} \bar U^*_\mathrm{old} F^+_\mathrm{new} \begin{bmatrix} \bar W_\mathrm{old} \\ W_Q[-1,:] \bar W'^\dagger_\mathrm{new} \end{bmatrix} \bar W'_\mathrm{old} W_Q[:-1,:] S_Q[:-1,:-1]^{-1} \\
&= U_Q^*[:-1,:-1] \begin{bmatrix} \tilde A_{r,\mathrm{old}} S_\mathrm{old} \\ \bar U'^*_\mathrm{old} \bar U^*_\mathrm{old} f_\mathrm{new} \end{bmatrix} W_Q[:-1,:] S_Q[:-1,:-1]^{-1}
\end{align*}

otherwise.

In both cases, the matrix at the center that contains $\tilde A_{r,\mathrm{old}}$ is constructed most efficiently. $\tilde A_{r,\mathrm{old}} S_\mathrm{old}$ is $O(r^2)$ as $S_\mathrm{old}$ is diagonal. The two matrix-vector multiplications, $U'^*_{\mathrm{old}} U^*_{\mathrm{old}} f_{\mathrm{new}}$ and $p^* F^+ W_\mathrm{old} W'_\mathrm{old}$, are $O(Nr+r^2) \approx O(Nr)$ and $O(Nm+rm+r^2) \approx O(Nm)$ given $N \gg m > r$. The inner product $p^* f_\mathrm{new}$ is $O(N)$. Once the central matrix is constructed, all the remaining matrix-matrix multiplications cost at most $O(r^3)$.

IncrementalDMD::updateDMD is a straightforward implementation of the prescribed algorithm. Some remarks:

• $F^\pm$ is never constructed. Snapshots are stored in std::vector<Vector*> d_snapshots instead.
• $\bar W$ is pre-allocated since otherwise, its size also increments and would degrade the efficiency. This requires information on the maximum number of snapshots and the maximum number of bases to be used.
• All the matrices from IncrementalSVDBrand, $\bar U, \bar U', S, \bar W, \bar W'$, are needed in IncrementalDMD. As such, a private struct object IncrementalDMDInternal is defined in IncrementalSVDBrand.h. The matrices can be accessed by IncrementalSVDBrand::getAllMatrices.

Demonstration

Incremental DMD can be used to construct a reduced-order model (ROM) along a full-order model (FOM) simulation. This has been demonstrated on a heat conduction simulation, which is already included in examples/dmd/heat_conduction.cpp. Whenever a new snapshot is added, it is fed into the DMD. Then DMD predicts the solution at the next time step and if that prediction is accurate enough, the DMD prediction is appended as the new snapshot. Whether the DMD prediction is accurate is checked by measuring the norm of the residual when the predicted solution is plugged into the governing equation of the FOM. This is problem-specific, and can be found in ConductionOperator::Residual in examples/dmd/heat_conduction_otf.cpp.

[chldkdtn]

It would be nice to include regression test for this example.