R API Reference

Complete reference for the koopman.dmd R package -- Rust-powered Koopman DMD with idiomatic S3 classes.

Installation

The koopman.dmd package is an R wrapper around the Rust koopman-dmd library, available on CRAN:

install.packages("koopman.dmd")

CRAN provides prebuilt binaries for Windows and macOS. Installing from source compiles the bundled Rust dependency tree locally, which requires a working Rust toolchain (cargo and rustc >= 1.85) in addition to the standard R development tools. To install from a checkout of the repository:

# Install from the local source directory
R CMD INSTALL koopman-dmd-r

Verify the installation by loading the package in an R session:

library(koopman.dmd)
cat("koopman.dmd loaded successfully\n")

DMD Functions

dmd

Perform Dynamic Mode Decomposition on a numeric matrix. This is the primary entry point for standard and extended DMD analysis. Returns an S3 object of class "dmd".

Return value. An S3 object of class "dmd" with the following components:

library(koopman.dmd)

# Generate sample data: 3 variables, 100 time steps
t <- seq(0, 2 * pi, length.out = 100)
data <- rbind(
  sin(t),
  cos(t),
  sin(2 * t)
)

# Standard DMD with auto-rank
result <- dmd(data)

# DMD with explicit rank and centering
result <- dmd(data, rank = 4, center = TRUE)

# Extended DMD with polynomial lifting
result <- dmd(data, lifting = "polynomial", lifting_param = 3L)

# Extended DMD with trigonometric lifting
result <- dmd(data, lifting = "trigonometric", lifting_param = 5L)

# Extended DMD with delay embedding
result <- dmd(data, lifting = "delay", lifting_param = 10L)

print.dmd

S3 print method for "dmd" objects. Displays a concise summary including the rank, number of variables, and number of time steps.

print(result)
# DMD Result
#   Rank:       4
#   Variables:  3
#   Time steps: 100

summary.dmd

S3 summary method for "dmd" objects. Prints a table of eigenvalues with their real parts, imaginary parts, moduli, and frequencies.

summary(result)
# DMD Eigenvalues:
#   Re(lambda)  Im(lambda)  |lambda|  Frequency
#   0.9998      0.0175      1.0000    0.0028
#   0.9998     -0.0175      1.0000   -0.0028
#   ...

predict.dmd

S3 predict method for "dmd" objects. Extrapolates the fitted DMD model forward in time.

Return value. A numeric matrix of dimensions n_vars x n_ahead containing the predicted states.

# Predict 50 future time steps
forecast <- predict(result, n_ahead = 50)
cat("Forecast dimensions:", dim(forecast), "\n")

# Predict from a custom initial condition
x0 <- c(0.5, 0.5, 0.5)
forecast <- predict(result, n_ahead = 50, x0 = x0)

dmd_reconstruct

Reconstruct the original data from the DMD decomposition. Produces the fitted values by combining modes, eigenvalues, and amplitudes. If centering was used, the stored mean is added back automatically.

Return value. A numeric matrix of dimensions n_vars x n_time containing the reconstructed data.

# Reconstruct the original 100 time steps
recon <- dmd_reconstruct(result, n_time = 100)
cat("Reconstruction dimensions:", dim(recon), "\n")

dmd_spectrum

Compute the spectral properties of the DMD decomposition. Returns a data frame summarizing the frequency content, growth rates, damping, and stability of each mode.

Return value. A data.frame with the following columns:

# Compute spectrum with default time step
spec <- dmd_spectrum(result)
print(spec)

# Compute spectrum with explicit time step
spec <- dmd_spectrum(result, dt = 0.01)
print(spec)

dmd_stability

Assess the stability of the DMD model by examining the eigenvalue spectrum.

Return value. A list with the following components:

stab <- dmd_stability(result)
cat("Spectral radius:", stab$spectral_radius, "\n")
cat("Stable:", stab$is_stable, "\n")
cat("Growing modes:", stab$n_growing, "\n")
cat("Decaying modes:", stab$n_decaying, "\n")
cat("Neutral modes:", stab$n_neutral, "\n")

dmd_error

Compute reconstruction error metrics for a DMD result against the original data.

Return value. A list with the following components:

err <- dmd_error(result, data)
cat("RMSE:          ", err$rmse, "\n")
cat("Relative error:", err$relative_error, "\n")
cat("Max error:     ", err$max_error, "\n")

dmd_dominant_modes

Extract the dominant DMD modes ranked by a specified criterion.

Return value. A data.frame with one row per dominant mode, including eigenvalue, frequency, amplitude, and growth rate columns.

# Get top 3 modes by amplitude
top <- dmd_dominant_modes(result, criterion = "amplitude", n = 3)
print(top)

dmd_residual

Compute residuals for each DMD mode, measuring how well each mode satisfies the linear operator relationship.

Return value. A list with the following components:

resid <- dmd_residual(result)
cat("Max residual: ", resid$max_residual, "\n")
cat("Mean residual:", resid$mean_residual, "\n")
print(resid$mode_residuals)

DMDc Functions

dmdc

Perform Dynamic Mode Decomposition with control (DMDc; Proctor, Brunton & Kutz 2016, doi:10.1137/15M1013857). Identifies the forced linear system xt+1 = A xt + B ut from snapshot pairs. Unlike dmd(), which takes one contiguous trajectory, dmdc() takes explicit pair matrices: X1 holds states at time t, X2 the states one step later, and U the control input applied during each transition. Columns may come from many concatenated trajectories. Returns an S3 object of class "dmdc".

Return value. An S3 object of class "dmdc" with the following components:

library(koopman.dmd)

# Simulate a forced linear system x_{t+1} = A0 x_t + B0 u_t
A0 <- matrix(c(0.9, 0, 0.1, 0.8), 2, 2)
B0 <- matrix(c(0.5, 1), 2, 1)
m <- 120
X1 <- matrix(0, 2, m); X2 <- matrix(0, 2, m); U <- matrix(0, 1, m)
x <- c(1, -0.5)
for (t in seq_len(m)) {
  u_t <- sin(0.7 * (t - 1)) + 0.5 * cos(2.3 * (t - 1) + 1)
  X1[, t] <- x
  U[, t] <- u_t
  x <- as.numeric(A0 %*% x + B0 * u_t)
  X2[, t] <- x
}

# Identify A and B from the snapshot pairs
fit <- dmdc(X1, X2, U, rank_input = 3)
round(fit$a, 6)   # recovers A0
round(fit$b, 6)   # recovers B0

# Known input matrix: estimate A only
fit_kb <- dmdc(X1, X2, U, rank_input = 2, known_B = B0)

# Autonomous multi-trajectory model from the same pairs
fit_auto <- dmdc(X1, X2)

print.dmdc

S3 print method for "dmdc" objects. Displays a concise one-line summary of the system dimensions and ranks.

print(fit)
# DMDc(n_states=2, n_inputs=1, rank_input=3, rank_output=2)

summary.dmdc

S3 summary method for "dmdc" objects. Prints the numbers of states and inputs, the truncation ranks, the time step, the singular values of the regression-input SVD, and the eigenvalue magnitudes.

summary(fit)
# DMDc(n_states=2, n_inputs=1, rank_input=3, rank_output=2)
#   dt:                    1
#   Singular values:       ...
#   Eigenvalue magnitudes: 0.9  0.8

predict.dmdc

S3 predict method for "dmdc" objects. Simulates the identified system xt+1 = A xt + B ut forward under a supplied input sequence.

Return value. A numeric matrix of dimensions n_states x k containing the predicted states x₁..x_k -- the successors of x0; x0 itself is not included.

# Replay the training inputs -- reproduces X2
pred <- predict(fit, U = U)

# Unforced (zero-input) response from a custom initial state
free <- predict(fit, x0 = c(1, 0), n_ahead = 20)

dmdc_stability

Assess the stability of the identified operator from the eigenvalues of a_tilde -- the spectrum of the unforced dynamics. The result has the same shape as dmd_stability().

Return value. A list with the following components:

stab <- dmdc_stability(fit)
cat("Spectral radius:", stab$spectral_radius, "\n")
cat("Stable:", stab$is_stable, "\n")

dmdc_spectrum

Compute per-mode frequency, growth rate, and stability for the identified operator, in the same shape as dmd_spectrum(). DMDc has no mode amplitudes, so the amplitude column is always 0.

Return value. A data.frame with the following columns:

# Spectrum with the stored time step
spec <- dmdc_spectrum(fit)
print(spec)

# Override the time step
spec <- dmdc_spectrum(fit, dt = 0.1)

Hankel-DMD Functions

hankel_dmd

Perform Hankel-DMD, which augments the data with time-delay embeddings before applying DMD. This is particularly effective for scalar time series or systems where a single observable does not fully capture the underlying dynamics. Returns an S3 object of class "hankel_dmd".

Return value. An S3 object of class "hankel_dmd" with components similar to the standard "dmd" class, plus delay embedding metadata.

library(koopman.dmd)

# Scalar time series: 1 x 500 matrix
t <- seq(0, 10 * pi, length.out = 500)
data <- matrix(sin(t) + 0.5 * sin(3 * t), nrow = 1)

# Hankel-DMD with 20 delays
hresult <- hankel_dmd(data, delays = 20)
print(hresult)

print.hankel_dmd

S3 print method for "hankel_dmd" objects. Displays a concise summary including the rank, delay embedding depth, and data dimensions.

print(hresult)
# Hankel-DMD Result
#   Rank:    10
#   Delays:  20
#   Vars:    1
#   Steps:   500

predict.hankel_dmd

S3 predict method for "hankel_dmd" objects. Extrapolates the fitted Hankel-DMD model forward in time.

Return value. A numeric matrix containing the predicted states in the original (non-embedded) variable space.

# Predict 100 steps ahead
hforecast <- predict(hresult, n_ahead = 100)
cat("Forecast dimensions:", dim(hforecast), "\n")

hankel_reconstruct

Reconstruct the original data from the Hankel-DMD decomposition, projecting back from the delay-embedded space to the original variable space.

Return value. A numeric matrix containing the reconstructed data.

hrecon <- hankel_reconstruct(hresult)
cat("Reconstruction dimensions:", dim(hrecon), "\n")

GLA Functions

gla

Perform Generalized Laplace Analysis (GLA) to extract Koopman eigenvalues and eigenfunctions from trajectory data. GLA is an iterative algorithm that converges to the dominant Koopman modes. Returns an S3 object of class "gla".

Return value. An S3 object of class "gla" containing the extracted Koopman eigenvalues and associated data.

library(koopman.dmd)

# Generate trajectory data
t <- seq(0, 4 * pi, length.out = 200)
data <- rbind(sin(t), cos(t))

# Extract 4 Koopman eigenvalues
gresult <- gla(data, n_eigenvalues = 4, tol = 1e-8)
print(gresult)

print.gla

S3 print method for "gla" objects. Displays a summary of the GLA result including the number of extracted eigenvalues and convergence information.

print(gresult)
# GLA Result
#   Eigenvalues: 4
#   Converged:   TRUE

predict.gla

S3 predict method for "gla" objects. Extrapolates the fitted GLA model forward in time using the extracted Koopman modes.

Return value. A numeric matrix containing the predicted states.

# Predict 50 steps ahead
gforecast <- predict(gresult, n_ahead = 50)
cat("Forecast dimensions:", dim(gforecast), "\n")

gla_reconstruct

Reconstruct the original data from the GLA decomposition.

Return value. A numeric matrix containing the reconstructed data.

grecon <- gla_reconstruct(gresult)
cat("Reconstruction dimensions:", dim(grecon), "\n")

Map Functions

generate_trajectory

Iterate a dynamical map from a given initial condition and return a matrix of states. Each column of the returned matrix is a snapshot of the system state at a single time step.

Return value. A numeric matrix of dimensions n_dim x n_steps where n_dim is the dimensionality of the map.

library(koopman.dmd)

# Standard map
traj <- generate_trajectory(
  ic = c(0.1, 0.2),
  map = "standard",
  n_steps = 1000,
  epsilon = 0.12
)
cat("Trajectory shape:", dim(traj), "\n")  # 2 x 1000

# Henon map
traj <- generate_trajectory(
  ic = c(0.0, 0.0),
  map = "henon",
  n_steps = 5000,
  a = 1.4,
  b = 0.3
)

# Logistic map
traj <- generate_trajectory(
  ic = 0.4,
  map = "logistic",
  n_steps = 500,
  r = 3.9
)

Helper Constructors

Convenience functions that return named lists of map parameters with sensible defaults. These can be expanded with do.call or used for documentation purposes.

# Using helper constructors
params <- standard_map(epsilon = 0.25)
print(params)
# $map
# [1] "standard"
# $epsilon
# [1] 0.25

params <- henon_map()  # uses defaults a=1.4, b=0.3
params <- logistic_map(r = 3.7)

# Use with do.call
traj <- do.call(generate_trajectory, c(
  list(ic = c(0.1, 0.2), n_steps = 1000),
  standard_map(epsilon = 0.3)
))

Harmonic Analysis Functions

harmonic_time_average

Compute the harmonic time average (HTA) of an observable along a trajectory generated by a dynamical map. The HTA is a key tool in Koopman operator theory for extracting frequency-specific information from orbits.

Return value. A list with the following components:

library(koopman.dmd)

# HTA of identity observable on standard map
hta <- harmonic_time_average(
  ic = c(0.1, 0.2),
  map = "standard",
  observable = "identity",
  omega = 0.0,
  n_iter = 10000,
  epsilon = 0.12
)
cat("Magnitude:", hta$magnitude, "\n")
cat("Phase:    ", hta$phase, "\n")

mesochronic_compute

Compute harmonic time averages over a grid of initial conditions, producing a mesochronic map of the phase space. This reveals the structure of Koopman eigenfunctions and distinguishes regular from chaotic regions.

Return value. A list with the following components:

# Mesochronic analysis of the standard map
meso <- mesochronic_compute(
  map = "standard",
  x_range = c(0, 1),
  y_range = c(0, 1),
  resolution = 200,
  observable = "identity",
  omega = 0.0,
  n_iter = 5000,
  epsilon = 0.12
)

# Plot with base R
image(
  meso$x_coords,
  meso$y_coords,
  meso$hta_matrix,
  col = hcl.colors(256, "viridis"),
  xlab = "x",
  ylab = "y",
  main = "Mesochronic Map"
)

classify_phase_space

Classify an initial condition as belonging to a regular or chaotic region of phase space, based on the convergence behavior of the harmonic time average.

Return value. A list with the following components:

# Classify a point on the standard map
cls <- classify_phase_space(
  ic = c(0.1, 0.2),
  map = "standard",
  omega = 0.0,
  n_iter = 10000,
  epsilon = 0.12
)
cat("Classification:", cls$classification, "\n")
cat("HTA magnitude: ", cls$hta_magnitude, "\n")

hta_convergence

Track the convergence of the harmonic time average as a function of the number of iterations. Useful for determining the required iteration depth and for visualizing the convergence behavior of regular versus chaotic orbits.

Return value. A data.frame with two columns:

# Track HTA convergence for a regular orbit
conv <- hta_convergence(
  ic = c(0.1, 0.2),
  map = "standard",
  observable = "identity",
  omega = 0.0,
  max_iter = 50000,
  step = 500,
  epsilon = 0.12
)
plot(conv$iteration, conv$magnitude,
     type = "l", xlab = "Iteration", ylab = "HTA Magnitude")

Examples

Complete DMD Workflow

library(koopman.dmd)

# -- Generate synthetic data --
t <- seq(0, 4 * pi, length.out = 200)
data <- rbind(
  sin(t) + 0.1 * sin(5 * t),
  cos(t) + 0.1 * cos(5 * t)
)

# -- Fit DMD --
result <- dmd(data, rank = 4)
print(result)
summary(result)

# -- Prediction --
forecast <- predict(result, n_ahead = 100)
matplot(t(forecast), type = "l", main = "DMD Forecast")

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

# -- Spectrum --
spec <- dmd_spectrum(result, dt = 0.01)
print(spec)

# -- Stability --
stab <- dmd_stability(result)
cat("Stable:", stab$is_stable, "\n")

# -- Dominant modes --
top <- dmd_dominant_modes(result, n = 3)
print(top)

Hankel-DMD for Scalar Time Series

library(koopman.dmd)

# Scalar time series with two frequencies
t <- seq(0, 20 * pi, length.out = 1000)
data <- matrix(sin(t) + 0.5 * sin(3 * t), nrow = 1)

# Fit Hankel-DMD
hresult <- hankel_dmd(data, delays = 30, rank = 4)
print(hresult)

# Predict and reconstruct
hforecast <- predict(hresult, n_ahead = 200)
hrecon <- hankel_reconstruct(hresult)

# Compare reconstruction to original
plot(data[1, ], type = "l", col = "black", ylab = "x",
     main = "Hankel-DMD Reconstruction")
lines(hrecon[1, ], col = "red", lty = 2)
legend("topright", c("Original", "Reconstructed"),
       col = c("black", "red"), lty = c(1, 2))

GLA Eigenvalue Extraction

library(koopman.dmd)

# Oscillatory data
t <- seq(0, 8 * pi, length.out = 500)
data <- rbind(sin(t), cos(t))

# Extract 4 Koopman eigenvalues with GLA
gresult <- gla(data, n_eigenvalues = 4, tol = 1e-8)
print(gresult)

# Predict and reconstruct
gforecast <- predict(gresult, n_ahead = 100)
grecon <- gla_reconstruct(gresult)

matplot(t(grecon), type = "l", main = "GLA Reconstruction")

Trajectory Generation and Phase Portraits

library(koopman.dmd)

# Standard map phase portrait
par(mfrow = c(1, 2))

# Regular orbit
traj1 <- generate_trajectory(
  ic = c(0.1, 0.2),
  map = "standard",
  n_steps = 5000,
  epsilon = 0.12
)
plot(traj1[1, ], traj1[2, ],
     pch = ".", xlab = "x", ylab = "y",
     main = "Standard Map (regular)")

# Chaotic orbit
traj2 <- generate_trajectory(
  ic = c(0.5, 0.5),
  map = "standard",
  n_steps = 50000,
  epsilon = 0.97
)
plot(traj2[1, ], traj2[2, ],
     pch = ".", xlab = "x", ylab = "y",
     main = "Standard Map (chaotic)")

# Henon attractor
par(mfrow = c(1, 1))
henon <- generate_trajectory(
  ic = c(0.0, 0.0),
  map = "henon",
  n_steps = 10000,
  a = 1.4,
  b = 0.3
)
plot(henon[1, ], henon[2, ],
     pch = ".", xlab = "x", ylab = "y",
     main = "Henon Attractor")

# Logistic map bifurcation diagram
rs <- seq(2.5, 4.0, by = 0.005)
plot(NULL, xlim = c(2.5, 4), ylim = c(0, 1),
     xlab = "r", ylab = "x", main = "Logistic Map Bifurcation")
for (r in rs) {
  traj <- generate_trajectory(
    ic = 0.4, map = "logistic", n_steps = 300, r = r
  )
  points(rep(r, 50), traj[1, 251:300], pch = ".", cex = 0.5)
}

Harmonic Time Averages and Mesochronic Analysis

library(koopman.dmd)

# -- Single-point HTA --
hta <- harmonic_time_average(
  ic = c(0.1, 0.2),
  map = "standard",
  observable = "sin_pi",
  omega = 0.0,
  n_iter = 20000,
  epsilon = 0.12
)
cat("HTA magnitude:", hta$magnitude, "\n")
cat("HTA phase:    ", hta$phase, "\n")

# -- Mesochronic map --
meso <- mesochronic_compute(
  map = "standard",
  x_range = c(0, 1),
  y_range = c(0, 1),
  resolution = 300,
  observable = "cos_pi",
  omega = 0.0,
  n_iter = 10000,
  epsilon = 0.12
)

image(
  meso$x_coords, meso$y_coords, meso$hta_matrix,
  col = hcl.colors(256, "viridis"),
  xlab = "x", ylab = "y",
  main = "Mesochronic Map (cos_pi, omega=0)"
)

# -- Phase space classification --
cls <- classify_phase_space(
  ic = c(0.5, 0.5),
  map = "standard",
  omega = 0.0,
  n_iter = 50000,
  epsilon = 0.97
)
cat("Classification:", cls$classification, "\n")

# -- Convergence tracking --
conv <- hta_convergence(
  ic = c(0.1, 0.2),
  map = "standard",
  observable = "quadratic",
  omega = 0.0,
  max_iter = 100000,
  step = 1000,
  epsilon = 0.12
)
plot(conv$iteration, conv$magnitude,
     type = "l", log = "x",
     xlab = "Iteration", ylab = "HTA Magnitude",
     main = "HTA Convergence")