Hankel-DMD

Time-delay embedding for scalar and low-dimensional time series

Overview

Hankel-DMD uses time-delay embedding to construct a Hankel matrix from scalar or low-dimensional time series. By stacking delayed copies of the signal as rows of a matrix, the method creates a Krylov-like subspace that captures temporal structure even from a single measurement channel. This lifts a one-dimensional observation into a higher-dimensional space where standard DMD can extract dynamic modes and frequencies.

The approach is grounded in Takens' embedding theorem, which guarantees that for a sufficient number of delays, the delay-coordinate map preserves the topology of the underlying attractor. In the Koopman framework, the Hankel matrix spans an approximation to a Koopman-invariant subspace, enabling spectral analysis of the operator from minimal measurement data.

Hankel Matrix Construction

Given a scalar time series [x_0, x_1, x_2, ..., x_{n-1}] with n samples, the Hankel matrix H is formed by arranging d delayed copies of the signal into rows. Each row is a window of consecutive samples shifted forward by one time step:

The matrix has d rows (the number of delays) and n - d + 1 columns. Each column is a delay vector representing the state of the system at a particular time, embedded in a d-dimensional space. Standard DMD is then applied to the column pairs of H to extract eigenvalues and modes.

HankelConfig

The HankelConfig struct controls how the Hankel matrix is built and how the subsequent DMD decomposition is performed.

Basic Usage

The following example creates a scalar sinusoidal signal, runs Hankel-DMD, and inspects the resulting eigenvalues. For a pure oscillation the eigenvalues should lie on (or very near) the unit circle in the complex plane.

use koopman_dmd::{HankelConfig, hankel_dmd};
use std::f64::consts::PI;

fn main() {
    // Generate a scalar oscillating signal: sin(2*pi*0.25*t)
    let n = 200;
    let dt = 0.1;
    let freq_true = 0.25;
    let signal: Vec<f64> = (0..n)
        .map(|i| (2.0 * PI * freq_true * i as f64 * dt).sin())
        .collect();

    // Configure Hankel-DMD
    let config = HankelConfig {
        delays: None,    // auto: n / 3
        rank: None,      // auto truncation
        dt,
    };

    // Run Hankel-DMD
    let result = hankel_dmd(&signal, &config).unwrap();

    // Inspect eigenvalues -- magnitudes near 1.0 for pure oscillation
    for eig in result.eigenvalues() {
        let mag = eig.norm();
        let freq = eig.arg() / (2.0 * PI * dt);
        println!("|lambda| = {:.4}, freq = {:.4} Hz", mag, freq);
    }
}

import numpy as np
import koopman_dmd as kdmd

# Generate a scalar oscillating signal: sin(2*pi*0.25*t)
n = 200
dt = 0.1
freq_true = 0.25
t = np.arange(n) * dt
signal = np.sin(2.0 * np.pi * freq_true * t)

# Configure and run Hankel-DMD
result = kdmd.hankel_dmd(signal, delays=None, rank=None, dt=dt)

# Inspect eigenvalues -- magnitudes near 1.0 for pure oscillation
for eig in result.eigenvalues():
    mag = np.abs(eig)
    freq = np.angle(eig) / (2.0 * np.pi * dt)
    print(f"|lambda| = {mag:.4f}, freq = {freq:.4f} Hz")

library(koopman.dmd)

# Generate a scalar oscillating signal: sin(2*pi*0.25*t)
n <- 200
dt <- 0.1
freq_true <- 0.25
t <- (0:(n - 1)) * dt
signal <- sin(2 * pi * freq_true * t)

# Configure and run Hankel-DMD
result <- hankel_dmd(signal, delays = NULL, rank = NULL, dt = dt)

# Inspect eigenvalues -- magnitudes near 1.0 for pure oscillation
eigs <- result$eigenvalues()
for (i in seq_along(eigs)) {
  mag <- Mod(eigs[i])
  freq <- Arg(eigs[i]) / (2 * pi * dt)
  cat(sprintf("|lambda| = %.4f, freq = %.4f Hz\n", mag, freq))
}

Reconstruction

After fitting a Hankel-DMD model, use hankel_reconstruct to obtain fitted values projected back into the original scalar space. Use hankel_predict to forecast future values beyond the training window.

// Reconstruct the fitted signal in original space
let fitted = hankel_reconstruct(&result).unwrap();
println!("Reconstruction length: {}", fitted.len());

// Forecast 50 steps into the future
let forecast = hankel_predict(&result, 50).unwrap();
println!("Forecast length: {}", forecast.len());

// Compute reconstruction error
let error: f64 = signal.iter()
    .zip(fitted.iter())
    .map(|(s, f)| (s - f).powi(2))
    .sum::<f64>()
    .sqrt() / signal.len() as f64;
println!("RMS reconstruction error: {:.6}", error);

# Reconstruct the fitted signal in original space
fitted = kdmd.hankel_reconstruct(result)
print(f"Reconstruction length: {len(fitted)}")

# Forecast 50 steps into the future
forecast = kdmd.hankel_predict(result, n_steps=50)
print(f"Forecast length: {len(forecast)}")

# Compute reconstruction error
error = np.sqrt(np.mean((signal - fitted) ** 2))
print(f"RMS reconstruction error: {error:.6f}")

# Reconstruct the fitted signal in original space
fitted <- hankel_reconstruct(result)
cat("Reconstruction length:", length(fitted), "\n")

# Forecast 50 steps into the future
forecast <- hankel_predict(result, n_steps = 50)
cat("Forecast length:", length(forecast), "\n")

# Compute reconstruction error
error <- sqrt(mean((signal - fitted)^2))
cat(sprintf("RMS reconstruction error: %.6f\n", error))

Frequency Recovery

A key application of Hankel-DMD is recovering oscillation frequencies from noisy scalar measurements. The following example embeds a known frequency into a noisy signal and demonstrates that Hankel-DMD recovers it accurately.

use koopman_dmd::{HankelConfig, hankel_dmd};
use std::f64::consts::PI;
use rand::distributions::{Distribution, Normal};

fn main() {
    let n = 500;
    let dt = 0.05;
    let freq_true = 1.5;  // 1.5 Hz signal
    let noise_std = 0.3;

    let mut rng = rand::thread_rng();
    let normal = Normal::new(0.0, noise_std);

    let signal: Vec<f64> = (0..n)
        .map(|i| {
            let t = i as f64 * dt;
            (2.0 * PI * freq_true * t).sin() + normal.sample(&mut rng)
        })
        .collect();

    let config = HankelConfig { delays: None, rank: Some(4), dt };
    let result = hankel_dmd(&signal, &config).unwrap();

    // Find the dominant frequency (largest-magnitude eigenvalue)
    let eigs = result.eigenvalues();
    let dominant = eigs.iter()
        .max_by(|a, b| a.norm().partial_cmp(&b.norm()).unwrap())
        .unwrap();
    let freq_recovered = dominant.arg().abs() / (2.0 * PI * dt);

    println!("True frequency:      {:.4} Hz", freq_true);
    println!("Recovered frequency: {:.4} Hz", freq_recovered);
}

import numpy as np
import koopman_dmd as kdmd

n = 500
dt = 0.05
freq_true = 1.5  # 1.5 Hz signal
noise_std = 0.3

t = np.arange(n) * dt
signal = np.sin(2.0 * np.pi * freq_true * t) + noise_std * np.random.randn(n)

result = kdmd.hankel_dmd(signal, rank=4, dt=dt)

# Find the dominant frequency (largest-magnitude eigenvalue)
eigs = result.eigenvalues()
dominant = eigs[np.argmax(np.abs(eigs))]
freq_recovered = np.abs(np.angle(dominant)) / (2.0 * np.pi * dt)

print(f"True frequency:      {freq_true:.4f} Hz")
print(f"Recovered frequency: {freq_recovered:.4f} Hz")

library(koopman.dmd)

n <- 500
dt <- 0.05
freq_true <- 1.5  # 1.5 Hz signal
noise_std <- 0.3

t <- (0:(n - 1)) * dt
signal <- sin(2 * pi * freq_true * t) + noise_std * rnorm(n)

result <- hankel_dmd(signal, rank = 4, dt = dt)

# Find the dominant frequency (largest-magnitude eigenvalue)
eigs <- result$eigenvalues()
dominant <- eigs[which.max(Mod(eigs))]
freq_recovered <- abs(Arg(dominant)) / (2 * pi * dt)

cat(sprintf("True frequency:      %.4f Hz\n", freq_true))
cat(sprintf("Recovered frequency: %.4f Hz\n", freq_recovered))

When to Use Hankel-DMD

Hankel-DMD is the right tool when your data does not naturally come as multi-variable snapshot matrices. The following table summarizes the key trade-offs between standard DMD and Hankel-DMD.

Use Hankel-DMD when you have a scalar time series and want to extract oscillation frequencies, when your measurement system provides a limited number of sensors, or when you want to apply DMD-based forecasting to one-dimensional data. For systems where you already have high-dimensional snapshot data from multiple spatial locations, standard DMD or Extended DMD will typically be more efficient and direct.