Krylov Dynamic Mode Decomposition (DMD)

• Krylov DMD is an optimized variant of dynamic mode decomposition that integrates time-delay embedding with Krylov subspace projections to capture oscillatory modes in high-dimensional data.
• It reduces computational complexity and memory usage by projecting large datasets onto a lower-dimensional subspace before applying SVD and spectral decomposition.
• Empirical benchmarks demonstrate that Krylov DMD achieves nearly identical modal accuracy as classic DMD while significantly lowering FLOPs and runtime for large-scale systems.

Krylov Dynamic Mode Decomposition (DMD) is an optimized variant of the Dynamic Mode Decomposition framework designed for efficient analysis of high-dimensional, highly oscillatory spatiotemporal datasets. By incorporating both time-delay coordinates (TDC) and Krylov-subspace-inspired projections, this methodology achieves significant computational and memory reductions with negligible loss in modal accuracy. Krylov DMD addresses bottlenecks present in standard DMD by projecting high-dimensional data onto a carefully constructed low-dimensional subspace before matrix factorization and spectral decomposition, enabling robust identification of oscillatory modes in large-scale systems (Murshed et al., 2020).

1. Standard DMD and Computational Challenges

Standard DMD operates on snapshot matrices $X \in \mathbb{R}^{M \times N}$, where $M$ is the spatial dimension and $N$ the number of temporal samples. The snapshot pair is split as $X_1 = [x_1, ..., x_{N-1}]$ and $X_2 = [x_2, ..., x_N]$, seeking the best-fit linear operator $A$ such that $X_2 \approx A X_1$. The core algorithm involves the following steps:

1. Economic-size singular value decomposition (SVD): $X_1 = U \Sigma V^*$, truncated to rank $r$.
2. Low-dimensional operator construction: $\tilde{S} = U^* X_2 V \Sigma^{-1}$.
3. Spectral decomposition: $\tilde{S} y_k = \mu_k y_k$, DMD modes $\phi_k = U y_k$, eigenvalues $\omega_k = \ln(\mu_k)/\Delta t$.
4. State reconstruction: $$x_{\text{DMD}}(t) = \Phi\,\mathrm{diag}(e^{\omega t})\,b$$, with $b = \Phi^\dagger x_1$.

However, when $M$ and/or $N$ are large, the SVD step induces computational complexity $O(\min(M,N)^2 \max(M,N))$ and memory usage $O(MN)$. Data from many relevant phenomena are “big” ($M \gg 10^3$) and “highly oscillatory,” requiring dense sampling ($N$ large), which amplifies these bottlenecks. Standard DMD may also fail to resolve oscillatory modes in non-Markovian or coarsely-sampled systems (Murshed et al., 2020).

2. Time-Delay Coordinates and the Need for Projection

Time-Delay Coordinates augment each snapshot $x_k$ by stacking $p$ past states, forming highly informative, tall Hankel-like matrices:

$X_{1,\text{aug}} = \begin{bmatrix} x_1 & \hdots & x_{q-1} \ x_2 & \hdots & x_q \ \vdots \ x_p & \hdots & x_{N-1} \end{bmatrix} \in \mathbb{R}^{(p+1)M \times (N-p)}$

This procedure effectively embeds hidden oscillatory dynamics into a higher-dimensional, approximately linear, manifold, enhancing the extractability of underlying coherent structures. However, it increases the row dimension to $qM$, drastically elevating the cost of subsequent SVD. Mitigation involves the application of projection operators $R \in \mathbb{R}^{a \times (qM)}$, with $a \ll qM$, resulting in compressed matrices $Z_{1,\text{aug}} = R X_{1,\text{aug}}$, $Z_{2,\text{aug}} = R X_{2,\text{aug}}$. The projection must preserve dominant DMD eigenvalues and modes while reducing both floating-point operations and memory requirements (Murshed et al., 2020).

3. Krylov-Subspace Projections

The Krylov subspace, $$\mathcal{K}_m(A, b) = \mathrm{span}\{b, Ab, A^2 b, ..., A^{m-1} b\}$$ for $A \in \mathbb{R}^{n \times n}$ and seed vector $b \in \mathbb{R}^n$, is a classical construct from numerical linear algebra exploited here to define low-dimensional projectors. Its construction uses the Arnoldi process (a modified Gram–Schmidt):

• Inputs: random matrix $A$, seed vector $b$.
• Iteratively generate orthonormal basis $V_m = [v_1, ..., v_m]$ for the subspace.
• The projection operator is then $R_K = V^* \in \mathbb{R}^{(a+1) \times M}$.

For data $X \in \mathbb{R}^{M \times N}$, $Z = R_K X = V^* X \in \mathbb{R}^{(a+1) \times N}$ yields a compressed representation. The orthogonal projector is $P = VV^*$, but $R_K$ suffices for right-multiplicative projection. This basis, once constructed, is data-agnostic and can be reused if $A$ and $b$ are fixed.

4. Krylov-Projected TDC-DMD: Algorithmic Workflow

The Krylov-projected time-delay DMD framework proceeds as follows:

1. TDC Embedding: Form augmented snapshot matrices $X_{1,\text{aug}}$, $X_{2,\text{aug}}$ of size $qM \times (N-q)$.
2. Projector Construction: Use Arnoldi on $A$, $b$ to obtain $V \in \mathbb{R}^{M \times (a+1)}$; set $R_K = V^*$.
3. Data Projection: $Z_1 = R_K X_{1,\text{aug}}$, $Z_2 = R_K X_{2,\text{aug}}$.
4. SVD: $Z_1 = U_Z \Sigma_Z V_Z^*$, truncated to $r \ll a+1$.
5. Koopman Matrix: $\tilde{S} = U_Z^* Z_2 V_Z \Sigma_Z^{-1}$.
6. Spectral Decomposition: $\tilde{S} y_k = \mu_k y_k$, $\phi_k = U_Z y_k$.
7. Amplitudes and Frequencies: $b = \Phi^\dagger x_1$, $\omega_k = \ln(\mu_k)/\Delta t$.
8. Reconstruction: $$x_{\rm DMD}(t) = \Phi \,\mathrm{diag}(e^{\omega t}) b$$.

This workflow compresses the effective row-dimension from $qM$ to $a+1$, where $a+1$ is the Krylov subspace dimension, delivering a dramatic reduction in resource consumption (Murshed et al., 2020).

5. Computational Complexity and Memory Analysis

A direct comparison of various DMD strategies for data of dimension $M$ and projection size $a$ is as follows:

Method	SVD Cost	Projection/Arnoldi Cost	Memory Usage
Standard TDC‐DMD	$O(\min(qM, N)^2 \max(qM, N))$	—	$O(qM \cdot N)$
Sampling/Gaussian Proj	$O(a^2 N)$ (after proj)	$O(a q M N)$	$O(a N)$
Krylov‐DMD	$O(a^2 N)$ (after proj)	$O(a M^2)$ (Arnoldi, once); $O(a M q N)$ (proj)	$O(M a)$ (basis) + $O(a N)$ (proj)

For $M = 10^4$–$10^6$ and $a = 10^2$–$10^3$, Krylov DMD achieves a reduction of one to two orders of magnitude in both FLOPs and memory. If the Krylov basis $V$ is constructed once and reused, the amortized cost is much lower than standard DMD. This enables practical application to very large datasets that would otherwise be infeasible with direct SVD (Murshed et al., 2020).

6. Experimental Benchmarks and Modal Accuracy

Empirical evaluation is performed on two canonical data sets:

• Double Gyre vorticity ($[0,2] \times [0,1]$, $M=10^4$, $N=200$).
• 2D compressible signal, two-frequency ($M=10^4$, $N=200$).

Time-delay embedding uses $q=2$. Projection dimensions for Krylov DMD and comparators are selected as $a=100$ (Double Gyre) and $a=50$ (Signal). The singular value truncation rank is $r=20$.

Results summary:

• Eigenvalue spectra (Im$(\lambda)$ vs Re$(\lambda)$) are nearly identical across projection variants, forming symmetric clusters near the unit circle, indicative of stable dynamics.
• Long-run reconstruction errors $\|x_{\text{true}}(t) - x_{\text{DMD}}(t)\|_2$ for Krylov DMD and sparse projection DMD are nearly indistinguishable from classic DMD when $a \cdot q \geq r$ (e.g., $100 \cdot 2 \geq 20$).
• Sampling-based methods can omit weak modes, resulting in occasional long-term drift.
• Runtime benchmarks (per trial): classic TDC-DMD $\sim3.2$s, sampling DMD $\sim0.12$s, Gaussian projection $\sim0.18$s, sparse projection $\sim0.15$s, Krylov DMD (amortized) $\sim0.20$s (Arnoldi step $\sim0.05$s).

7. Theoretical Considerations, Limitations, and Extensions

Krylov DMD offers substantive advantages:

• High data reduction: compresses $M q$ to $a$, with $a \ll M q$, without significant loss in modal accuracy.
• Robust mode identification in highly oscillatory regimes through joint TDC and Krylov projection.
• Compatibility with sparse storage where appropriate.

Potential limitations include:

• Initial Arnoldi cost $O(a M^2)$, which can be burdensome for very large $M$ unless $V$ is reused.
• The accuracy and stability depend on the choice of $a$, $A$, and $b$ (projection dimension and seeding), necessitating problem-dependent cross-validation.
• Error bounds for combined TDC/Krylov projection are not fully characterized.

Possible directions for further development include streaming/online Arnoldi for time-varying data, adaptive selection of $a$ for energy target capture, application to $\ell_1$-based compressive DMD, and extension to PDE-constrained inverse problems or control-oriented settings (Murshed et al., 2020).