Model reduction
Overview
A simulation field — temperature over a mesh, velocity in a fluid — can have millions of degrees of freedom, yet the interesting behaviour often lives in a handful of dominant patterns. Model reduction finds those patterns and replaces the expensive high-dimensional system with a cheap low-dimensional one.
The two tools in this chapter are the workhorses of the field. Proper orthogonal decomposition (POD) compresses a family of fields into a few basis vectors; dynamic mode decomposition (DMD) reads the dynamics of a system straight out of its snapshots. Both reduce the cost of solving, storing, and predicting — at the price of a small, controllable approximation error.
Compressing snapshots with POD
Rasmah.ModelReduction.pod — Function
pod(S; num_basis=nothing, energy=nothing, method=:svd, seed=0, oversampling=10, power=2) -> PODResultProper Orthogonal Decomposition of the snapshot matrix S (columns are snapshots, $n \times m$).
Theory
A family of fields — snapshots of a simulation, frames of a video — usually lives near a low-dimensional subspace: the columns of $S$ are nearly a linear combination of a few dominant directions. POD finds those directions as the leading left singular vectors of $S$. The retained energy $\sum_{i \le k} \sigma_i^2 / \sum_i \sigma_i^2$ measures how much of the data the $k$-mode basis explains, which is why POD is the workhorse of reduced-order modeling: project a large field onto a handful of modes and solve in $k$ dimensions instead of $n$.
Arguments
S: the snapshot matrix (columns are snapshots).
Keyword options
num_basis = nothing: fix the rank $k$;energyselects instead the smallest rank retaining that fraction of $\sum \sigma^2$.method = :svd::svd(full),:randomized(Halko–Martinsson–Tropp), or:lanczos(svdsfor very tallS).
Returns
A PODResult with basis (n×k orthonormal columns), singular_values, energy (retained fraction), and num_basis.
Example
julia> using Rasmah
julia> xs = range(-1, 1; length = 50);
julia> S = hcat([exp.(-((xs .- c) ./ 0.2) .^ 2) for c in range(-0.5, 0.5; length = 10)]...);
julia> pod(S; energy = 0.999).num_basis
7xs = range(-1, 1; length = 50)
S = hcat([exp.(-((xs .- c) ./ 0.2) .^ 2) for c in range(-0.5, 0.5; length = 10)]...)
r = pod(S; energy = 0.999)
(r.num_basis, round(r.energy; digits = 4))(7, 0.9994)A family of ten 50-dimensional snapshots collapses to 7 modes while retaining 99.9% of the energy — a seven-fold reduction. Projecting onto fewer modes and measuring the leftover is one line:
Rasmah.ModelReduction.pod_error — Function
pod_error(S, Φ) -> RealRelative projection error ‖S − ΦΦᵀS‖/‖S‖ of the basis Φ on S (Frobenius).
Φ = pod_basis(S; num_basis = r.num_basis)
round(pod_error(S, Φ); digits = 4)0.0248Reading dynamics with DMD
Rasmah.ModelReduction.dmd — Function
dmd(X, Y; rank=nothing) -> DMDResultExact dynamic mode decomposition of the snapshot pair (X, Y) (Y is X shifted one step). Returns Ritz eigenvalues plus exact DMD modes $\Phi = Y V \Sigma^{-1} W$ of the rank-rank (default full) subspace.
Theory
DMD assumes the snapshots evolve under an unknown linear operator, $y_k = A x_k$, and finds the best low-rank $A$ from the data. The eigenvalues of the reduced operator $\tilde A = U^T Y V \Sigma^{-1}$ are the growth/decay rates and frequencies of the dynamics, and the DMD modes are the associated spatial structures — so a pair of shifted snapshots is enough to read off the spectrum of the system that produced them.
Arguments
X: snapshots $[x_1 \ \cdots \ x_m]$.Y: the same snapshots shifted one step $[x_2 \ \cdots \ x_{m+1}]$.
Keyword options
rank = nothing: truncate to arank-dimensional subspace (default: the numerical rank ofX).
Returns
A DMDResult with modes, eigenvalues, singular_values, amplitudes, and rank.
Example
julia> using Rasmah
julia> A = [0.0 -1.0; 1.0 0.0];
julia> X = [1.0 0.0 -1.0; 0.0 1.0 0.0];
julia> round.(dmd(X, A * X; rank = 2).eigenvalues; digits = 4) == ComplexF64[0.0 - 1.0im, 0.0 + 1.0im]
trueSee also: pod, dmd_predict.
Rasmah.ModelReduction.dmd_predict — Function
dmd_predict(d::DMDResult, b, t::Integer) -> VectorPredict the state at (integer) step t from amplitudes b: x(t) = Σᵢ φᵢ λᵢᵗ bᵢ.
Rasmah.ModelReduction.dmd_reconstruct — Function
dmd_reconstruct(d::DMDResult, x0, m::Integer) -> MatrixReconstruct the first m snapshots [x₀ … x_{m-1}] from the initial state x₀ using the DMD model.
Rasmah.ModelReduction.dmd_amplitudes — Function
dmd_amplitudes(d::DMDResult, x0) -> VectorLeast-squares amplitudes b = Φ⁺ x₀ so that x₀ = Φ b.
A = [0.0 -1.0; 1.0 0.0] # a rotation: eigenvalues ±i
X = [1.0 0.0 -1.0; 0.0 1.0 0.0] # three snapshots of the orbit
d = dmd(X, A * X; rank = 2)
round.(d.eigenvalues; digits = 4)2-element Vector{ComplexF64}:
0.0 - 1.0im
0.0 + 1.0imThe eigenvalues land exactly on $\pm i$ — DMD recovered the rotation frequency of the system from three snapshots alone, without ever seeing the matrix A. That is the whole promise of data-driven model reduction: the snapshots are enough.
Non-intrusive surrogates
When you cannot project (a black-box code, an experiment), fit a surrogate instead — a function that reproduces the input-output map from samples:
Rasmah.ModelReduction.radial_basis_interpolation — Function
radial_basis_interpolation(X, Y; kernel=:gaussian, shape=1.0) -> RBFInterpolantFit s(x) = Σⱼ wⱼ φ(‖x − xⱼ‖) to the training data X (d×m, parameters as columns) and Y (q×m, outputs as columns), solving W Φ = Y for the interpolation matrix Φᵢⱼ = φ(‖xᵢ − xⱼ‖).
X = [0.0 0.5 1.0 1.5 2.0]
Y = [0.0 0.25 1.0 2.25 4.0] # y = x^2
rbf = radial_basis_interpolation(X, Y; kernel = :thin_plate)
round.(rbf([1.25]); digits = 4)1-element Vector{Float64}:
1.5783The radial-basis interpolant reproduces the samples exactly and interpolates smoothly between them, so it can stand in for $x^2$ where no closed form exists.
Next steps
Reduced models still need a gradient and an optimizer to be useful — see Optimization — and a neural network is one of the most flexible surrogates of all, covered in Neural networks.