RasmahRasmah

Advanced Model Reduction

Overview

POD and DMD are the workhorses; this tutorial covers the layer on top — the reduced models that make a ROM accurate, cheap and data-driven:

  • Quadratic manifolds — a nonlinear (x = Vx̃ + V̄(x̃⊗x̃)) reduced manifold.
  • Tensor trains — a multi-dimensional array factorised into a chain of small cores.
  • RBF interpolation — a mesh-free surrogate s(x) = Σ wⱼ φ(‖x−xⱼ‖).
  • Hyper-reduction — nonnegative least-squares (ECSW/EQP) sampling weights.
  • Balanced truncation — an H₂-optimal linear ROM via Hankel singular values.
  • DMD variants — total-least-squares and forward-backward DMD.
using Rasmah
using LinearAlgebra

Quadratic manifold

Fit a quadratic manifold to snapshots and check decode ∘ encode reproduces them:

X = [sin.(range(0, 2π; length = 50))'; cos.(range(0, 2π; length = 50))']
qm = quadratic_manifold(X, 2)
x = X[:, 1]
println("quadratic-manifold round-trip: ", norm(decode(qm, encode(qm, x)) - x) < 1e-8)
quadratic-manifold round-trip: true

Tensor train

tt_svd factorises a tensor into a chain of cores; tt_full reconstructs it:

A = rand(4, 5, 6)
tt = tt_svd(A; tol = 1e-12)
println("tensor-train round-trip: ", norm(tt_full(tt) - A) < 1e-10)
println("TT ranks: ", tt_rank(tt))
tensor-train round-trip: true
TT ranks: [4, 6]

RBF interpolation

Fit a Gaussian radial-basis interpolant to a few samples of sin and evaluate it at a new point:

Xr = reshape(collect(0.0:0.5:2.0), 1, :)
Yr = reshape(sin.(Xr), 1, :)
rbf = radial_basis_interpolation(Xr, Yr; kernel = :gaussian, shape = 1.0)
println("RBF at x=1.0: ", round(rbf([1.0])[1]; digits = 4), "  (sin(1) = ", round(sin(1.0); digits = 4), ")")
RBF at x=1.0: 0.8415  (sin(1) = 0.8415)

Hyper-reduction (ECSW / EQP weights)

Find sparse, nonnegative weights w on the rows of a basis Φ whose weighted sum reproduces the exact integral Φᵀ1:

Φ = rand(20, 3)
w = eqp_weights(Φ)
println("EQP reproduces integrals: ", norm(Φ' * w - Φ' * ones(20)) < 1e-8)
println("weights nonnegative & sparse: ", minimum(w) >= 0, ", ", count(>(0), w), " nonzero of 20")
EQP reproduces integrals: true
weights nonnegative & sparse: true, 3 nonzero of 20

Balanced truncation

Balanced truncation returns an H₂-optimal reduced linear system plus its Hankel singular values (the energy coupling of each state):

A = [-0.5 1.0 0.0; -1.0 -0.5 0.0; 0.0 0.0 -2.0]
B = reshape([1.0, 0.0, 1.0], 3, 1)
C = reshape([1.0 0.0 1.0], 1, 3)
res = balanced_truncation(A, B, C; r = 2)
println("Hankel singular values: ", round.(res.hsv; digits = 3))
println("reduced A is ", size(res.Ar))
Hankel singular values: [0.761, 0.349, 0.038]
reduced A is (2, 2)

DMD variants

Total-least-squares and forward-backward DMD recover the Ritz eigenvalues of a linear system more robustly than plain DMD:

Ad = [-0.2 3.0; -3.0 -0.2]
Xd = zeros(2, 60); Xd[:, 1] = [1.0, 0.0]
for k in 1:59
    Xd[:, k + 1] = Xd[:, k] + 0.05 * (Ad * Xd[:, k])
end
Yd = Xd[:, 2:end]
λt = tls_dmd(Xd[:, 1:end-1], Yd; rank = 2).eigenvalues
λf = fb_dmd(Xd[:, 1:end-1], Yd; rank = 2).eigenvalues
λtrue = eigvals(I + 0.05 * Ad)                     # discrete-time Koopman eigenvalues
println("true eigenvalues:    ", round.(λtrue; digits = 3))
println("TLS-DMD eigenvalues: ", round.(λt; digits = 3))
println("FB-DMD eigenvalues:  ", round.(λf; digits = 3))
true eigenvalues:    ComplexF64[0.99 - 0.15im, 0.99 + 0.15im]
TLS-DMD eigenvalues: ComplexF64[0.99 - 0.15im, 0.99 + 0.15im]
FB-DMD eigenvalues:  ComplexF64[0.99 - 0.15im, 0.99 + 0.15im]

Next steps