High Order Dynamic Mode Decomposition (HODMD)

• HODMD is a data‐driven modal decomposition technique that extends classical DMD by incorporating time-delay embeddings to capture coherent spatio-temporal dynamics.
• It reformulates the evolution as a higher‐order autoregressive or companion matrix problem, increasing the effective state dimension for improved noise robustness.
• HODMD has been applied across diverse fields—from digital twins in agriculture to cardiovascular flow modeling—demonstrating enhanced forecasting accuracy and computational speed.

Searching arXiv for recent and relevant HODMD papers to ground the article. High Order Dynamic Mode Decomposition (HODMD) is a data-driven modal decomposition and reduced-order modelling technique that extends classical Dynamic Mode Decomposition (DMD) by incorporating time-delay embedding, or equivalently a higher-order linear recurrence, to extract coherent spatio-temporal structures from noisy, transient, and spectrally rich measurements. Its distinctive regime is the one in which the temporal or spectral complexity of the signal exceeds what the available spatial measurements can resolve, so that delayed coordinates are used to increase the effective dimension of the observable space without changing the sensor network. Across recent work, HODMD has been used for reduced-order forecasting, modal identification, feature extraction, clustering, and surrogate modelling in underground-farm digital twins, large-scale circuits, echocardiography, intraventricular hemodynamics, power systems, reacting flows, bubble columns, mechanical vibrations, and time-parallel fluid simulation (Conti et al., 2023, Tuor et al., 2023, Liu et al., 5 Aug 2025, Li et al., 10 Feb 2025, Groun et al., 2024).

1. Classical DMD and the higher-order extension

Classical DMD starts from uniformly sampled snapshots $x_k \in \mathbb{R}^n$ and forms the data matrices

$$X = [x_1, x_2, \dots, x_{m-1}], \qquad Y = [x_2, x_3, \dots, x_m].$$

After a thin or truncated SVD,

$X = U \Sigma V^*,$

one defines the reduced operator

$\tilde{A} = U^* Y V \Sigma^{-1},$

solves

$\tilde{A} W = W \Lambda,$

and recovers DMD modes

$\Phi = Y V \Sigma^{-1} W.$

With discrete-time eigenvalues $\lambda_i$, the continuous-time rates are

$$\omega_i = \frac{\log(\lambda_i)}{\Delta t} = \alpha_i + i 2\pi f_i,$$

so that

$$x(t) \approx \sum_{i=1}^{r} b_i \phi_i e^{\omega_i t}.$$

This representation yields oscillation frequencies, growth or decay rates, and spatial modes in a compact spectral model (Conti et al., 2023, Liu et al., 5 Aug 2025, Corrochano et al., 2022).

HODMD replaces the one-step closure of classical DMD by a $d$-step memory. In the Hankel or delay-embedding view, one forms

$$X = [x_1, x_2, \dots, x_{m-1}], \qquad Y = [x_2, x_3, \dots, x_m].$$0

builds shifted augmented matrices from the $$X = [x_1, x_2, \dots, x_{m-1}], \qquad Y = [x_2, x_3, \dots, x_m].$$1, and then applies the standard DMD pipeline in the augmented space. In the autoregressive view, HODMD assumes

$$X = [x_1, x_2, \dots, x_{m-1}], \qquad Y = [x_2, x_3, \dots, x_m].$$2

which can be written as a first-order companion evolution for the stacked state. Setting $$X = [x_1, x_2, \dots, x_{m-1}], \qquad Y = [x_2, x_3, \dots, x_m].$$3 recovers standard DMD. The central benefit is that time-delay embedding increases the effective dimensionality from $$X = [x_1, x_2, \dots, x_{m-1}], \qquad Y = [x_2, x_3, \dots, x_m].$$4 to $$X = [x_1, x_2, \dots, x_{m-1}], \qquad Y = [x_2, x_3, \dots, x_m].$$5, allowing more oscillatory content to be resolved when few sensors are available or when standard DMD is underdetermined (Conti et al., 2023, Liu et al., 5 Aug 2025, Corrochano et al., 2022).

This “higher order” designation refers to temporal memory, not to a fundamentally different spectral objective. HODMD still produces modes with exponential time dependence, but it estimates their eigenvalues and amplitudes in an augmented coordinate system that is better conditioned for noisy, multiscale, or quasi-periodic data.

2. Equivalent formulations and theoretical generalizations

The most common implementation of HODMD is the block-companion formulation. If

$$X = [x_1, x_2, \dots, x_{m-1}], \qquad Y = [x_2, x_3, \dots, x_m].$$6

then

$$X = [x_1, x_2, \dots, x_{m-1}], \qquad Y = [x_2, x_3, \dots, x_m].$$7

with $$X = [x_1, x_2, \dots, x_{m-1}], \qquad Y = [x_2, x_3, \dots, x_m].$$8 a companion operator whose upper blocks shift the delay window and whose last block row encodes the $$X = [x_1, x_2, \dots, x_{m-1}], \qquad Y = [x_2, x_3, \dots, x_m].$$9-step recurrence. The eigenpairs of this augmented operator determine the HODMD modes and their discrete-time eigenvalues, and the original-state reconstruction is obtained by selecting the first block of each augmented mode. This equivalence between delayed embedding and higher-order autoregression explains why HODMD is often described simultaneously as a Koopman-inspired delay method and as a data-driven AR($X = U \Sigma V^*,$0) model (Liu et al., 5 Aug 2025, Lazpita et al., 29 Jul 2025).

Several variants extend this core construction. In multidimensional HODMD, or mdHODMD, the raw field data are first compressed by HOSVD and HODMD is then applied to the temporal coefficients rather than directly to the full tensor. In intraventricular-flow modelling, this tensor-based formulation was combined with spectral constraints that retained harmonics near the cardiac fundamental and pruned modes using a growth-rate threshold, after which the retained modes were enforced to have $X = U \Sigma V^*,$1 for long-horizon stability (Lazpita et al., 29 Jul 2025). In reacting-flow analysis, multi-dimensional HODMD similarly used HOSVD to exploit the variable-by-space-by-time structure of combustion datasets before applying the delayed Koopman model (Corrochano et al., 2022).

A different line of theory embeds HODMD into constrained Koopman regression. Invariant Consistent DMD shows that HODMD can be recovered by imposing linear functional constraints that encode the block-shift relationships among delayed observables, thereby placing higher-order and affine DMD in a common constrained framework (Seenivasaharagavan et al., 2023). An even more abstract formulation analyzes higher-order dynamical systems directly through higher-order Liouville operators acting on signal-valued RKHSs. In that setting, systems such as $X = U \Sigma V^*,$2 are handled without augmenting the measured output by numerical differentiation of the full trajectory, and the resulting computations remain structurally close to first-order occupation-kernel DMD (Rosenfeld et al., 2021).

These formulations suggest that HODMD is not a single algorithmic variant but a family of closely related spectral estimators built around the same principle: enrich temporal observability by incorporating delayed coordinates or higher-order evolution constraints.

3. Computational workflow and parameterization

Most practical HODMD pipelines share a common sequence: preprocess the data, construct delayed snapshots, reduce dimension by SVD or HOSVD, estimate the reduced operator, compute eigenpairs, select modes, fit amplitudes by least squares, and reconstruct or forecast. What varies across applications is how aggressively the preprocessing, truncation, and mode-selection stages are tuned (Conti et al., 2023, Liu et al., 5 Aug 2025, Bell-Navas et al., 2024, Corrochano et al., 2023).

Preprocessing is often decisive. In the underground-farm study, occasional missing samples were linearly interpolated and the signals were mean-centered; the authors report that standardization to zero mean significantly improved reconstruction and forecasts (Conti et al., 2023). In echocardiography pipelines, cropping, grayscale conversion, temporal homogenization, and sequence splitting to preserve uniform $X = U \Sigma V^*,$3 were required before HODMD was applied to vectorized frames or to image tensors (Bell-Navas et al., 2024). In reacting-flow databases, centering and scaling were essential because temperature and species concentrations differed by many orders of magnitude; range scaling and Pareto-type scalings changed both clustering behaviour and reconstruction accuracy (Corrochano et al., 2023).

The delay order $X = U \Sigma V^*,$4 is usually selected empirically. Underground-farm forecasting used a grid search on three-day RRMSE and selected $X = U \Sigma V^*,$5, $X = U \Sigma V^*,$6, and $X = U \Sigma V^*,$7 for three monthly scenarios sampled every ten minutes (Conti et al., 2023). Circuit MOR used the embedding size $X = U \Sigma V^*,$8 and recommended the empirical rule $X = U \Sigma V^*,$9 to ensure sufficient augmented dimension, with errors on the IBMPG1t benchmark dropping from $\tilde{A} = U^* Y V \Sigma^{-1},$0 at $\tilde{A} = U^* Y V \Sigma^{-1},$1 to $\tilde{A} = U^* Y V \Sigma^{-1},$2 at $\tilde{A} = U^* Y V \Sigma^{-1},$3 (Liu et al., 5 Aug 2025). In echocardiography, sequence-dependent choices such as $\tilde{A} = U^* Y V \Sigma^{-1},$4, $\tilde{A} = U^* Y V \Sigma^{-1},$5, or $\tilde{A} = U^* Y V \Sigma^{-1},$6 were used to stabilize cardiac and respiratory frequencies (Bell-Navas et al., 2024). In power systems, Li, Wen, and Schäfer selected $\tilde{A} = U^* Y V \Sigma^{-1},$7 for IEEE 14-bus and $\tilde{A} = U^* Y V \Sigma^{-1},$8 for WECC after inspecting reconstruction error and the recovery of oscillatory decay (Li et al., 10 Feb 2025).

Rank truncation is handled by singular-value thresholds, fixed retained ranks, or sparsity-promoting DMD. Circuit MOR retained enough singular values to capture at least $\tilde{A} = U^* Y V \Sigma^{-1},$9 energy (Liu et al., 5 Aug 2025). Mechanical-vibration analysis recommends optimal hard-thresholding and, in noisy measurements, TLS-type operator estimation instead of ordinary least squares (Tuor et al., 2023). In underground-farm modelling, forward-backward DMD was used to reduce bias due to sensor noise, while mode selection relied on an integral contribution criterion rather than on raw amplitude alone, excluding numerically transient modes with very large $\tilde{A} W = W \Lambda,$0 (Conti et al., 2023).

Mode selection is application-dependent. Some studies retain modes by amplitude tolerance, others by cumulative contribution, spectral windows, or cluster-specific reconstruction error. In h-HODMD, variables are iteratively grouped according to which cluster-specific HODMD reconstruction best fits them, producing dedicated mode sets for chemically or dynamically coherent subsets rather than a single global basis (Corrochano et al., 2023). In spectrally constrained mdHODMD for ventricular flow, modes were pruned by growth-rate threshold $\tilde{A} W = W \Lambda,$1 and by closeness to harmonic bands near $\tilde{A} W = W \Lambda,$2 (Lazpita et al., 29 Jul 2025).

4. Empirical uses and reported performance

Recent applications indicate that HODMD is used both as a reduced-order predictor and as a modal-analysis tool, with quantitative gains reported in forecasting stability, reconstruction fidelity, speedup, and feature quality across domains (Conti et al., 2023, Liu et al., 5 Aug 2025, Groun et al., 2024, Lazpita et al., 29 Jul 2025, Li et al., 10 Feb 2025, Mendez et al., 2023, Corrochano et al., 2022, Tuor et al., 2023, Liu, 5 Mar 2025).

Domain	Role of HODMD	Reported outcome
Underground agricultural farm	Interpretable reduced-order digital-twin model for VPD	Three physically interpretable mode pairs; three-day forecasts; average RMSE $\tilde{A} W = W \Lambda,$3 and $\tilde{A} W = W \Lambda,$4 in typical scenarios
Large-scale circuits	Data-driven model order reduction	Speedups of $\tilde{A} W = W \Lambda,$5–$\tilde{A} W = W \Lambda,$6 with sub-percent $\tilde{A} W = W \Lambda,$7 errors
Echocardiography classification	Feature extraction and data augmentation	Accuracy improvement of up to $\tilde{A} W = W \Lambda,$8 percentage points on unseen data
Intraventricular flow	mdHODMD reduced-order model with spectral constraints	Errors below $\tilde{A} W = W \Lambda,$9 for one geometry and below $\Phi = Y V \Sigma^{-1} W.$0 for the second; speed-up factor of at least $\Phi = Y V \Sigma^{-1} W.$1
Power-system frequency dynamics	Reconstruction, forecasting, and modal identification	Test RRMSE $\Phi = Y V \Sigma^{-1} W.$2 on IEEE 14-bus and $\Phi = Y V \Sigma^{-1} W.$3 on WECC; local and global modes identified
Mechanical vibrations	High-resolution spectral decomposition	Frequency error below $\Phi = Y V \Sigma^{-1} W.$4 Hz in simulated close-mode tests; sharp KDS spectra beyond FFT limits

In building and environmental modelling, HODMD identified three recurring conjugate mode pairs in an underground hydroponic farm: a mean-trend pair at $\Phi = Y V \Sigma^{-1} W.$5 Hz, a lighting-cycle pair at $\Phi = Y V \Sigma^{-1} W.$6 Hz, and a miscellaneous operational pair at $\Phi = Y V \Sigma^{-1} W.$7 Hz. Six dominant modes sufficed for two typical seasonal scenarios, whereas twelve were retained in a heatwave case with stronger transients (Conti et al., 2023). In bubble-column hydrodynamics, HODMD consistently recovered low-frequency plume oscillation modes near $\Phi = Y V \Sigma^{-1} W.$8–$\Phi = Y V \Sigma^{-1} W.$9 Hz and higher-frequency content near $\lambda_i$0–$\lambda_i$1 Hz, and six measurement points in top, centre, and bottom regions were sufficient to recover the dominant frequency across all tested superficial velocities (Mendez et al., 2023).

In large-scale circuits, HODMD was reported to overcome reconstruction failures caused by insufficient spatial resolution. For the linear transmission line, the reported $\lambda_i$2 error was $\lambda_i$3 with an $\lambda_i$4 speedup; for the IBMPG1t power-grid benchmark the method achieved $\lambda_i$5 error, a maximum amplitude difference below $\lambda_i$6 V, and a $\lambda_i$7 speedup; for IBMPG6t, the reported $\lambda_i$8 error was $\lambda_i$9 with a $$\omega_i = \frac{\log(\lambda_i)}{\Delta t} = \alpha_i + i 2\pi f_i,$$0 speedup (Liu et al., 5 Aug 2025).

In cardiovascular flow ROMs, mdHODMD learned from as few as the first three simulated cycles and predicted ten cycles beyond the training window. The reported reconstruction and prediction errors remained below $$\omega_i = \frac{\log(\lambda_i)}{\Delta t} = \alpha_i + i 2\pi f_i,$$1 for one idealized ventricle and below $$\omega_i = \frac{\log(\lambda_i)}{\Delta t} = \alpha_i + i 2\pi f_i,$$2 for another, while the speed-up factor relative to full CFD was approximately $$\omega_i = \frac{\log(\lambda_i)}{\Delta t} = \alpha_i + i 2\pi f_i,$$3 and $$\omega_i = \frac{\log(\lambda_i)}{\Delta t} = \alpha_i + i 2\pi f_i,$$4 for the two geometries (Lazpita et al., 29 Jul 2025). In power systems, Li, Wen, and Schäfer reported that HODMD identified both local and inter-area oscillations in the WECC system, including local modes at $$\omega_i = \frac{\log(\lambda_i)}{\Delta t} = \alpha_i + i 2\pi f_i,$$5 Hz and $$\omega_i = \frac{\log(\lambda_i)}{\Delta t} = \alpha_i + i 2\pi f_i,$$6 Hz and global modes at $$\omega_i = \frac{\log(\lambda_i)}{\Delta t} = \alpha_i + i 2\pi f_i,$$7 Hz and $$\omega_i = \frac{\log(\lambda_i)}{\Delta t} = \alpha_i + i 2\pi f_i,$$8 Hz, with damping percentages of $$\omega_i = \frac{\log(\lambda_i)}{\Delta t} = \alpha_i + i 2\pi f_i,$$9, $$x(t) \approx \sum_{i=1}^{r} b_i \phi_i e^{\omega_i t}.$$0, $$x(t) \approx \sum_{i=1}^{r} b_i \phi_i e^{\omega_i t}.$$1, and $$x(t) \approx \sum_{i=1}^{r} b_i \phi_i e^{\omega_i t}.$$2 for selected modes (Li et al., 10 Feb 2025).

Medical-imaging work uses HODMD differently. In one study, dominant HODMD modes from echocardiography videos were turned into images and appended to the training database, increasing unseen four-class accuracy from $$x(t) \approx \sum_{i=1}^{r} b_i \phi_i e^{\omega_i t}.$$3 to $$x(t) \approx \sum_{i=1}^{r} b_i \phi_i e^{\omega_i t}.$$4 and five-class unseen accuracy from $$x(t) \approx \sum_{i=1}^{r} b_i \phi_i e^{\omega_i t}.$$5 to $$x(t) \approx \sum_{i=1}^{r} b_i \phi_i e^{\omega_i t}.$$6 (Groun et al., 2024). In a related ViT pipeline for small datasets, HODMD-derived reconstructions and mode images improved sequence-level classification, with the best configuration reaching per-sequence accuracy $$x(t) \approx \sum_{i=1}^{r} b_i \phi_i e^{\omega_i t}.$$7 (Bell-Navas et al., 2024).

5. Relation to neighboring methods and common misconceptions

HODMD is often contrasted first with standard DMD. The distinction is not merely that HODMD produces more modes; rather, it changes the regression problem by fitting the evolution in a delay-embedded space. This is why HODMD can recover frequencies that standard DMD misses when the sensor network is sparse or when the dominant dynamics are distributed across several time scales (Conti et al., 2023, Liu et al., 5 Aug 2025, Li et al., 10 Feb 2025).

It is also frequently contrasted with classical spectral tools such as FFT, STFT, and PSD. In bubble columns and mechanical vibrations, the reported limitation of FFT-based methods is not that they fail universally, but that they are local, point-dependent, and constrained by time-bandwidth trade-offs. HODMD instead yields global coherent modes together with growth or decay rates and can reconstruct damped oscillations with high spectral resolution from multi-point or full-field data (Mendez et al., 2023, Tuor et al., 2023). The kernel density spectrum used in vibration analysis is a post-processing device for HODMD poles, not an alternative dynamical model (Tuor et al., 2023).

Compared with black-box predictors, HODMD retains an explicitly interpretable modal structure. The underground-farm study positions it against SARIMA and ANN-type models, arguing that HODMD is more stable under atypical conditions because it filters noise by emphasizing coherent global dynamics instead of fitting local fluctuations (Conti et al., 2023). In circuits, HODMD is contrasted with projection-based MOR methods such as PRIMA, Krylov, and balanced truncation because it does not require explicit circuit matrices after snapshot acquisition and does not impose constraints on circuit topology or component type (Liu et al., 5 Aug 2025). In power systems, HODMD outperformed standard DMD and also outperformed or matched nonlinear baselines such as LANDO and SINDy on the reported benchmarks (Li et al., 10 Feb 2025).

A recurring misconception is that HODMD is intrinsically nonlinear. The evidence across these studies suggests a more precise statement: HODMD uses a linear superposition in delayed or lifted coordinates and therefore still inherits Koopman-style linearization assumptions. This is why papers repeatedly list limitations such as sensitivity to parameter selection, extrapolation risk outside the training regime, and degraded performance under strong nonstationarity or abrupt regime changes (Conti et al., 2023, Lazpita et al., 29 Jul 2025, Groun et al., 2024). Another misconception, especially in medical applications, is to treat HODMD as a classifier. In the echocardiography studies, HODMD is a feature extractor and augmentation mechanism; the actual classifiers are a CNN or a ViT trained on HODMD-derived images (Groun et al., 2024, Bell-Navas et al., 2024).

6. Limitations, variants, and current directions

The limitations reported across domains are consistent. HODMD assumes a modal expansion of the form

$$x(t) \approx \sum_{i=1}^{r} b_i \phi_i e^{\omega_i t}.$$8

so model quality depends on how well the underlying dynamics can be represented by delayed linear recurrence. External forcings, abrupt transients, moving boundaries, sparse sensor placement, and regime changes can all force the inclusion of additional modes, reduce interpretability, or destabilize extrapolation (Conti et al., 2023, Lazpita et al., 29 Jul 2025, Liu, 5 Mar 2025). The memory footprint also grows with delay order because the stacked matrices scale with the augmented dimension, a practical issue emphasized in circuit simulation and time-parallel PDE solvers (Liu et al., 5 Aug 2025, Liu, 5 Mar 2025).

Several variants attempt to address these issues. h-HODMD replaces a single global model by iteratively clustered models and, in reacting-flow databases, improved reconstruction for most variables relative to standard HODMD; under Range scaling, improvements of approximately $$x(t) \approx \sum_{i=1}^{r} b_i \phi_i e^{\omega_i t}.$$9 in RRMSE were reported for key radicals such as OH, O, H$d$0O$d$1, and CH$d$2 (Corrochano et al., 2023). Spectrally constrained mdHODMD uses HOSVD denoising, harmonic filtering, and explicit control of growth rates to stabilize long-horizon predictions in periodic cardiac flows (Lazpita et al., 29 Jul 2025). Invariant Consistent DMD incorporates linear constraints that preserve functional relationships among observables and consistency along geometric invariants, recovering HODMD as a special case while also enabling eigenfunctions at eigenvalue $d$3 associated with invariant sets (Seenivasaharagavan et al., 2023).

Open directions are explicitly identified in several studies. For non-periodic cardiac dynamics, proposed remedies include adaptive windows, time-varying operators, and piecewise HODMD across phases (Lazpita et al., 29 Jul 2025). For deployment in large-scale time-parallel simulation, randomized SVD and parallelized HODMD steps are proposed to reduce the SVD bottleneck (Liu, 5 Mar 2025). In echocardiography, future work points toward larger datasets, view normalization, and the use of numerical spectral features alongside images (Bell-Navas et al., 2024). A plausible implication is that current research is shifting from HODMD as a stand-alone modal estimator toward HODMD as a modular component inside broader pipelines for ROM, control, feature learning, and structure-preserving Koopman regression.