RasmahRasmah

Model Reduction: POD and DMD

println("Julia ", VERSION)
Julia 1.12.6

Overview

A fine finite-element or finite-volume simulation produces a state vector with millions of unknowns at thousands of time steps. Analysing, optimising, or running that model in real time is often impossible — which is where reduced-order modelling (ROM) comes in. The idea: even though the state lives in a huge space, the physically relevant dynamics usually live in a low-dimensional subspace. If you can find that subspace, you can project the problem onto it and get a small, fast model with almost the same accuracy.

Rasmah implements the two workhorse methods natively:

  • Proper orthogonal decomposition (POD) — find an optimal low-rank basis of snapshots (spatial fields), by a singular-value decomposition (SVD).
  • Dynamic mode decomposition (DMD) — find the best-fit linear operator that maps one snapshot to the next, and read frequencies / growth rates off its eigenvalues.

Proper orthogonal decomposition (POD)

Take a collection of snapshots S (each column is one spatial field, n×m). POD computes the SVD S = U Σ Vᵀ and keeps the first k left singular vectors as an orthonormal basis Φ of the dominant directions of variation. The squares of the singular values are the energy (variance) carried by each mode, so the sum of the retained σᵢ² over the total Σ σᵢ² is the fraction of energy captured.

Here the snapshots are a smooth, low-rank family — Gaussian bumps whose centre μ slides from -1 to 1 — so a handful of modes should capture almost all the energy:

using Rasmah
using LinearAlgebra

x = range(-2, 2, length=200)
μs = range(-1, 1, length=30)
S = reduce(hcat, [exp.(-((x .- μ) ./ 0.35).^2) for μ in μs])   # 200 × 30

r = pod(S; energy=0.9999)
println("modes for 99.99% energy = ", r.num_basis, "  (out of ", size(S, 2), ")")
modes for 99.99% energy = 9  (out of 30)

The energy spectrum drops off fast for a smooth family, and the retained fraction is exactly what we asked for:

println("retained energy fraction = ", round(r.energy, digits=6))
println("relative projection error = ", round(pod_error(S, r.basis), digits=5))
println("analytic  √(1 − energy)   = ", round(sqrt(1 - r.energy), digits=5))
retained energy fraction = 0.999908
relative projection error = 0.00957
analytic  √(1 − energy)   = 0.00957

The projection error ‖S − ΦΦᵀS‖/‖S‖ equals √(1 − energy) for an SVD basis — so the check is exact, not approximate.

Dynamic mode decomposition (DMD)

DMD models the sequence X as x_{k+1} = A x_k and finds the best-fit A = Y X⁺ (where Y is X shifted one step) restricted to a rank-r POD subspace. The eigenvalues of A are the Ritz values λ: their angle gives the frequency, their magnitude the growth/decay rate. For a signal sampled from a continuous-time system with eigenvalue λ_c, the discrete eigenvalue is exactly λ = exp(λ_c Δt).

A decaying, rotating two-state system d/dt [c; s] = [[σ,−ω],[ω,σ]] [c; s] is the perfect test, because we know the answer in closed form:

σ, ω, dt = -0.1, 5.0, 0.05
ts = collect(0:dt:5)
X = reduce(hcat, [[exp(σ*t)*cos(ω*t), exp(σ*t)*sin(ω*t)] for t in ts])  # 2 × m

d = dmd(X[:, 1:end-1], X[:, 2:end])
println("DMD eigenvalues = ", round.(d.eigenvalues, digits=5))
println("analytic  exp((σ ± iω)·dt) = ", round(exp((σ + im*ω)*dt), digits=5))
DMD eigenvalues = ComplexF64[0.96408 - 0.24617im, 0.96408 + 0.24617im]
analytic  exp((σ ± iω)·dt) = 0.96408 + 0.24617im

The modes reconstruct the trajectory exactly:

rec = dmd_reconstruct(d, X[:, 1], size(X, 2))
println("reconstruction error = ", round(norm(rec - X) / norm(X), digits=6))
reconstruction error = 0.0

Hyper-reduction (DEIM)

Once you have a POD basis you often need to evaluate a nonlinear field only at a few points rather than everywhere. DEIM picks those points greedily so that interpolating from them reproduces the whole field through the basis:

p = deim(r.basis)                              # sample-point indices
f = S[:, 15]                                   # one snapshot, in full
frec = deim_interpolate(r.basis, p, f[p])      # reconstructed from |p| points
println("DEIM points used = ", length(p))
println("DEIM interpolation error = ", round(norm(frec - f) / norm(f), digits=5))
DEIM points used = 9
DEIM interpolation error = 0.01002

greedy_basis, optdmd, balanced truncation, operator inference, and the non-intrusive surrogates build on these two primitives — see the API reference for the full catalogue.

Next steps