Getting started with koopman.dmd

What this package does

Dynamic Mode Decomposition (DMD) takes a sequence of measurements from a dynamical system and fits a linear operator that advances the state one step in time. Diagonalising that operator decomposes the data into modes, each with its own frequency and growth or decay rate.

The connection to Koopman operator theory is what makes this more than a linear method. The Koopman operator describes how observables of a system evolve, and it is linear even when the underlying dynamics are not — at the cost of acting on an infinite-dimensional space. DMD and its extensions here are finite-dimensional approximations of that operator, which is why they can capture nonlinear behaviour that a naive linear fit would miss.

The numerics are implemented in Rust and reached through extendr. No BLAS or LAPACK installation is needed.

Data layout

Every function takes data as a matrix with one row per variable and one column per time step. This is the transpose of the usual “tidy” layout, and it is the most common source of confusing results.

t <- seq(0, 10, length.out = 200)
X <- rbind(sin(t), cos(t))
dim(X)   # 2 variables, 200 time steps
#> [1]   2 200

Core DMD

d <- dmd(X, rank = 2, dt = t[2] - t[1])
summary(d)
#> Dynamic Mode Decomposition
#>   Rank: 2
#>   Data dimensions: 2 variables x 200 time steps
#>   Centered: FALSE
#>   Time step: 0.0502513
#>   Singular values: 10.2428, 9.6997
#>   Eigenvalue magnitudes: 1, 1

The spectrum reports each mode’s frequency, growth rate and amplitude. For this pure oscillation the modes sit essentially on the unit circle, meaning neither growth nor decay:

dmd_spectrum(d)
#>   index magnitude       phase frequency   period   growth_rate amplitude
#> 1     0         1  0.05025126 0.1591549 6.283185 -2.209344e-15       0.5
#> 2     1         1 -0.05025126 0.1591549 6.283185 -2.209344e-15       0.5
#>   stability
#> 1   neutral
#> 2   neutral

dmd_stability() summarises that directly. A spectral radius at 1 is a system that neither grows nor decays:

dmd_stability(d)
#> $is_stable
#> [1] TRUE
#> 
#> $is_unstable
#> [1] FALSE
#> 
#> $is_marginal
#> [1] TRUE
#> 
#> $spectral_radius
#> [1] 1

pred <- predict(d, n_ahead = 50)
dim(pred)
#> [1]  2 50

plot(t, X[1, ], type = "l", xlab = "time", ylab = "x1",
     xlim = c(0, 13), main = "DMD forecast")
t_future <- seq(max(t) + (t[2] - t[1]), by = t[2] - t[1], length.out = 50)
lines(t_future, pred[1, ], col = "red", lwd = 2)
legend("bottomleft", c("observed", "forecast"), col = c("black", "red"),
       lty = 1, bty = "n")

dmd_error(d)
#> $rmse
#> [1] 1.002686e-14
#> 
#> $mae
#> [1] 8.01883e-15
#> 
#> $mape
#> [1] 3.017636e-12
#> 
#> $relative_error
#> [1] 1.418013e-14
dmd_residual(d)
#> $residual_norm
#> [1] 7.38836e-15
#> 
#> $residual_relative
#> [1] 5.23747e-16

DMD with control: when the system is driven

Standard DMD assumes the system evolves on its own. If the data comes from a system driven by a measured input — an actuated mechanical system, a circuit with an applied voltage — plain DMD folds the forcing into the identified operator, biasing it. dmdc() (Proctor, Brunton and Kutz 2016) separates the two by identifying the forced linear system x_{t+1} = A x_t + B u_t.

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 input applied during each transition, so columns may come from many concatenated trajectories.

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 (i in seq_len(m)) {
  u_i <- sin(0.7 * (i - 1)) + 0.5 * cos(2.3 * (i - 1) + 1)
  X1[, i] <- x
  U[, i] <- u_i
  x <- as.numeric(A0 %*% x + B0 * u_i)
  X2[, i] <- x
}

fit <- dmdc(X1, X2, U, rank_input = 3)
round(fit$a, 6)   # recovers A0
#>      [,1] [,2]
#> [1,]  0.9  0.1
#> [2,]  0.0  0.8
round(fit$b, 6)   # recovers B0
#>      [,1]
#> [1,]  0.5
#> [2,]  1.0

Joint identification needs the input to be persistently exciting and exogenous — if u is computed from the state (feedback), the regression cannot separate A from the feedback path. In that case, or whenever the input coupling is known by construction, pin it with known_B and only A is estimated:

fit2 <- dmdc(X1, X2, U, rank_input = 2, known_B = B0)
round(fit2$a, 6)
#>      [,1] [,2]
#> [1,]  0.9  0.1
#> [2,]  0.0  0.8

The eigenvalues of the identified operator describe the unforced dynamics, and predict() steps the fitted system under any input sequence:

dmdc_stability(fit)
#> $is_stable
#> [1] TRUE
#> 
#> $is_unstable
#> [1] FALSE
#> 
#> $is_marginal
#> [1] FALSE
#> 
#> $spectral_radius
#> [1] 0.9
pred <- predict(fit, U = U)      # replaying the training input reproduces X2
max(abs(pred - X2))
#> [1] 1.776357e-15

With U = NULL, dmdc() fits an autonomous model from explicit pairs — useful when the data is many short trajectories from different initial conditions, which dmd() cannot digest.

Extended DMD: when the dynamics are nonlinear

Standard DMD fits a linear operator to the measured coordinates. If the dynamics are nonlinear in those coordinates, that fit will be poor. Lifting maps the data into a richer set of observables where the dynamics are closer to linear — a direct, finite-dimensional stand-in for what the Koopman operator does exactly.

Consider a signal whose second component is a quadratic function of the first:

t2 <- seq(0, 10, length.out = 200)
Xn <- rbind(sin(t2), sin(t2)^2)

plain  <- dmd(Xn, rank = 2)
lifted <- dmd(Xn, lifting = "polynomial", lifting_param = 2)

plain_err  <- dmd_error(plain)
lifted_err <- dmd_error(lifted)
c(plain = plain_err$rmse, lifted = lifted_err$rmse)
#>     plain    lifted 
#> 0.6436273 0.6436273

Available lifting functions are "polynomial", "polynomial_cross", "trigonometric" and "delay", with lifting_param giving the degree, the number of harmonics, or the number of delays.

Hankel-DMD: one measured variable

Often only a single scalar signal is observed. Time-delay embedding reconstructs a higher-dimensional state from that one series, which is enough for DMD to work with — this is Takens’ embedding idea applied to the Koopman setting.

tt <- seq(0, 4 * pi, length.out = 200)
y  <- matrix(sin(tt), nrow = 1)   # a single row

h <- hankel_dmd(y, delays = 20)
h
#> HankelDMD(rank=2, delays=20, n_obs=1)

Leaving rank = NULL lets the rank be chosen automatically. Asking for more modes than the signal supports — a pure sinusoid supports two — leaves the reduced operator near-singular and the eigendecomposition can fail to converge.

pred_h <- predict(h, n_ahead = 20)
dim(pred_h)
#> [1]  1 20

Generalized Laplace Analysis

GLA computes Koopman eigenfunctions directly, through weighted time averages, rather than by diagonalising a fitted operator. It is a useful cross-check on a DMD result, since the two arrive at the spectrum by different routes.

g <- gla(X, n_eigenvalues = 2)
g
#> GLA(n_eigenvalues=2, n_obs=2, n_time=200)

Phase space analysis

The second half of the package addresses a different question. For area-preserving maps, the interest is not forecasting but classifying orbits: which initial conditions lead to regular motion, and which to chaos.

Several standard maps are built in:

traj <- generate_trajectory("standard", c(0.1, 0.2), 2000, epsilon = 0.9)
dim(traj)
#> [1]    2 2001

plot(traj[1, ], traj[2, ], pch = ".", xlab = "x", ylab = "p",
     main = "Chirikov standard map, epsilon = 0.9")

regular <- harmonic_time_average("standard", c(0.5, 0.0), "sin_pi", 0.5, 2000,
                                 epsilon = 0.9)
chaotic <- harmonic_time_average("standard", c(0.1, 0.2), "sin_pi", 0.5, 2000,
                                 epsilon = 0.9)
c(regular = regular$magnitude, chaotic = chaotic$magnitude)
#>      regular      chaotic 
#> 0.0074984158 0.0008363106

conv <- hta_convergence("standard", c(0.5, 0.0), "sin_pi", 0.5, 2000,
                        epsilon = 0.9)
conv$dynamics_type
#> [1] "resonating_periodic"

mhp <- mesochronic_compute("standard", c(0, 1), c(0, 1), 40, "sin_pi", 0.5, 300,
                           epsilon = 0.9)
image(mhp$x_coords, mhp$y_coords, mhp$hta_matrix,
      col = hcl.colors(64, "YlGnBu", rev = TRUE),
      xlab = "x", ylab = "p", main = "Mesochronic harmonic plot")

labels <- classify_phase_space(as.vector(mhp$hta_matrix))
table(factor(labels, levels = 1:3,
             labels = c("resonating", "chaotic", "non-resonating")))
#> 
#>     resonating        chaotic non-resonating 
#>            721            866             13

Where to go next

• ?dmd, ?dmdc and ?hankel_dmd for the full argument lists
• ?mesochronic_compute for the phase space tools
• The Rust crate koopman-dmd and the Python package koopman-dmd expose the same functionality

References

Schmid, P.J. (2010). Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics, 656, 5–28. https://doi.org/10.1017/S0022112010001217

Proctor, J.L., Brunton, S.L., & Kutz, J.N. (2016). Dynamic Mode Decomposition with Control. SIAM Journal on Applied Dynamical Systems, 15(1), 142–161. https://doi.org/10.1137/15M1013857

Kutz, J.N., Brunton, S.L., Brunton, B.W., & Proctor, J.L. (2016). Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems. SIAM. https://doi.org/10.1137/1.9781611974508

Mezić, I. (2020). Spectrum of the Koopman operator, spectral expansions in functional spaces, and state-space geometry. https://doi.org/10.48550/arXiv.2009.05883

Levnajić, Z. & Mezić, I. (2014). Ergodic theory and visualization. https://doi.org/10.48550/arXiv.0808.2182