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.

[dylan-copeland]

There are significant changes, especially in the cpp file, in the existing class IncrementalSVDBrand, to accommodate incremental DMD. Can this PR be refactored to avoid changes to IncrementalSVDBrand? It is understandable that the data from IncrementalSVDBrand needs to be accessed for incremental DMD, so it makes sense to add access functions. Hopefully it is possible to avoid changing so much of the class's implementation, by using access functions, making a new class derived from IncrementalSVDBrand, or some other way.

[swsuh28]

I've somehow fixed this problem in 06b4d6b. The lengthy IncrementalDMDInternal things have been removed, but the thing is that we still need to intercept some of the intermediate matrices from IncrementalSVDBrand for IncrementalDMD. Specifically, matrices before IncrementalSVDBrand::addNewSample or IncrementalSVDBrand::addLinearlyDependentSample are accessed by IncrementalDMD. For that, I have added extra pointers pointing to the matrices before being updated. (See Line 176-190 of lib/linalg/svd/IncrementalSVDBrand.h. IncrementalDMD is also declared as a friend class of IncrementalSVDBrand so it can access those private pointers. I find that this makes minimal changes to IncrementalSVDBrand, but if you have better ideas, please let me know.

[dylan-copeland]

Thanks for the improvement. There is still a lot of new code in IncrementalSVDBrand.cpp, and it is hard to tell whether it affects the existing incremental SVD or just adds additional intermediate matrices needed by IncrementalDMD. If it should not change the SVD algorithm, and if it only stores intermediate matrices, can this be done optionally by setting a flag? My concerns are (i) verifying that the old incremental SVD is not changed, and (ii) not adding extra computation or memory usage when IncrementalDMD is not being used.

[dylan-copeland]

One way to separate the new code from the existing IncrementalSVDBrand might be to make a class derived from IncrementalSVDBrand, containing the additional data members. Functions like buildIncrementalSVD, updateTemporalBasis, and addNewSample could become virtual. Some code from the class IncrementalSVDBrand would be reused, while the new code is added only in the derived class.