Generalized Laplace Analysis

Extract Koopman eigenfunctions directly from trajectory data using weighted time averages.

Overview

Generalized Laplace Analysis (GLA) computes Koopman eigenfunctions directly via weighted time averages, without forming the DMD matrix. Based on the theoretical framework introduced by Mezic (2020), GLA provides a fundamentally different approach to spectral analysis of dynamical systems. Rather than constructing a finite-dimensional approximation of the Koopman operator, GLA iteratively refines eigenvalue estimates and extracts the corresponding eigenfunctions from trajectory data using spectral methods.

This approach is particularly valuable when specific eigenvalues are of interest or when forming the full DMD matrix is computationally expensive. GLA operates directly on observable time series, making it well-suited for streaming data and high-dimensional systems.

Mathematical Formulation

The Koopman operator U acts on observables of a dynamical system. A Koopman eigenfunction phi_lambda associated with eigenvalue lambda satisfies the functional equation:

where T is the dynamical system map. In other words, the eigenfunction transforms multiplicatively under the dynamics. GLA computes this eigenfunction via a weighted time average along trajectories:

Here, g is an observable function evaluated along the orbit of x. The weighting by lambda^{-k} isolates the component of the observable that evolves according to the eigenvalue lambda. When the eigenvalue estimate is correct, this sum converges to the projection of g onto the eigenfunction. The iterative refinement procedure adjusts the eigenvalue estimate to maximize convergence.

GlaConfig

The GlaConfig struct controls GLA computation. All parameters have sensible defaults, and the configuration can be customized as needed.

Basic Usage

The following example creates an oscillating two-variable signal, configures GLA to extract two eigenvalues, and inspects the results. For oscillatory data, the eigenvalues should have magnitude near 1.

use koopman_dmd::{GlaConfig, Matrix, gla};
use std::f64::consts::PI;

fn main() {
    // Create an oscillating 2-variable signal
    let n = 500;
    let dt = 0.02;
    let omega = 2.0 * PI * 0.5;
    let mut data = Matrix::zeros(2, n);
    for i in 0..n {
        let t = i as f64 * dt;
        data[(0, i)] = (omega * t).cos();
        data[(1, i)] = (omega * t).sin();
    }

    // Configure GLA with 2 eigenvalues
    let config = GlaConfig {
        eigenvalues: None,
        n_eigenvalues: 2,
        tol: 1e-6,
        max_iter: None,
    };

    // Run GLA
    let result = gla(&data, &config).unwrap();

    // Inspect eigenvalues -- magnitude should be near 1
    for (i, ev) in result.eigenvalues().iter().enumerate() {
        println!("eigenvalue {}: {} (|lambda| = {:.6})",
            i, ev, ev.norm());
    }

    // Check convergence info
    let info = result.convergence_info();
    println!("Converged: {}", info.converged);
    println!("Iterations: {}", info.iterations);
}

import koopman_dmd as kdmd
import numpy as np

# Create an oscillating 2-variable signal
n = 500
dt = 0.02
omega = 2.0 * np.pi * 0.5
t = np.arange(n) * dt
data = np.vstack([np.cos(omega * t), np.sin(omega * t)])

# Configure GLA with 2 eigenvalues
config = kdmd.GlaConfig(
    eigenvalues=None,
    n_eigenvalues=2,
    tol=1e-6,
    max_iter=None,
)

# Run GLA
result = kdmd.gla(data, config)

# Inspect eigenvalues -- magnitude should be near 1
for i, ev in enumerate(result.eigenvalues()):
    print(f"eigenvalue {i}: {ev} (|lambda| = {abs(ev):.6f})")

# Check convergence info
info = result.convergence_info()
print(f"Converged: {info.converged}")
print(f"Iterations: {info.iterations}")

library(koopman.dmd)

# Create an oscillating 2-variable signal
n <- 500
dt <- 0.02
omega <- 2 * pi * 0.5
t <- (0:(n - 1)) * dt
data <- rbind(cos(omega * t), sin(omega * t))

# Configure GLA with 2 eigenvalues
config <- gla_config(
  eigenvalues = NULL,
  n_eigenvalues = 2,
  tol = 1e-6,
  max_iter = NULL
)

# Run GLA
result <- gla(data, config)

# Inspect eigenvalues -- magnitude should be near 1
evs <- result$eigenvalues()
for (i in seq_along(evs)) {
  cat(sprintf("eigenvalue %d: %s (|lambda| = %.6f)\n",
      i, evs[i], Mod(evs[i])))
}

# Check convergence info
info <- result$convergence_info()
cat("Converged:", info$converged, "\n")
cat("Iterations:", info$iterations, "\n")

Prediction and Reconstruction

Once GLA has been fitted, you can use the extracted eigenfunctions and eigenvalues for forecasting with gla_predict and for computing fitted values with gla_reconstruct. These results can be compared against standard DMD to evaluate the quality of the GLA decomposition.

// Forecast 50 steps ahead using GLA
let forecast = gla_predict(&result, 50).unwrap();
println!("Forecast shape: {} x {}",
    forecast.nrows(), forecast.ncols());

// Reconstruct fitted values over the training window
let reconstruction = gla_reconstruct(&result).unwrap();

// Compare with standard DMD
let dmd_result = DMD::new(2).fit(&x, &y).unwrap();
let dmd_recon = dmd_result.reconstruct();

// Compute reconstruction error for both methods
let gla_error = (&reconstruction - &data).norm() / data.norm();
let dmd_error = (&dmd_recon - &data).norm() / data.norm();
println!("GLA reconstruction error: {:.6e}", gla_error);
println!("DMD reconstruction error: {:.6e}", dmd_error);

# Forecast 50 steps ahead using GLA
forecast = kdmd.gla_predict(result, n_steps=50)
print(f"Forecast shape: {forecast.shape}")

# Reconstruct fitted values over the training window
reconstruction = kdmd.gla_reconstruct(result)

# Compare with standard DMD
dmd = kdmd.DMD(rank=2)
dmd_result = dmd.fit(data[:, :-1], data[:, 1:])
dmd_recon = dmd_result.reconstruct()

# Compute reconstruction error for both methods
gla_error = np.linalg.norm(reconstruction - data) / np.linalg.norm(data)
dmd_error = np.linalg.norm(dmd_recon - data) / np.linalg.norm(data)
print(f"GLA reconstruction error: {gla_error:.6e}")
print(f"DMD reconstruction error: {dmd_error:.6e}")

# Forecast 50 steps ahead using GLA
forecast <- gla_predict(result, n_steps = 50)
cat("Forecast dimensions:", dim(forecast), "\n")

# Reconstruct fitted values over the training window
reconstruction <- gla_reconstruct(result)

# Compare with standard DMD
x <- data[, 1:(n - 1)]
y <- data[, 2:n]
dmd_result <- dmd_fit(x, y, rank = 2)
dmd_recon <- dmd_result$reconstruct()

# Compute reconstruction error for both methods
gla_error <- norm(reconstruction - data, "F") / norm(data, "F")
dmd_error <- norm(dmd_recon - data, "F") / norm(data, "F")
cat(sprintf("GLA reconstruction error: %.6e\n", gla_error))
cat(sprintf("DMD reconstruction error: %.6e\n", dmd_error))

When to Use GLA vs DMD

GLA and standard DMD approach the same underlying problem from different angles. The right choice depends on your specific use case and computational constraints.

GLA is better for

• Extracting specific eigenvalues -- when you know which frequencies or growth rates are important and want to target them directly, GLA can isolate individual Koopman eigenvalues without computing the full spectrum.
• Large systems where the full DMD matrix is expensive -- GLA avoids constructing and decomposing the DMD matrix, operating instead on time averages along trajectories. This makes it more scalable for very high-dimensional data.
• Direct eigenfunction access -- GLA naturally produces eigenfunctions evaluated along the trajectory, providing a direct connection to Koopman operator theory.

DMD is better for

• General-purpose decomposition -- when you want a complete picture of the system dynamics without prior knowledge of which eigenvalues matter, standard DMD computes the full set of modes simultaneously.
• Modes and amplitudes -- DMD provides spatial modes and their amplitudes in a single computation, making it straightforward to rank modes by energy and identify dominant structures.
• Fast computation on moderate-sized problems -- for problems where the snapshot matrices fit comfortably in memory, DMD's SVD-based approach is typically faster than GLA's iterative refinement.

Convergence

GLA uses an iterative refinement procedure to improve eigenvalue estimates. At each iteration, the weighted time average is computed with the current eigenvalue guess, and the eigenvalue is updated based on the resulting eigenfunction estimate. The process continues until the change between successive eigenvalue estimates falls below the specified tolerance tol, or until the maximum number of iterations max_iter is reached.

The tol parameter controls the precision of the final eigenvalue estimates. A tolerance of 1e-6 (the default) is suitable for most applications. For high-precision work, values of 1e-10 or smaller can be used, though this will increase computation time. The convergence rate depends on the spectral gap between eigenvalues: well-separated eigenvalues converge quickly, while closely spaced eigenvalues may require more iterations.