Extended DMD

Lifting functions for nonlinear dynamics

Overview

Standard DMD fits a linear model to observed data, which works well when the underlying dynamics are already close to linear. When the dynamics are genuinely nonlinear, the linear approximation can break down. Extended DMD addresses this by lifting the observables into a higher-dimensional space where the nonlinear dynamics become approximately linear under the Koopman operator.

Given state vector x, the lifted state is constructed as:

g(x) = [ x, g₁(x), g₂(x), ... ]

The original state variables are always retained as the first elements of the lifted vector. Additional terms — polynomials, trigonometric functions, or time-delayed copies — expand the observation space so that a linear operator in the lifted space can capture nonlinear relationships in the original space.

Key idea: Instead of finding a linear map A such that x(t+1) = A x(t), Extended DMD finds a linear map K such that g(x(t+1)) = K g(x(t)), where K operates on the richer lifted space.

Lifting Types

Polynomial Lifting

Polynomial lifting adds all monomial terms up to a specified degree d. For a two-dimensional state [x, y] with degree 2, the lifted vector becomes:

g([x, y]) = [ x, y, x², xy, y² ]

Higher degrees capture more complex nonlinearities but increase the dimensionality of the lifted space combinatorially. Degree 2 or 3 is typically sufficient for moderate nonlinearities.

// Rust: polynomial lifting configuration
let lifting = LiftingConfig::Polynomial { degree: 2 };

Trigonometric Lifting

Trigonometric lifting adds sine and cosine terms at integer harmonics of each state variable. With harmonics: h, it appends sin(k * x) and cos(k * x) for k = 1, 2, ..., h for each component of x. This is effective for systems with periodic or oscillatory structure.

// Rust: trigonometric lifting configuration
let lifting = LiftingConfig::Trigonometric { harmonics: 3 };

Delay Lifting

Delay lifting appends time-delayed copies of the state to each snapshot. With delays: d, snapshot t becomes [x(t), x(t-1), ..., x(t-d)]. This is closely related to Hankel DMD and is useful when the system has memory or when only partial state measurements are available.

// Rust: delay lifting configuration
let lifting = LiftingConfig::Delay { delays: 5 };

Lifting Type	Best For	Lifted Dimension Growth
`Polynomial`	Algebraic nonlinearities (quadratic drag, cubic stiffness)	Combinatorial in degree and state dimension
`Trigonometric`	Periodic or oscillatory dynamics	2 * harmonics * state dimension
`Delay`	Partial observations, systems with memory	delays * state dimension

Usage

The following example applies polynomial lifting to a nonlinear signal generated from sin²(t). After fitting, predictions are automatically projected back to the original observation space.

use koopman_dmd::{DmdConfig, LiftingConfig};
use ndarray::Array2;

fn main() {
    // Generate nonlinear signal: sin^2(t)
    let n = 200;
    let dt = 0.05;
    let t: Vec<f64> = (0..n).map(|i| i as f64 * dt).collect();
    let signal: Vec<f64> = t.iter().map(|&ti| ti.sin().powi(2)).collect();

    // Arrange into snapshot matrix (1 x n)
    let data = Array2::from_shape_vec((1, n), signal).unwrap();

    // Configure Extended DMD with polynomial lifting
    let config = DmdConfig {
        lifting: Some(LiftingConfig::Polynomial { degree: 2 }),
        dt: Some(dt),
        ..DmdConfig::default()
    };

    // Fit the model
    let result = koopman_dmd::fit(&data, &config).unwrap();

    // Predict -- output is in the original 1-d observation space
    let pred = result.predict(n);
    println!("Prediction shape: {:?}", pred.shape());
    // [1, 200] -- automatically back-projected
}

import numpy as np
from koopman_dmd import DmdConfig, LiftingConfig, fit

# Generate nonlinear signal: sin^2(t)
n = 200
dt = 0.05
t = np.arange(n) * dt
signal = np.sin(t) ** 2

# Arrange into snapshot matrix (1 x n)
data = signal.reshape(1, -1)

# Configure Extended DMD with polynomial lifting
config = DmdConfig(
    lifting=LiftingConfig.polynomial(degree=2),
    dt=dt,
)

# Fit and predict
result = fit(data, config)
pred = result.predict(n)

print(f"Prediction shape: {pred.shape}")
# (1, 200) -- automatically back-projected

library(koopmandmd)

# Generate nonlinear signal: sin^2(t)
n <- 200
dt <- 0.05
t <- seq(0, by = dt, length.out = n)
signal <- sin(t)^2

# Arrange into snapshot matrix (1 x n)
data <- matrix(signal, nrow = 1)

# Configure Extended DMD with polynomial lifting
config <- dmd_config(
  lifting = lifting_polynomial(degree = 2),
  dt = dt
)

# Fit and predict
result <- dmd_fit(data, config)
pred <- dmd_predict(result, n)

cat("Prediction shape:", dim(pred), "\n")
# 1 200 -- automatically back-projected

When to Use Lifting

Extended DMD with lifting functions is most effective when:

Known polynomial structure — the dynamics involve quadratic, cubic, or higher-order terms (e.g., fluid drag proportional to velocity squared).
Known trigonometric structure — the signal contains products of sinusoids or periodic components at known harmonics (e.g., sin²(t) = (1 - cos(2t)) / 2).
Standard DMD error is high — when fitting standard DMD produces large reconstruction or prediction error, lifting into a richer basis often improves accuracy.
Partial state observations — delay embedding reconstructs the full state space from scalar or low-dimensional measurements (Takens' theorem).

Caution: Lifting increases the dimension of the snapshot matrix. Very high polynomial degrees or large numbers of harmonics can lead to overfitting and increased computation time. Start with low settings and increase only if reconstruction error justifies it.

Comparison: Standard vs. Extended DMD

The following example computes DMD with and without polynomial lifting on a sin² signal and compares the reconstruction error (RMSE).

use koopman_dmd::{DmdConfig, LiftingConfig};
use ndarray::Array2;

fn main() {
    let n = 200;
    let dt = 0.05;
    let t: Vec<f64> = (0..n).map(|i| i as f64 * dt).collect();
    let signal: Vec<f64> = t.iter().map(|&ti| ti.sin().powi(2)).collect();
    let data = Array2::from_shape_vec((1, n), signal.clone()).unwrap();

    // Standard DMD (no lifting)
    let config_std = DmdConfig {
        dt: Some(dt),
        ..DmdConfig::default()
    };
    let result_std = koopman_dmd::fit(&data, &config_std).unwrap();
    let pred_std = result_std.reconstruct();

    // Extended DMD (polynomial degree 2)
    let config_ext = DmdConfig {
        lifting: Some(LiftingConfig::Polynomial { degree: 2 }),
        dt: Some(dt),
        ..DmdConfig::default()
    };
    let result_ext = koopman_dmd::fit(&data, &config_ext).unwrap();
    let pred_ext = result_ext.reconstruct();

    // Compute RMSE for each
    let rmse_std = rmse(&data, &pred_std);
    let rmse_ext = rmse(&data, &pred_ext);

    println!("Standard DMD RMSE: {:.6}", rmse_std);
    println!("Extended DMD RMSE: {:.6}", rmse_ext);
    // Extended DMD RMSE should be significantly lower
}

fn rmse(a: &Array2<f64>, b: &Array2<f64>) -> f64 {
    let diff = a - b;
    (diff.mapv(|v| v * v).mean().unwrap()).sqrt()
}

import numpy as np
from koopman_dmd import DmdConfig, LiftingConfig, fit

n = 200
dt = 0.05
t = np.arange(n) * dt
signal = np.sin(t) ** 2
data = signal.reshape(1, -1)

# Standard DMD
config_std = DmdConfig(dt=dt)
result_std = fit(data, config_std)
pred_std = result_std.reconstruct()

# Extended DMD with polynomial lifting
config_ext = DmdConfig(
    lifting=LiftingConfig.polynomial(degree=2),
    dt=dt,
)
result_ext = fit(data, config_ext)
pred_ext = result_ext.reconstruct()

# Compare RMSE
rmse_std = np.sqrt(np.mean((data - pred_std) ** 2))
rmse_ext = np.sqrt(np.mean((data - pred_ext) ** 2))

print(f"Standard DMD RMSE: {rmse_std:.6f}")
print(f"Extended DMD RMSE: {rmse_ext:.6f}")
# Extended DMD RMSE should be significantly lower

library(koopmandmd)

n <- 200
dt <- 0.05
t <- seq(0, by = dt, length.out = n)
signal <- sin(t)^2
data <- matrix(signal, nrow = 1)

# Standard DMD
config_std <- dmd_config(dt = dt)
result_std <- dmd_fit(data, config_std)
pred_std <- dmd_reconstruct(result_std)

# Extended DMD with polynomial lifting
config_ext <- dmd_config(
  lifting = lifting_polynomial(degree = 2),
  dt = dt
)
result_ext <- dmd_fit(data, config_ext)
pred_ext <- dmd_reconstruct(result_ext)

# Compare RMSE
rmse_std <- sqrt(mean((data - pred_std)^2))
rmse_ext <- sqrt(mean((data - pred_ext)^2))

cat(sprintf("Standard DMD RMSE: %.6f\n", rmse_std))
cat(sprintf("Extended DMD RMSE: %.6f\n", rmse_ext))
# Extended DMD RMSE should be significantly lower

Back-Projection

When Extended DMD is used, the internal computation operates on the lifted state vector, which has higher dimension than the original observations. However, all user-facing outputs — reconstructions, predictions, and mode shapes — are automatically projected back to the original observation space. This is handled through the LiftingInfo struct stored in the result.

LiftingInfo records:

original_dim — the number of rows in the original (un-lifted) snapshot matrix.
lifted_dim — the total number of rows after lifting.
config — the LiftingConfig variant that was applied, enabling the inverse mapping.

Because the original state variables always occupy the first original_dim rows of the lifted vector, back-projection is a simple extraction of those rows. This means:

result.predict(n) returns a matrix with the same number of rows as the input data, not the lifted dimension.
result.reconstruct() likewise returns data in the original space.
If you need the full lifted reconstruction for analysis, use result.reconstruct_lifted() to obtain the complete lifted state trajectory.

// Access lifting metadata
if let Some(ref info) = result.lifting_info {
    println!("Original dim: {}", info.original_dim);
    println!("Lifted dim:   {}", info.lifted_dim);
}

// Predictions are in original space
let pred = result.predict(100);
assert_eq!(pred.nrows(), 1); // original dimension

// Full lifted reconstruction if needed
let lifted = result.reconstruct_lifted();
assert_eq!(lifted.nrows(), 3); // [x, x^2] for degree-2 poly on 1-d input