Analysis Tools

Post-decomposition diagnostics for understanding system dynamics, stability, and reconstruction quality.

Overview

After computing a DMD decomposition, the analysis module provides a suite of diagnostic tools for interpreting and validating the results. These tools cover:

Spectrum analysis -- frequency, growth rate, and damping ratio for each mode.
Stability analysis -- classification of eigenvalues relative to the unit circle.
Error metrics -- RMSE, relative error, and maximum pointwise error against original data.
Residuals -- per-mode residual norms measuring how well each mode satisfies the DMD eigenvalue equation.
Dominant modes -- extraction of the most significant modes by amplitude or energy.
Pseudospectrum -- resolvent norm computation for sensitivity analysis of non-normal operators.
Convergence -- evaluation of reconstruction error as a function of rank truncation.

Spectrum Analysis

The dmd_spectrum function extracts spectral information from a completed DMD result. Given the time step dt, it converts discrete-time eigenvalues into continuous-time quantities and classifies the stability of each mode.

For each mode, the returned ModeInfo contains:

Field	Description
`frequency`	Oscillation frequency in Hz, derived from the argument of the eigenvalue.
`magnitude`	Absolute value of the eigenvalue, \|lambda\|.
`growth_rate`	Continuous-time growth rate, computed as log(\|lambda\|) / dt.
`damping_ratio`	Ratio quantifying the rate of amplitude decay relative to the oscillation frequency.
`stability`	Classification: `Growing` (\|lambda\| > 1), `Decaying` (\|lambda\| < 1), or `Neutral` (\|lambda\| = 1 within tolerance).

use koopman_dmd::analysis::{dmd_spectrum, Stability};

let dt = 0.01;
let spectrum = dmd_spectrum(&result, dt);

for mode in &spectrum {
    println!(
        "freq={:.3} Hz  |lambda|={:.4}  growth={:.4}  {:?}",
        mode.frequency,
        mode.magnitude,
        mode.growth_rate,
        mode.stability,
    );
}

// Filter for growing modes only
let growing: Vec<_> = spectrum
    .iter()
    .filter(|m| matches!(m.stability, Stability::Growing))
    .collect();

import koopman_dmd as kdmd

dt = 0.01
spectrum = kdmd.dmd_spectrum(result, dt)

for mode in spectrum:
    print(
        f"freq={mode.frequency:.3f} Hz  "
        f"|lambda|={mode.magnitude:.4f}  "
        f"growth={mode.growth_rate:.4f}  "
        f"{mode.stability}"
    )

# Filter for growing modes only
growing = [m for m in spectrum if m.stability == "Growing"]

library(koopmandmd)

dt <- 0.01
spectrum <- dmd_spectrum(result, dt)

# spectrum is a data.frame with columns:
# frequency, magnitude, growth_rate, damping_ratio, stability
print(spectrum)

# Filter for growing modes only
growing <- spectrum[spectrum$stability == "Growing", ]

Stability Analysis

The dmd_stability function provides a summary assessment of the system's stability based on eigenvalue locations relative to the unit circle. An optional tolerance parameter controls the width of the neutral band around |lambda| = 1.

The returned StabilityResult contains:

Field	Description
`spectral_radius`	Maximum eigenvalue magnitude, max(\|lambda\|).
`is_stable`	True if all eigenvalues lie within the unit circle (within tolerance).
`is_unstable`	True if any eigenvalue lies strictly outside the unit circle.
`n_growing`	Number of modes with \|lambda\| > 1 + tol.
`n_decaying`	Number of modes with \|lambda\| < 1 - tol.
`n_neutral`	Number of modes with \|lambda\| within the tolerance band around 1.

use koopman_dmd::analysis::dmd_stability;

let tol = 1e-6;
let stab = dmd_stability(&result, tol);

println!("Spectral radius: {:.6}", stab.spectral_radius);
println!("Stable: {}", stab.is_stable);
println!(
    "Modes -- growing: {}, decaying: {}, neutral: {}",
    stab.n_growing, stab.n_decaying, stab.n_neutral
);

import koopman_dmd as kdmd

tol = 1e-6
stab = kdmd.dmd_stability(result, tol)

print(f"Spectral radius: {stab.spectral_radius:.6f}")
print(f"Stable: {stab.is_stable}")
print(
    f"Modes -- growing: {stab.n_growing}, "
    f"decaying: {stab.n_decaying}, "
    f"neutral: {stab.n_neutral}"
)

library(koopmandmd)

tol <- 1e-6
stab <- dmd_stability(result, tol)

cat("Spectral radius:", stab$spectral_radius, "\n")
cat("Stable:", stab$is_stable, "\n")
cat("Modes -- growing:", stab$n_growing,
    "decaying:", stab$n_decaying,
    "neutral:", stab$n_neutral, "\n")

Error Metrics

The dmd_error function compares the DMD reconstruction against the original data matrix and returns an ErrorMetrics struct summarizing the approximation quality.

Field	Description
`rmse`	Root mean square error between reconstruction and original data.
`relative_error`	RMSE normalized by the Frobenius norm of the original data.
`max_error`	Maximum absolute pointwise error across all entries.

use koopman_dmd::analysis::dmd_error;

let err = dmd_error(&result, &original_data);

println!("RMSE:           {:.6e}", err.rmse);
println!("Relative error: {:.6e}", err.relative_error);
println!("Max error:      {:.6e}", err.max_error);

import koopman_dmd as kdmd

err = kdmd.dmd_error(result, original_data)

print(f"RMSE:           {err.rmse:.6e}")
print(f"Relative error: {err.relative_error:.6e}")
print(f"Max error:      {err.max_error:.6e}")

library(koopmandmd)

err <- dmd_error(result, original_data)

cat(sprintf("RMSE:           %.6e\n", err$rmse))
cat(sprintf("Relative error: %.6e\n", err$relative_error))
cat(sprintf("Max error:      %.6e\n", err$max_error))

Residuals

The dmd_residual function computes per-mode residual norms. For each mode, the residual measures how well the mode-eigenvalue pair satisfies the fundamental DMD equation:

r_i = || A * phi_i - lambda_i * phi_i ||

where A is the approximated linear operator, phi_i is the i-th DMD mode, and lambda_i is the corresponding eigenvalue. Small residuals indicate that the mode is a faithful eigenvector of the underlying dynamics. Large residuals may signal numerical issues, insufficient rank, or strongly nonlinear behavior that the linear operator cannot capture.

use koopman_dmd::analysis::dmd_residual;

let residuals = dmd_residual(&result);

for (i, r) in residuals.iter().enumerate() {
    println!("Mode {}: residual = {:.6e}", i, r);
}

import koopman_dmd as kdmd

residuals = kdmd.dmd_residual(result)

for i, r in enumerate(residuals):
    print(f"Mode {i}: residual = {r:.6e}")

library(koopmandmd)

residuals <- dmd_residual(result)

for (i in seq_along(residuals)) {
    cat(sprintf("Mode %d: residual = %.6e\n", i, residuals[i]))
}

Dominant Modes

The dmd_dominant_modes function extracts the top n most significant modes according to a specified criterion. This is useful for reducing a high-rank decomposition to its most physically meaningful components.

Two selection criteria are available:

Criterion	Formula	Interpretation
`DominantCriterion::Amplitude`	\|b_i\|	Modes with the largest initial amplitude contribution. Ranks by how strongly each mode participates in the initial condition.
`DominantCriterion::Energy`	\|b_i\|^2 * \|lambda_i\|	Modes carrying the most energy, accounting for both amplitude and persistence. Modes that are large but rapidly decaying are penalized.

use koopman_dmd::analysis::{dmd_dominant_modes, DominantCriterion};

// Top 5 modes by amplitude
let top_amp = dmd_dominant_modes(
    &result,
    DominantCriterion::Amplitude,
    5,
);

// Top 5 modes by energy
let top_energy = dmd_dominant_modes(
    &result,
    DominantCriterion::Energy,
    5,
);

for (idx, score) in &top_energy {
    println!("Mode {}: energy score = {:.6}", idx, score);
}

import koopman_dmd as kdmd

# Top 5 modes by amplitude
top_amp = kdmd.dmd_dominant_modes(result, criterion="amplitude", n=5)

# Top 5 modes by energy
top_energy = kdmd.dmd_dominant_modes(result, criterion="energy", n=5)

for idx, score in top_energy:
    print(f"Mode {idx}: energy score = {score:.6f}")

library(koopmandmd)

# Top 5 modes by amplitude
top_amp <- dmd_dominant_modes(result, criterion = "amplitude", n = 5)

# Top 5 modes by energy
top_energy <- dmd_dominant_modes(result, criterion = "energy", n = 5)

# top_energy is a data.frame with columns: index, score
print(top_energy)

Pseudospectrum

The dmd_pseudospectrum function computes the resolvent norm over a grid in the complex plane. The pseudospectrum reveals the sensitivity of eigenvalues to perturbations, which is especially important for non-normal operators where eigenvalues alone can be misleading.

For a given complex point z, the resolvent norm is defined as:

sigma(z) = || (zI - A)^{-1} ||

The function accepts vectors defining the real and imaginary axes of the grid and returns a 2D array of resolvent norm values. Regions of large resolvent norm indicate where small perturbations to the operator could shift eigenvalues, highlighting potential transient growth even in nominally stable systems.

use koopman_dmd::analysis::dmd_pseudospectrum;
use ndarray::Array1;

// Define a grid over the complex plane
let grid_re = Array1::linspace(-1.5, 1.5, 200);
let grid_im = Array1::linspace(-1.5, 1.5, 200);

// Compute resolvent norm on the grid
let sigma = dmd_pseudospectrum(&result, &grid_re, &grid_im);

// sigma is a 200x200 Array2<f64> of resolvent norms
println!("Pseudospectrum shape: {:?}", sigma.dim());

import numpy as np
import koopman_dmd as kdmd

# Define a grid over the complex plane
grid_re = np.linspace(-1.5, 1.5, 200)
grid_im = np.linspace(-1.5, 1.5, 200)

# Compute resolvent norm on the grid
sigma = kdmd.dmd_pseudospectrum(result, grid_re, grid_im)

# sigma is a (200, 200) ndarray of resolvent norms
print(f"Pseudospectrum shape: {sigma.shape}")

library(koopmandmd)

# Define a grid over the complex plane
grid_re <- seq(-1.5, 1.5, length.out = 200)
grid_im <- seq(-1.5, 1.5, length.out = 200)

# Compute resolvent norm on the grid
sigma <- dmd_pseudospectrum(result, grid_re, grid_im)

# sigma is a 200x200 matrix of resolvent norms
cat("Pseudospectrum dimensions:", dim(sigma), "\n")

Convergence

The dmd_convergence function evaluates how the DMD reconstruction error changes as a function of rank truncation. This is useful for selecting the optimal number of modes: too few modes underfit the data, while too many modes risk capturing noise.

Given a range of ranks, the function performs a DMD fit at each rank and returns a vector of (rank, error) pairs. The error is the relative reconstruction error at each truncation level. A characteristic "elbow" in the resulting curve typically indicates the point where additional modes contribute primarily noise rather than signal.

use koopman_dmd::analysis::dmd_convergence;

// Evaluate DMD error for ranks 1 through 20
let rank_range = 1..=20;
let convergence = dmd_convergence(&data, rank_range);

for (rank, error) in &convergence {
    println!("rank={:2}  relative_error={:.6e}", rank, error);
}

// Find the rank with the steepest drop in error
let best_rank = convergence
    .windows(2)
    .enumerate()
    .max_by(|(_, a), (_, b)| {
        let drop_a = a[0].1 - a[1].1;
        let drop_b = b[0].1 - b[1].1;
        drop_a.partial_cmp(&drop_b).unwrap()
    })
    .map(|(i, _)| i + 2);
println!("Suggested rank: {:?}", best_rank);

import koopman_dmd as kdmd
import numpy as np

# Evaluate DMD error for ranks 1 through 20
rank_range = range(1, 21)
convergence = kdmd.dmd_convergence(data, rank_range)

for rank, error in convergence:
    print(f"rank={rank:2d}  relative_error={error:.6e}")

# Find the elbow point
errors = np.array([e for _, e in convergence])
drops = np.diff(errors)
best_rank = np.argmin(drops) + 2
print(f"Suggested rank: {best_rank}")

library(koopmandmd)

# Evaluate DMD error for ranks 1 through 20
rank_range <- 1:20
convergence <- dmd_convergence(data, rank_range)

# convergence is a data.frame with columns: rank, error
print(convergence)

# Plot the convergence curve
plot(convergence$rank, convergence$error,
     type = "b", log = "y",
     xlab = "Rank", ylab = "Relative Error",
     main = "DMD Convergence")