Stochastic Hankel DMD Overview

• Stochastic Hankel DMD is a technique for approximating the stochastic Koopman operator from delay-embedded data, enabling robust spectral analysis in noisy systems.
• It constructs Hankel matrices from time-delay observables and applies low-rank SVD to extract eigenvalues and eigenfunctions from random dynamical systems.
• Bayesian and noise-robust variants improve model validation, forecasting, and reduced order modeling in applications like wind turbines and nonlinear oscillators.

Stochastic Hankel Dynamic Mode Decomposition (Stochastic Hankel DMD, sHankel-DMD, or SHDMD) is a family of operator-theoretic numerics for extracting spectral information—particularly Koopman eigenvalues and eigenfunctions—from stochastic or random dynamical systems (RDS) via time-delay–embedded data matrices. By generalizing classical (deterministic) Hankel DMD to accommodate either intrinsic dynamical noise or epistemic uncertainty in hyperparameters, these algorithms rigorously approximate the stochastic Koopman operator, enabling robust spectral analysis, forecasting, and model reduction in noisy, real-world systems (Črnjarić-Žic et al., 2017, Wanner et al., 2020, Palma et al., 2024).

1. Mathematical Foundations for Stochastic Koopman Operator Approximation

Let $M$ denote a measurable state space (typically a subset of $\mathbb{R}^d$), with probability space $(\Omega,\mathcal F,P)$ governing random variables $\omega$. A random dynamical system (RDS) is described by a cocycle map $\varphi(t,\omega,x)$, with Markov/semigroup properties in the time parameter. The stochastic Koopman operator family $(\mathcal{K}^t)_{t\ge0}$ acts on scalar observables $f:M\rightarrow\mathbb C$ via

$$\mathcal{K}^t f(x) = \mathbb{E}[f(\varphi(t,\omega,x))]$$

where the expectation is with respect to randomness in the dynamics.

For time-homogeneous (Markovian) RDS, $\{\mathcal{K}^k\}$ forms a semigroup: $\mathcal{K}^{s+t}=\mathcal{K}^s\circ\mathcal{K}^t$. The spectral decomposition (eigenvalues, eigenfunctions) of $\mathbb{R}^d$0 encodes the linear evolution of nonlinear random flows in infinite dimensions and serves as the target for sHankel-DMD algorithms (Črnjarić-Žic et al., 2017, Wanner et al., 2020). Under ergodic assumptions, time averages along trajectories can be used to approximate ensemble spectral quantities.

2. Stochastic Hankel Matrix Construction

Extending classical delay-embedding ideas, the stochastic Hankel matrix is built from expectations of delayed observables under stochastic flow:

• For a fixed observable $\mathbb{R}^d$1, initial $\mathbb{R}^d$2, and integer $\mathbb{R}^d$3, the stochastic trajectory vector

$\mathbb{R}^d$4

is defined, with $\mathbb{R}^d$5 as the time-one map.

• Its $\mathbb{R}^d$6-th shifted, expectation-valued version is

$\mathbb{R}^d$7

• Assemble the $\mathbb{R}^d$8 stochastic Hankel matrix

$\mathbb{R}^d$9

The columns correspond to repeated application of $(\Omega,\mathcal F,P)$0 to $(\Omega,\mathcal F,P)$1, evaluated along the stochastic trajectory. This construction underlies sHankel-DMD for Koopman spectral approximation (Črnjarić-Žic et al., 2017).

In alternative formalism, delay-augmented snapshots of a multi-variate process $(\Omega,\mathcal F,P)$2 form the Hankel matrices

$(\Omega,\mathcal F,P)$3

for the embedding order $(\Omega,\mathcal F,P)$4, supporting application to both stochastic and Bayesianized settings (Palma et al., 2024).

3. Stochastic Hankel DMD Algorithms and Convergence

3.1 Residual-based sHankel-DMD

The canonical algorithm proceeds as follows (Črnjarić-Žic et al., 2017):

1. Data Preparation: Estimate $(\Omega,\mathcal F,P)$5 for $(\Omega,\mathcal F,P)$6 via ensemble or time/ergodic averaging from observed or simulated RDS data.
2. Snapshot Matrices: Form $(\Omega,\mathcal F,P)$7 and $(\Omega,\mathcal F,P)$8.
3. Low-rank SVD: $(\Omega,\mathcal F,P)$9 with rank $\omega$0 selected such that $\omega$1 for a fixed tolerance $\omega$2.
4. Rayleigh Quotient: $\omega$3.
5. Eigenproblem: Solve $\omega$4; form Ritz eigenvectors $\omega$5.
6. Residual Filtering: Retain modes with $\omega$6.
7. Eigenfunction Reconstruction: Project DMD modes as approximations to Koopman eigenfunctions.

Convergence Guarantee

Given ergodicity, Markovianity (semigroup property), and the existence of a finite-dimensional $\omega$7-invariant Krylov subspace spanned by $\omega$8, the eigenvalues/eigenfunctions computed using sHankel-DMD converge almost surely to those of $\omega$9 as trajectory length $\varphi(t,\omega,x)$0. The proof leverages Birkhoff’s ergodic theorem for empirical inner products, reduction to a companion matrix representation on invariant subspaces, and spectral continuity (Črnjarić-Žic et al., 2017). Under uniform Gramian conditioning and mixing, the convergence rate satisfies $\varphi(t,\omega,x)$1 with high probability (Wanner et al., 2020).

3.2 Noise-robust and Bayesian/Stochastic Variants

Random dynamical systems and noisy observables introduce bias in standard DMD/Hankel-DMD by contaminating covariance estimates. The noise-robust approach introduces statistically independent “dual” observables $\varphi(t,\omega,x)$2 to build cross-covariance matrices $\varphi(t,\omega,x)$3, $\varphi(t,\omega,x)$4, $\varphi(t,\omega,x)$5, and forms debiased estimators: $\varphi(t,\omega,x)$6 This approach provably removes the leading noise bias in the limit of large data and under suitable independence (Wanner et al., 2020).

In applications where DMD hyperparameters (e.g., time-delay order, trajectory window) are uncertain or system parameters vary, a Bayesian formulation treats these as random variables. Under this “SHDMD” paradigm:

• Uniform priors are placed on the training window length ($\varphi(t,\omega,x)$7) and delay horizon ($\varphi(t,\omega,x)$8).
• Monte Carlo sampling over $\varphi(t,\omega,x)$9 yields an ensemble of DMD operators, modes, and spectra.
• Posterior predictive means and variances provide point forecasts and credible confidence bands for future states or observables (Palma et al., 2024).

4. Error Metrics, Model Validation, and Numerical Examples

Quantitative error metrics, central for validation and model selection, include:

• Normalized Root Mean Square Error (NRMSE):

$(\mathcal{K}^t)_{t\ge0}$0

• Normalized Average Min/Max Absolute Error (NAMMAE):

$(\mathcal{K}^t)_{t\ge0}$1

• Jensen-Shannon Divergence (JSD): Empirical PDFs of test and forecast trajectories are binned, and

$(\mathcal{K}^t)_{t\ge0}$2

Posterior uncertainty in error metrics is accessible in the Bayesian/MC setting (Palma et al., 2024).

Representative test cases include:

• Noisy circle rotations: sHankel-DMD recovers Koopman eigenvalues $(\mathcal{K}^t)_{t\ge0}$3 under uniform-noise perturbations; standard Hankel DMD yields biased results (Črnjarić-Žic et al., 2017, Wanner et al., 2020).
• Noisy Van der Pol oscillator: sHankel-DMD extracts the full lattice of Koopman eigenvalues and corresponding spectral components, matching theoretical predictions (Črnjarić-Žic et al., 2017).
• Floating offshore wind turbine dynamics: Bayesian SHDMD delivers both forecasts and posterior uncertainty bands for real-world multi-channel time series (Palma et al., 2024).
• Linear maps with additive/multiplicative noise and stochastic Stuart-Landau oscillators: sHankel-DMD achieves unbiased eigenvalue recovery and accurately captures spectral decay/oscillations (Wanner et al., 2020).

5. Implementation Notes and Computational Aspects

Efficient sHankel-DMD relies on well-conditioned linear algebra on large Hankel matrices. Key aspects include:

• Singular Value Decomposition (SVD): SVD of $(\mathcal{K}^t)_{t\ge0}$4 matrices at $(\mathcal{K}^t)_{t\ge0}$5 complexity; truncated versions promote stability and allow low-rank model reduction.
• Residual-based trimming (DMD RRR): Filtering out high-residual modes ameliorates the influence of sampling noise.
• Hyperparameter choices: Embedding dimension $(\mathcal{K}^t)_{t\ge0}$6 must exceed the number of modes; window $(\mathcal{K}^t)_{t\ge0}$7, rank cutoff $(\mathcal{K}^t)_{t\ge0}$8 trades statistical power versus overfitting.
• Normalization and filtering: Z-score normalization and low-pass signal filtering improve numerical conditioning and noise suppression (Palma et al., 2024).
• Monte Carlo/Bayesian loop: In stochastic settings, ensemble runs over randomly drawn hyperparameters allow credible intervals and uncertainty quantification with negligible additional computational overhead (Palma et al., 2024).

The methodology is compatible with fast prototyping in MATLAB (economy SVD), NumPy/SciPy, or randomized SVD variants for larger-scale or real-time deployment (Palma et al., 2024).

6. Generalization and Relations to Other Approaches

Stochastic Hankel DMD extends deterministic Hankel DMD, enabling operator-theoretic spectral analysis—and thus equation-free, reduced modeling—even under stochastic forcing or observational noise. Unlike standard DMD, which suffers from bias and instability in the presence of dynamical or measurement randomness, sHankel-DMD and its robust or Bayesian variants directly target the expected evolution in the sense of the stochastic Koopman operator, providing strong theoretical convergence and empirical reliability (Črnjarić-Žic et al., 2017, Wanner et al., 2020).

A plausible implication is that sHankel-DMD offers advantages for systems identification, forecasting under uncertainty, and digital-twin model reduction, where quantification of predictive confidence is essential (Palma et al., 2024).

7. Applications and Practical Impact

sHankel-DMD and SHDMD have been validated across a broad spectrum of systems:

• Stochastic differential equations (Van der Pol, Stuart–Landau)
• Randomly perturbed maps (circle rotations, linear systems with additive and multiplicative noise)
• Experimental, industrial datasets (floating offshore wind turbines with high-dimensional, multi-modal, physical sensor data)

Empirical results consistently demonstrate that sHankel-DMD recovers theoretically correct Koopman spectra (in both decay rates and frequencies), effectively filters noise, quantifies uncertainty, and supports robust system identification and short-term forecasting. The Bayesian variant enables rigorous propagation of hyperparameter uncertainty into prediction intervals, supporting continuous, real-time digital twin applications where online retraining and resampling are feasible (Palma et al., 2024).

The method thus generalizes classical spectral decomposition for deterministic systems to the stochastic regime, preserving computational tractability, extensibility to control or exogenous inputs, and compatibility with large-scale, real-time scenarios.