Core DMD

Standard Dynamic Mode Decomposition for linear dynamical systems

Mathematical Background

Dynamic Mode Decomposition (DMD) extracts spatiotemporal coherent structures from time-series data by approximating the best-fit linear operator that advances the state forward in time. Given a sequence of snapshots arranged as columns of a data matrix X, the algorithm proceeds as follows.

Step 1 -- Split the data

Partition X into two overlapping matrices: X' containing the first n-1 columns and X'' containing the last n-1 columns. The relationship we seek is:

Step 2 -- Compute the SVD

Take the (optionally truncated) singular value decomposition of X':

where U is the left singular matrix, Σ is the diagonal matrix of singular values, and V* is the conjugate transpose of the right singular matrix.

Step 3 -- Project onto the POD basis

Rather than constructing the full high-dimensional operator A, project it onto the reduced basis spanned by the columns of U:

This small r x r matrix Ã captures the essential dynamics, where r is the truncation rank.

Step 4 -- Eigendecomposition

Compute the eigendecomposition of the projected matrix:

The DMD eigenvalues λi lie on the diagonal of Λ. The high-dimensional DMD modes are recovered as:

Step 5 -- Compute amplitudes

Solve for the amplitude vector b that best matches the initial snapshot x0:

The amplitudes weight the contribution of each mode to the overall dynamics.

DmdConfig Options

The DmdConfig struct controls how the decomposition is performed. The following fields are available:

use koopman_dmd::{DmdConfig};

// Auto-rank, no centering, no lifting
let config = DmdConfig::default();

// Explicit rank with centering
let config = DmdConfig {
    rank: Some(4),
    center: true,
    lifting: None,
};

from koopman_dmd import DmdConfig

# Auto-rank, no centering, no lifting
config = DmdConfig()

# Explicit rank with centering
config = DmdConfig(rank=4, center=True)

library(koopmandmd)

# Auto-rank, no centering, no lifting
config <- dmd_config()

# Explicit rank with centering
config <- dmd_config(rank = 4, center = TRUE)

Basic Usage

The core workflow is: build a configuration, pass your data matrix, and call the fit function. The result contains eigenvalues, modes, amplitudes, and metadata about the decomposition.

use koopman_dmd::{DmdConfig, dmd_fit};
use nalgebra as na;

// Build a 3x100 data matrix (3 state variables, 100 snapshots)
let data: Mat<f64> = generate_signal(3, 100);

let config = DmdConfig::default();
let result = dmd_fit(&data, &config);

// Inspect the result
println!("Rank:        {}", result.rank);
println!("Data dim:    {}", result.data_dim);
println!("Eigenvalues: {:?}", result.eigenvalues);
println!("Modes shape: {:?}", result.modes.shape());
println!("Amplitudes:  {:?}", result.amplitudes);

import numpy as np
from koopman_dmd import DmdConfig, dmd_fit

# Build a 3x100 data matrix (3 state variables, 100 snapshots)
data = generate_signal(3, 100)

config = DmdConfig()
result = dmd_fit(data, config)

# Inspect the result
print("Rank:       ", result.rank)
print("Data dim:   ", result.data_dim)
print("Eigenvalues:", result.eigenvalues)
print("Modes shape:", result.modes.shape)
print("Amplitudes: ", result.amplitudes)

library(koopmandmd)

# Build a 3x100 data matrix (3 state variables, 100 snapshots)
data <- generate_signal(3, 100)

config <- dmd_config()
result <- dmd_fit(data, config)

# Inspect the result
cat("Rank:       ", result$rank, "\n")
cat("Data dim:   ", result$data_dim, "\n")
print(result$eigenvalues)
cat("Modes dim:  ", dim(result$modes), "\n")
print(result$amplitudes)

Mean Centering

Setting center = true subtracts the temporal mean of each row before the decomposition and stores it for later reconstruction. This is useful when your data has a non-zero equilibrium or steady state.

Effects on eigenvalues:

• Without centering, you may see a dominant eigenvalue near 1 + 0i that corresponds to the steady state.
• With centering, that mode is absorbed into the stored mean and the remaining eigenvalues capture genuine oscillatory or decaying dynamics.
• Continuous-time eigenvalues (obtained via logarithm) shift accordingly: the near-zero continuous eigenvalue disappears when centering is enabled.

Rank Selection

Automatic rank

When rank is None, the library selects rank automatically by examining the singular values of X'. Singular values below a threshold (typically relative to the largest singular value) are discarded. This is a good default for exploratory analysis.

Manual rank

Setting an explicit rank lets you control the trade-off between reconstruction accuracy and model complexity. A lower rank produces a smoother, more interpretable model but may miss fine-grained dynamics. A higher rank captures more detail but can overfit to noise.

Guidelines for choosing rank:

• Plot the singular values in descending order. Look for an "elbow" where the values drop sharply -- this is often a good truncation point.
• Start with auto-rank and examine the reconstruction error. If the error is acceptable, the automatic choice is sufficient.
• For systems with known physics, set the rank to the expected number of dynamic modes (e.g., two for a simple harmonic oscillator).
• Compare relative reconstruction error across several rank values to find the point of diminishing returns.

Prediction

The library provides two prediction strategies. Both extrapolate the fitted model forward (or backward) in time, but they differ in how they compute subsequent states.

predict_modes

Reconstruct each time step directly from the DMD eigenvalues, modes, and amplitudes:

This method is non-recursive: every time step is computed independently from the spectral decomposition. It does not accumulate numerical error over long horizons, making it well-suited for long-range forecasts.

predict_matrix

Advance the state one step at a time by multiplying by the reconstructed operator:

This is a recursive approach: each step depends on the previous one. It can be faster for short horizons and naturally respects the linear map structure, but numerical errors may accumulate over many steps.

When to use each

• predict_modes -- preferred for long-range extrapolation, frequency analysis, or when you need arbitrary time indices (including non-integer steps in continuous time).
• predict_matrix -- preferred for short-range forecasting, real-time streaming applications, or when the full operator A is needed for control design.

use koopman_dmd::{dmd_fit, predict_modes, predict_matrix, DmdConfig};

let config = DmdConfig::default();
let result = dmd_fit(&data, &config);

// Predict 50 future time steps via spectral formula
let forecast_modes = predict_modes(&result, 50);

// Predict 50 future time steps via recursive matrix multiply
let forecast_matrix = predict_matrix(&result, 50);

from koopman_dmd import dmd_fit, predict_modes, predict_matrix, DmdConfig

config = DmdConfig()
result = dmd_fit(data, config)

# Predict 50 future time steps via spectral formula
forecast_modes = predict_modes(result, 50)

# Predict 50 future time steps via recursive matrix multiply
forecast_matrix = predict_matrix(result, 50)

library(koopmandmd)

config <- dmd_config()
result <- dmd_fit(data, config)

# Predict 50 future time steps via spectral formula
forecast_modes <- predict_modes(result, 50)

# Predict 50 future time steps via recursive matrix multiply
forecast_matrix <- predict_matrix(result, 50)

Reconstruction and Error

After fitting, you typically want to assess how well the DMD model reproduces the original data. The library provides two utilities for this purpose.

dmd_reconstruct

Produces the fitted values by reconstructing each snapshot from the DMD modes, eigenvalues, and amplitudes. If centering was used, the stored mean is added back automatically. The result has the same dimensions as the original data matrix.

dmd_error

Computes error metrics comparing the reconstruction to the original data. Returns:

• RMSE -- root mean squared error across all entries of the data matrix.
• Relative error -- the Frobenius norm of the residual divided by the Frobenius norm of the original data, expressed as a fraction.

use koopman_dmd::{dmd_fit, dmd_reconstruct, dmd_error, DmdConfig};

let config = DmdConfig { rank: Some(4), ..DmdConfig::default() };
let result = dmd_fit(&data, &config);

// Reconstruct the fitted values
let fitted = dmd_reconstruct(&result);

// Compute error metrics
let err = dmd_error(&data, &result);
println!("RMSE:           {:.6}", err.rmse);
println!("Relative error: {:.6}", err.relative);

from koopman_dmd import dmd_fit, dmd_reconstruct, dmd_error, DmdConfig

config = DmdConfig(rank=4)
result = dmd_fit(data, config)

# Reconstruct the fitted values
fitted = dmd_reconstruct(result)

# Compute error metrics
err = dmd_error(data, result)
print(f"RMSE:           {err.rmse:.6f}")
print(f"Relative error: {err.relative:.6f}")

library(koopmandmd)

config <- dmd_config(rank = 4)
result <- dmd_fit(data, config)

# Reconstruct the fitted values
fitted <- dmd_reconstruct(result)

# Compute error metrics
err <- dmd_error(data, result)
cat("RMSE:          ", err$rmse, "\n")
cat("Relative error:", err$relative, "\n")