Dynamic Mode Decomposition (DMD)

DMD is a purely data-driven spectral method (Schmid, 2010). It decomposes an IBI time series into dynamic modes and temporal eigenvalues, then uses them to reconstruct or forecast the signal. No mechanistic model of the heart is assumed.

The implementation is the mean-centered, unit-circle-projected, time-delay (Hankel) DMD with a global least-squares amplitude solve ("V2 / UnitCircleDMD").

Theory

The fitting pipeline has five steps:

1. Mean centering

The temporal mean $\mu$ of the series is subtracted before the embedding (affine / constant-mode handling) and added back, never rescaled, at reconstruction.

2. Hankel Embedding

The centered scalar IBI series $\mathbf{x} = [x_1, x_2, \ldots, x_n]$ is lifted into a data matrix via delay embedding with a fixed delay dimension $d$ (default 50, clamped to $\lfloor n/2 \rfloor$ for short series):

\[X = \begin{bmatrix} x_1 & x_2 & \cdots & x_{n-d+1} \\ x_2 & x_3 & \cdots & x_{n-d+2} \\ \vdots & & & \vdots \\ x_d & x_{d+1}& \cdots & x_n \end{bmatrix}\]

3. DMD Decomposition with energy-based rank selection

The data matrix is split into shifted snapshots $X_1 = X[:, 1\!:\!end\!-\!1]$ and $X_2 = X[:, 2\!:\!end]$. Exact DMD finds the best-fit linear operator $A$ such that $X_2 \approx A X_1$ via an SVD of $X_1$:

\[X_1 \approx U_r \Sigma_r V_r^T\]

The retained rank $r$ is not fixed by hand: it is the smallest rank whose squared singular values capture the energy fraction (default 0.99) of the total, capped at rank (the rank cap $r_{\max}$, default 20). The reduced operator $\tilde{A} = U_r^T X_2 V_r \Sigma_r^{-1}$ is diagonalized to obtain eigenvalues $\lambda_i$ and eigenvectors $W$. The dynamic modes are:

\[\Phi = X_2 V_r \Sigma_r^{-1} W\]

4. Unit-circle eigenvalue projection

Each retained eigenvalue is projected onto the unit circle ($\lambda \to \lambda / |\lambda|$), preserving its frequency while making the free run marginally stable: oscillations are sustained rather than decaying or exploding.

5. Global least-squares amplitudes and reconstruction

The mode amplitudes $\mathbf{b}$ are fitted by least squares against the full first-row Hankel trajectory (not a single initial snapshot), so the reconstruction variance tracks the data. The signal is then reconstructed as:

\[x(t) \approx \mu + \operatorname{Re}\!\left(\sum_{i=1}^{r} b_i \, \lambda_i^{t-1}\right)\]

The physiological 300–2000 ms clip is applied only as an output guard in simulate: the mean is always the data mean $\mu$.

Honest limit

RR tachograms are broadband / aperiodic. A low-rank linear model recovers the mean and the dominant LF oscillation, but cannot reproduce the full broadband variance without over-shooting: the corresponding variance check in the test suite is an explicit @test_broken. DMD is deliberately kept as the weaker baseline generative model in HeartRateLab; its real strength for RR is short-horizon prediction via forecast.

Parameters

Parameter	Default	Description
`rank`	20	Rank cap $r_{\max}$. The retained rank $r$ is chosen by the energy threshold, then clamped to at most `rank`.
`d`	50	Delay / embedding dimension of the Hankel matrix (clamped to $\lfloor n/2 \rfloor$). Use ≈100 for faster-mean records.
`energy`	0.99	SVD energy fraction used to pick the retained rank $r$.

DMD has no parameter_space: it is entirely data-driven and has no priors to sample. rank, d, and energy are structural knobs of the decomposition, not fitted physiological parameters.

Workflow

DMD must be fitted before simulating. fit decomposes the training data; simulate reconstructs using the learned modes.

dmd = DMD()   →   result = fit(dmd, ibis)   →   simulate(result.model, nothing, n)

Note that simulate is called on result.model (the fitted DMD with the mean, modes, eigenvalues, and amplitudes filled in), not on the original dmd object.

Examples

Fit and Reconstruct

using HeartRateLab

ibis = read_txt("data.txt")

dmd = DMD()   # rank cap 20, d = 50, energy = 0.99
result = fit(dmd, ibis)

println("Retained rank r:       ", result.diagnostics["rank"])
println("Embedding dimension:   ", result.diagnostics["embedding_dimension"])
println("Fitted mean (ms):      ", round(result.diagnostics["mean"]; digits=1))

# Reconstruct the training signal
reconstructed = simulate(result.model, nothing, length(ibis))

Effect of the Rank Cap on Reconstruction Quality

using HeartRateLab

ibis = read_txt("data.txt")

for rmax in [2, 5, 10, 20]
    result = fit(DMD(rank=rmax), ibis)
    recon  = simulate(result.model, nothing, length(ibis))
    rmse   = sqrt(sum((ibis .- recon).^2) / length(ibis))
    println("rank cap=", rmax, "  retained r=", result.diagnostics["rank"],
            "  RMSE=", round(rmse; digits=1), " ms")
end

Short-Horizon Forecasting

using HeartRateLab

hist = read_txt("data.txt")[1:2000]
ŷ = forecast(DMD(d=50), hist, 5)   # next 5 beats (richer rank by default)

API Reference

HeartRateLab.Models.DMD — Type

DMD <: AbstractHRVModel

Dynamic Mode Decomposition model for RR-interval (IBI) dynamics.

This is the V2 / UnitCircleDMD formulation validated in docs/dmd-rr-modeling-research.md: a mean-centered, unit-circle-projected, time-delay (Hankel) DMD with a global least-squares amplitude solve. It replaces the earlier intentionally-weak model that forced the mean to ~800 ms and collapsed the variance.

The pipeline:

Center. Subtract the temporal mean μ of the data before the embedding (affine / constant-mode handling, Hirsh et al. 2019). The mean is added back at reconstruction — it is never rescaled.
Embed. Build a d-delay Hankel matrix of the centered signal so DMD has a state space to act on (Takens / Hankel-DMD, Brunton et al. 2017).
Decompose. Exact DMD on the Hankel snapshots; retain r modes by an SVD energy threshold (default 0.99) capped at rmax.
Stabilize. Project every retained eigenvalue onto the unit circle (λ → λ/|λ|), preserving each mode's frequency while making the free run marginally stable (sustained, neither decaying nor exploding) — the mean-subtraction≈DFT regime (Chen/Tu/Rowley 2012).
Global amplitudes. Fit the mode amplitudes b by least squares against the full first-row Hankel trajectory (the optDMD-style global amplitude step), so reconstruction variance matches the data instead of being pinned by a single initial column.

The physiological 300–2000 ms clip is applied only as an output guard in simulate — never as a mean rescale (that was the core bug of the old model).

Honest limit

RR tachograms are broadband / aperiodic. A low-rank linear model recovers the mean and the dominant LF oscillation but cannot reproduce the full broadband variance without over-shooting. DMD is a strong short-horizon forecaster (see forecast) and a fair generative model for RR.

Parameters / fields

rank : rank cap rmax (kept as the public knob; DMD(rank=…) works).
d : delay / embedding dimension (default 50; ≈100 for faster-mean records, e.g. NSRDB 16265 at ~620 ms).
energy: SVD energy fraction used to pick the retained rank r (default 0.99).
μ : fitted temporal mean (0 until fitted).
modes : exact-DMD modes Φ (d × r, empty until fitted).
evals : retained, unit-circle-projected eigenvalues λ (empty until fitted).
b : global least-squares amplitudes (empty until fitted).
r : number of retained modes actually used (0 until fitted).

source

HeartRateLab.Models.fit — Method

fit(model::DMD, data::Vector{Float64}; kwargs...) -> ModelFitResult

Fit the V2 Hankel-DMD operator to IBI data: center on the data mean, build the model.d-delay Hankel matrix, run exact DMD, retain r modes by the energy threshold (model.energy, capped at model.rank), project the eigenvalues onto the unit circle, and solve for global least-squares amplitudes.

Returns a ModelFitResult whose model is a fitted DMD carrying the temporal mean, modes, unit-circle eigenvalues and amplitudes. DMD stores no continuous params, so params is an empty NamedTuple and the method is reported as :gradient (deterministic, like the other data-driven fits).

source

HeartRateLab.Models.simulate — Method

simulate(model::DMD, params::Union{NamedTuple,Nothing}, n_beats::Int) -> Vector{Float64}

Free-running reconstruction from a fitted V2 DMD model:

x_c(t) = Re( Σ_j b_j · λ_j^{t-1} ),    x(t) = x_c(t) + μ.

The amplitudes b were fit against the basis Ψ[t,j] = φ1_j · λ_j^{t-1}, so they already fold in the first-row mode weight φ1. params is ignored (DMD is data-driven). The physiological 300–2000 ms clip is applied as an output guard only; the mean is the data mean μ, not a forced constant.

source

HeartRateLab.Models.forecast — Method

forecast(model::DMD, history::AbstractVector, h::Int; d=model.d, rank=max(model.rank, 40), energy=1.0) -> Vector{Float64}

h-step-ahead forecast of an IBI series from its history.

This is DMD's real strength for RR: with a richer rank, the one-step out-of-sample NRMSE ≈ 0.19, beating both the persistence (≈0.35) and mean (≈1.04) baselines at every horizon (see docs/dmd-rr-modeling-research.md §5). A fresh DMD operator is fit on history (centered + Hankel-embedded). Unlike the generative simulate path, the eigenvalues are kept at their raw magnitudes — predictive accuracy, not free-running stability, is the goal.

The most recent d-beat delay vector seeds the rank-r Hankel state; the exact-DMD propagator A = Uᵣ Ã Uᵣ' is applied h times and the newest Hankel coordinate (z[end] + μ) is read off at each step, giving the closed-loop multi-step prediction. The mean used is history's mean (added back).

For forecasting the report recommends a richer rank than for generation (r ≈ 20–40), so this defaults to energy=1.0 (use the full rank cap, ≥ 40) rather than the 0.99 energy threshold that the generative model uses.

Returns the h predicted beats following history (length h). Forecasting does not mutate model.

Example

hist = read_txt("example.txt")[1:2000]
ŷ = forecast(DMD(d=50), hist, 5)   # next 5 beats, richer rank by default

source