The Reduced Basis Method: Snapshots, Offline/Online, and Certified Errors

The last two posts built the theory: well-posedness, coercivity, inf-sup stability, supremizers, and quasi-optimal truth approximations. This post is where it all pays off — turning that theory into an actual reduced-order model that’s both fast and certified.

From a truth system to a genuinely small one

The truth approximation \(u_\mu^\delta\) from the last post lives in a finite-dimensional space \(\mathcal U^\delta\) of dimension \(N\) — but \(N\) can easily be in the hundreds of thousands for a fine mesh. The whole point of the reduced basis method (RBM) is to replace this with a space of dimension \(N_{\text{rb}} \ll N\) — often just a few dozen — built not from generic mesh functions, but from actual solutions of the problem at a handful of carefully chosen parameter values.

Building the reduced spaces

The reduced trial space is easy to describe: solve the truth problem at a set of training parameters \(\mu^{(1)}, \dots, \mu^{(N_{\text{rb}})}\), keep each solution as a snapshot \(\xi_n := u_{\mu^{(n)}}^\delta\), and set

\[\mathcal U_N := \mathrm{span}\{\xi_1, \dots, \xi_{N_{\text{rb}}}\}.\]

The test space \(\mathcal V_N(\mu)\) is where the danger hides. The naive choice is Galerkin, \(\mathcal V_N = \mathcal U_N\) — reuse the same snapshots as test functions. For non-coercive problems, this can destroy inf-sup stability entirely.

🤔 Why does plain Galerkin fail here?

For a coercive problem, \(b(u,u) \ge \alpha\|u\|^2\) means the best possible test partner for any \(u\) is \(u\) itself. Reduce the trial space to 10 snapshots, set the test space equal to the same 10 snapshots, and stability is inherited for free — whatever \(u_N \in \mathcal U_N\) you pick, its perfect partner (\(u_N\) itself) is trivially already in the test space. This is exactly Remark 1.3.4 from the last post, now one dimension smaller.

For a non-coercive problem, there’s no such shortcut: \(b(u,u)\) can be zero, or even negative. Each \(u\) needs a genuinely different partner, its supremizer \(s_u\), which can point in a completely different direction — e.g. because the operator involves derivatives that shift or skew the function. There is no mathematical reason the supremizer of a snapshot should be expressible as a linear combination of the snapshots themselves. If you set \(\mathcal V_N = \mathcal U_N\) anyway, you’re implicitly demanding that every trial function find its test partner inside a space that was never built to contain it. The supremum collapses:

\[\sup_{v_N \in \mathcal V_N} \frac{b(u_N, v_N)}{\|v_N\|} \approx 0 \quad\Longrightarrow\quad \beta_N \approx 0.\]

On the computer, this shows up as the tiny reduced matrix \(\boldsymbol B_N\) becoming singular or wildly ill-conditioned — even though the underlying physical problem was perfectly well-behaved. The instability was introduced entirely by a careless choice of subspace.

The fix is exactly the supremizer machinery from the earlier post: for every snapshot \(\xi_n\), explicitly compute its supremizer \(\eta_n(\mu) := s_{\xi_n}(\mu)\), and define

\[\mathcal V_N(\mu) := \mathrm{span}\{\eta_1(\mu), \dots, \eta_{N_{\text{rb}}}(\mu)\}.\]

Since \(\mathcal U_N \ne \mathcal V_N\), this is a Petrov–Galerkin method. Every trial function’s perfect test partner is now guaranteed to already be sitting in the test space by construction, so \(\beta_N \ge \beta_\mu^\circ > 0\) is preserved, the reduced matrix stays invertible, and Xu–Zikatanov (from the last post) guarantees a quasi-optimal reduced solution.

Keeping everything affine: the empirical interpolation method

Offline/online efficiency (precompute once, reuse cheaply for every new \(\mu\)) depends on the operator being affine in the parameter,

\[b(u,v;\mu) = \sum_{q=1}^{Q_b} \vartheta_q(\mu)\, b_q(u,v),\]

exactly as in the first post. But what if it isn’t? Suppose a coefficient looks like \(b(u,v; e^{\mu x})\) — here \(\mu\) and \(x\) are tangled together inside an exponential, and there’s no way to split them into a sum of (parameter function) \(\times\) (spatial function). Nothing here can be precomputed offline as-is.

Definition 1.4.3 — Empirical Interpolation Method (EIM)

Given a non-affine function \(g(x;\mu)\), the EIM approximates it by an affine expansion

\[g(x;\mu) \;\approx\; \sum_{q=1}^{Q} \vartheta_q(\mu)\, g_q(x),\]

where both the basis functions \(g_q\) and the interpolation points \(x_q\) (“magic points”) are chosen by a greedy procedure: at each step, the function currently making the largest interpolation error determines the next magic point and basis function.

The key difference from classical polynomial interpolation: EIM doesn’t fix the interpolation points in advance — it discovers them, adapted to the specific function being approximated.

Theorems 1.4.10 & 1.4.12

The EIM approximation error is bounded by the best possible interpolation error, amplified by the Lebesgue constant \(\Lambda_M\) (where \(M\) is the number of magic points). This makes EIM quasi-optimal, and for analytic functions \(g\), the convergence is exponential.

Notice the shape of that guarantee: best possible error, amplified by a stability-type constant. It’s the same pattern as Xu–Zikatanov’s quasi-optimality bound from the last post — there the amplifier was \(C_\mu/\beta_\mu^\circ\), here it’s \(\Lambda_M\). Whenever you see “best approximation times a controllable stability factor,” it’s this same idea showing up again.

Assembling the reduced system, offline and online

🤔 How do you actually get from the weak form to a matrix?

The reduced solution lives in \(\mathcal U_N\), so it’s exactly a linear combination of the snapshots:

\[u_\mu^N = \sum_{n=1}^{N_{\text{rb}}} \alpha_n(\mu)\, \xi_n.\]

The snapshots \(\xi_n\) are already known (computed offline). The only unknowns left to find online are the coefficients \(\alpha_n(\mu)\).

Substitute this ansatz into the reduced weak formulation \(b_\mu(u_\mu^N, v_N) = f_\mu(v_N)\) for all \(v_N \in \mathcal V_N(\mu)\):

\[b_\mu\Big(\sum_{n=1}^{N_{\text{rb}}} \alpha_n(\mu)\,\xi_n,\; v_N\Big) = f_\mu(v_N).\]

Since \(b_\mu\) is linear in its first argument, pull the sum and the scalars \(\alpha_n(\mu)\) out front:

\[\sum_{n=1}^{N_{\text{rb}}} \alpha_n(\mu) \cdot b_\mu(\xi_n, v_N) = f_\mu(v_N).\]

This has to hold for every test function \(v_N\) — but it’s enough to check it on the basis of the test space. Testing against each \(\eta_{n'}(\mu)\) in turn (\(n' = 1,\dots, N_{\text{rb}}\)) gives \(N_{\text{rb}}\) separate scalar equations:

\[\sum_{n=1}^{N_{\text{rb}}} \alpha_n(\mu) \cdot b_\mu(\xi_n, \eta_{n'}(\mu)) = f_\mu(\eta_{n'}(\mu)) \qquad \text{for } n' = 1,\dots,N_{\text{rb}}.\]

That’s precisely a matrix-vector multiplication: the entry in row \(n'\), column \(n\) of a matrix \(\boldsymbol B_N(\mu)\) is \(b_\mu(\xi_n, \eta_{n'}(\mu))\), the entry in row \(n'\) of a vector \(\mathbf f_N(\mu)\) is \(f_\mu(\eta_{n'}(\mu))\), and the unknown vector is \(\boldsymbol\alpha_N = (\alpha_1,\dots,\alpha_{N_{\text{rb}}})^\top\). So the reduced problem is exactly the linear system

\[\boldsymbol B_N(\mu)\, \boldsymbol\alpha_N = \mathbf f_N(\mu).\]

This is the exact same construction as assembling a finite element system — the only difference is what the basis functions are. In FEM, \(\xi_n\) and \(\eta_{n'}\) are tiny “hat” functions on a mesh, and the system is enormous (think \(100{,}000 \times 100{,}000\)). In RBM, \(\xi_n\) are global snapshot solutions and \(\eta_{n'}\) are their supremizers, and the system is tiny (think \(20 \times 20\)).

Because \(\eta_{n'}(\mu)\) depends on \(\mu\), and because everything (the bilinear form, the snapshots, the supremizers) is affine or made affine via EIM, \(\boldsymbol B_N(\mu)\) itself decomposes affinely:

\[\boldsymbol B_N(\mu) = \sum_{q=1}^{Q_b} \vartheta_q(\mu)\, \boldsymbol B_N^{(q)}, \qquad \boldsymbol B_N^{(q)} := \big[b_q(\xi_n, \eta_{n'})\big]_{n,n'}.\]

• Offline (once, expensive): evaluate every \(b_q(\xi_n, \eta_{n',q'})\) and store the small \(N_{\text{rb}} \times N_{\text{rb}}\) matrices \(\boldsymbol B_N^{(q)}\). This costs roughly \(\mathcal O(N_{\text{rb}}^2 N)\) — expensive, but paid exactly once.
• Online (per query, cheap): for a new \(\mu\), just weight and sum the precomputed small matrices by the scalars \(\vartheta_q(\mu)\), costing \(\mathcal O(Q_b N_{\text{rb}}^2)\), then solve the tiny system at \(\mathcal O(N_{\text{rb}}^3)\).

A direct truth solve costs \(\mathcal O(N)\) — linear, but in a huge dimension \(N\). RBM online costs \(\mathcal O(N_{\text{rb}}^3)\) — cubic, but in a tiny dimension \(N_{\text{rb}} \ll N\). Drag the slider below to see how many repeated queries it takes for the one-time offline investment to pay off.

The error estimator: the heart of the reduced basis method

All of this speed is worthless without a guarantee. We want a computable bound

\[\|u_\mu^\delta - u_\mu^N\|_{\mathcal U} \;\le\; \Delta_N(\mu),\]

and, echoing Lemma 1.3.2 from the last post almost exactly one level down,

\[\Delta_N(\mu) = \frac{\|r_N^\delta(\mu)\|_{\mathcal V'}}{\beta_\mu^\circ}.\]

Two challenges stand between this formula and something you can actually compute online.

Challenge A — the numerator. The dual norm \(\|r_N^\delta(\mu)\|_{\mathcal V'}\) is itself defined as a supremum over an infinite (or at least huge) set of test functions — not something you evaluate directly. The fix is the Riesz representation theorem: find the representer \(\hat r_N(\mu) \in \mathcal V^\delta\) solving

\[(\hat r_N(\mu), v)_{\mathcal V} = r_N^\delta(\mu)(v) \qquad \text{for all } v \in \mathcal V^\delta\]

— note this is exactly the same shape of equation that defined the supremizer earlier, and the same Cauchy–Schwarz argument gives \(\|r_N^\delta(\mu)\|_{\mathcal V'} = \|\hat r_N(\mu)\|_{\mathcal V}\). Because the residual \(r_N^\delta(\mu) = f_\mu - B_\mu u_\mu^N\) is affine in \(\mu\) (both \(f_\mu\) and \(B_\mu\) are affine, and \(u_\mu^N\) depends on \(\mu\) only through the scalar coefficients \(\alpha_n(\mu)\)), the Riesz representer decomposes affinely too — its \(\mu\)-independent pieces can be precomputed offline, and online you only assemble a small linear combination of them.

💡 Remark 1.4.17 — the square-root effect

To get the actual norm \(\|\hat r_N(\mu)\|_{\mathcal V}\), you assemble a small \(Q_r \times Q_r\) inner-product computation from the precomputed pieces, and then take a square root at the very end. That square root isn’t free numerically: it costs you roughly half your significant digits of floating-point precision. Worth remembering if an error estimate ever looks suspiciously imprecise — the culprit is often this last square root, not a bug elsewhere.

Challenge B — the denominator. We also need a cheap, online, lower bound for \(\beta_\mu^\circ\) — inverting the full operator to get it exactly would cost \(\mathcal O(N^3)\), defeating the whole purpose. The standard fix is the Successive Constraint Method (SCM): pose finding a relaxed lower bound \(\sigma_{\text{LB}}(\mu)\) as a small linear program, built from bounding/interpolating values at neighboring training points. Crucially, the size of this LP depends only on the number of affine terms \(Q\) and the number of neighboring training points used — never on \(N\).

Bounding the stability constant online: The Successive Constraint Method (SCM)

We have a fast way to compute the numerator of our error estimator (the dual norm of the residual) using the Riesz representer. But the denominator—the LBB or inf-sup constant \(\beta_\mu^\circ\)—poses a massive problem:

\[\|u_\mu - u_\mu^N\|_{\mathcal U} \le \Delta_N(\mu) = \frac{\|r_N(\mu)\|_{\mathcal V'}}{\beta_\mu^\circ}.\]

To compute the exact \(\beta_\mu^\circ\), we would have to find the smallest eigenvalue of a huge \(N \times N\) generalized eigenvalue problem. Doing this online for every new parameter \(\mu\) would completely destroy our computational efficiency, making the reduced model uselessly slow.

The Goal: We need a sharp, guaranteed lower bound \(\sigma_{\text{LB}}(\mu) \le \beta_\mu^\circ\) (or its square) that is strictly greater than zero and can be computed online in $\mathcal O(1)$ complexity independent of the high dimension \(N\).

The Brilliant Trick: Reformulation as an Optimization Problem

Thanks to Lemma 1.2.10, we know that the (squared) inf-sup or coercivity constant can be written as an infimum over a combination of the affine parameter functions \(\vartheta_q(\mu)\) and parameter-independent quotients \(h_q(v)\).

We can define a linear objective function \(J\):

\[J(\mu; \mathbf{z}) := \sum\limits_{q=1}^Q \vartheta_q(\mu) z_q,\]

where the vector \(\mathbf{z} = (z_1, \dots, z_Q)^T\) represents values coming from the high-dimensional space (i.e., \(z_q = h_q(v)\)).

The true value of our stability constant \(\sigma(\mu)\) is exactly the minimum of this function over the true set \(\mathcal S_0\) of all possible vectors \(\mathbf{z}\) generated by the infinite-dimensional space:

\[\sigma(\mu) = \inf\limits_{\mathbf{z} \in \mathcal S_0} J(\mu; \mathbf{z}).\]

The Sandwich: Upper and Lower Bounds ($\mathcal S_{\text{UB}}$ and $\mathcal S_{\text{LB}}$)

Because the true set \(\mathcal S_0\) is high-dimensional and impossible to handle online, we replace it with two surrogate sets to “sandwich” the true value.

1. The Inner Set (Upper Bound $\mathcal S_{\text{UB}}$): Offline, we pick a few training parameters \(\omega_k\) and compute the exact, true vectors \(\mathbf{z}^*(\omega_k)\) using the expensive eigenvalue solver. Our set \(\mathcal S_{\text{UB}}\) consists only of these few precomputed vectors. Since \(\mathcal S_{\text{UB}} \subset \mathcal S_0\), we are restricting our choices. A minimum over a smaller set is always greater than or equal to the true minimum, giving us our upper bound: \(\sigma_{\text{UB}}(\mu) \ge \sigma(\mu)\).

2. The Outer Set (Lower Bound $\mathcal S_{\text{LB}}$): This is the crucial part! We construct a bounding box \(\mathcal B_Q\) (by finding the absolute minimum and maximum possible values for each \(z_q\)). Then, we intersect this box with constraints: we force the points \(\mathbf{z}\) to satisfy the condition that the objective function \(J\) evaluated at our previously chosen training parameters \(\omega_k\) must be greater than or equal to the exact, known values there. Since \(\mathcal S_0 \subset \mathcal S_{\text{LB}}\), we have expanded our choices. A minimum over a larger set is always less than or equal to the true minimum, yielding our desired lower bound: \(\sigma_{\text{LB}}(\mu) \le \sigma(\mu)\).

Why is $\sigma_{\text{LB}}(\mu)$ a stroke of genius? The relaxed set \(\mathcal S_{\text{LB}}\) is nothing but a simple polytope in \(\mathbb{R}^Q\). Finding the infimum over this set is just a Linear Program (LP) of the tiny dimension \(Q\) (the number of affine terms). Solving an LP of dimension \(Q = 10\) or \(20\) online takes mere milliseconds using the Simplex algorithm!

The Offline Greedy Algorithm for SCM (Algorithm 1.4.22)

How do we decide which parameters \(\omega_k\) to evaluate offline to build our constraints? We use another Greedy algorithm!

We define a tolerance for the quality of our bounds using the relative gap (Relative Training Bound):

\[\varepsilon_K^*(\mu) := \frac{\sigma_{\text{UB}}(\mu) - \sigma_{\text{LB}}(\mu)}{\sigma_{\text{UB}}(\mu)}.\]

During the offline phase, we scan our training set for the parameter \(\mu\) where this “gap” between the upper and lower bound is the largest. We then fire up the expensive eigenvalue solver for this worst-case \(\mu\), find its exact \(\mathbf{z}^*\), and add it to our sets as a new constraint. We repeat this process until the lower bound is sufficiently tight (the gap is below our tolerance) everywhere in the parameter domain.

🎙️ Elevator Pitch for the Oral Exam

“To certify our reduced model, we need a computable error estimator. The denominator of this estimator is the inf-sup stability constant \(\beta_\mu^\circ\). Computing this exactly requires solving a high-dimensional eigenvalue problem, which would completely destroy our online efficiency.

The Successive Constraint Method (SCM) solves this by reformulating the computation as an affine optimization problem. Instead of minimizing over the intractable, high-dimensional space online, we relax the feasible set into a simple polytope \(\mathcal S_{\text{LB}}\) in the tiny space \(\mathbb{R}^Q\) (where \(Q\) is the number of affine terms). Minimizing over this enlarged set guarantees a strict lower bound \(\sigma_{\text{LB}}(\mu) \le \beta_\mu^\circ\), which can be computed online as a tiny Linear Program in \(\mathcal O(Q)\) operations.

The constraints that shape this polytope are built incrementally during the offline phase via a Greedy algorithm, which always selects the parameter point where the gap between a cheaply computable upper bound and our lower bound is maximally large.”

The weak greedy algorithm

With a fast online solve \(u_\mu^N\) and an equally fast, certified error bound \(\Delta_N(\mu)\) in hand, training the reduced space becomes a search problem: scan the parameter domain \(\mathcal P\) for the point where the current model performs worst, add a snapshot there, and repeat.

It’s called weak greedy because step 5 of the training algorithm doesn’t maximize the true (unknown, expensive) error \(\|u_\mu^\delta - u_\mu^N\|\) — it maximizes the computable estimator \(\Delta_N(\mu)\) instead. This closes the loop back to everything built above: if the error estimator is loose or unreliable, the greedy algorithm will dutifully hunt for the wrong point, and the resulting reduced space won’t actually be trained where it matters.

How good can a reduced model get? The Kolmogorov N-width

We have seen how the weak greedy algorithm constructs a reduced space \(\mathcal U_N\) by iteratively picking worst-case parameter points. But this raises a fundamental theoretical question: What is the absolute best error we could ever hope to achieve with any \(N\)-dimensional linear space?

To answer this, we need an absolute benchmark from approximation theory.

Definition 1.5.1 — Kolmogorov N-width

Let \(\mathcal P \subset \mathbb{R}^P\) be the parameter domain and \(\mathcal U\) the high-dimensional solution space. The Kolmogorov N-width of the solution manifold \(\mathcal U_{\mathcal P} := \{ u_\mu : \mu \in \mathcal P \}\) is defined as:

\(d_N(\mathcal P) := \inf\limits_{\substack{\mathcal U_N \subset \mathcal U \\ \dim(\mathcal U_N) = N}} \sup\limits_{\mu \in \mathcal P} \inf\limits_{w_N \in \mathcal U_N} \|u_\mu - w_N\|_{\mathcal U}.\)

Reading this formula from the inside out reveals its elegant logical structure:

1. \(\inf\limits_{w_N \in \mathcal U_N} \|u_\mu - w_N\|_{\mathcal U}\): For a fixed parameter \(\mu\) and a fixed subspace \(\mathcal U_N\), this is simply the best approximation error (the orthogonal projection of the true solution \(u_\mu\) onto your space).
2. \(\sup\limits_{\mu \in \mathcal P}\): We want a model that works reliably everywhere. So we search the entire parameter domain for the parameter \(\mu\) that yields the worst-case error for that given space \(\mathcal U_N\).
3. \(\inf\limits_{\substack{\mathcal U_N \subset \mathcal U \\ \dim(\mathcal U_N) = N}}\): Here is the “oracle.” We optimize over all possible \(N\)-dimensional linear subspaces in the universe. Which of these infinitely many spaces pushes this worst-case error down the furthest?

In practice, finding this optimal space is NP-hard. The Kolmogorov \(N\)-width is strictly a theoretical benchmark—the absolute smallest worst-case error any linear space of dimension \(N\) could possibly achieve.

Exponential convergence and the curse of dimensionality

This next theorem is the mathematical engine that explains why RBM is so extraordinarily successful for certain classes of equations.

Theorem 1.5.2 — The Exponential Miracle

Suppose the bilinear form \(b_\mu(\cdot,\cdot)\) is inf-sup stable (well-posed) and depends affinely on the parameter \(\mu\) with \(Q_b\) affine terms. Then there exist positive constants \(C, c > 0\) such that the Kolmogorov \(N\)-width satisfies:

\(d_N(\mathcal P) \le C e^{-c\, N^{1/Q_b}}.\)

What this means in practice:

• Exponential decay: The best possible error drops exponentially with the basis size \(N\). This is why \(N = 20\) or \(N = 30\) basis vectors are often enough to reach accuracies of \(10^{-6}\). If the error decayed only algebraically (like \(1/N\)), we would need millions of vectors.
• The hidden curse of dimensionality (\(Q_b\)): Look closely at the exponent \(N^{1/Q_b}\). If you have a simple parameterization (\(Q_b = 2\)), the rate is strongly exponential. But if you have a highly complex geometry approximated via EIM resulting in \(Q_b = 100\) affine terms, then \(1/100 \approx 0\). The term \(N^{1/100}\) barely grows as \(N\) increases, meaning the exponential decay is virtually flattened out. The more affine terms EIM needs, the worse the theoretical convergence rate of the reduced model becomes.

The Proof Idea

Why does the error fall exponentially? Where does this property come from?

The answer lies entirely in the regularity of the parameter-to-solution map: \(\mu \mapsto u_\mu\).

Because the underlying problem is inf-sup stable and affine, this mapping is incredibly “smooth”. Theorem 1.5.3 (Complex Implicit Function Theorem) uses functional analysis to show that this mapping can actually be extended into the complex numbers, where it is holomorphic (complex differentiable / analytic).

Whenever a function is holomorphic, it can be perfectly approximated by a Taylor series (or polynomials). And from approximation theory, we know that approximating analytic functions with polynomials yields exponential convergence. That is exactly where the \(e^{-c N^{1/Q_b}}\) comes from!

💡 Deeper Dive: Analytic functions vs. Wave equations

Since analyticity is the core requirement for this exponential miracle, what exactly is an analytic function, and what happens when we lose it?

1. The Mathematical Definition A function \(f(x)\) is analytic (class \(C^\omega\)) at a point if it is infinitely differentiable and can be exactly represented by its Taylor series in a neighborhood around that point:

\[f(x) = \sum_{n=0}^\infty \frac{f^{(n)}(x_0)}{n!}(x - x_0)^n.\]

2. Why is this more than just being infinitely differentiable (\(C^\infty\))? Imagine a physical wave approaching you. Before the wave arrives (e.g., \(x < 0\)), the water is perfectly flat (\(f(x) = 0\)). At \(x = 0\), the water smoothly begins to rise, for instance following the function \(f(x) = e^{-1/x^2}\) for \(x > 0\).

This function is infinitely differentiable (\(C^\infty\)) everywhere—there is no kink! But it is not analytic. If you construct the Taylor series in the flat water at \(x = -1\), all derivatives are zero. The Taylor series predicts “the water remains flat forever.” It completely fails to foresee the incoming wave. Analytic functions, by contrast, are absolutely “rigid” and predictable; knowing them on a tiny interval determines their behavior everywhere.

(Note: In the complex plane \(\mathbb{C}\), if a function is differentiable just once, it is automatically infinitely differentiable and analytic. This is called being holomorphic).

3. The Connection to MOR and the Wave Equation This explains the fundamental struggle of Model Order Reduction with transport and wave equations:

• For heat conduction or elliptic problems, the map \(\mu \mapsto u_\mu\) is holomorphic. Changing the parameter changes the solution in a perfectly predictable, “Taylor-series-like” way. Our reduced basis behaves like a polynomial expansion, capturing the solution exponentially fast.
• But real transport problems (a shockwave, a moving front, a sudden switch) are perfectly localized waves. Like our water example, they have a region of “nothing” followed by “something,” which cannot be analytic. Because analyticity breaks down, the parameter-to-solution map is no longer holomorphic, the Taylor-series proof dies, and the Kolmogorov \(N\)-width crashes from exponential decay to a painfully slow algebraic rate of \(N^{-1/2}\).

Closing the gap: Theory vs. Practice (Weak-Greedy Convergence)

In the previous section, we established the Kolmogorov \(N\)-width as the theoretical limit for model reduction. But since finding that ideal space is NP-hard, we use the Weak-Greedy algorithm to construct our space iteratively.

This sets up a duel between two protagonists:

1. \(d_N(\mathcal P)\) (The Oracle): The Kolmogorov \(N\)-width. The absolute smallest worst-case error that any space of dimension \(N\) could possibly achieve.
2. \(\sigma_N(\mathcal P)\) (The Reality): The Greedy Error. The actual worst-case error of the specific reduced space \(\mathcal U_N\) that our Weak-Greedy algorithm has cobbled together after \(N\) steps.

The ultimate goal of the convergence theorems is to bound our practical reality \(\sigma_N\) from above by the theoretical limit \(d_N\).

Theorems 1.6.1 & 1.6.2 — Weak-Greedy Convergence

The Algebraic Case (e.g., Wave/Transport): If the Kolmogorov \(N\)-width decays algebraically (\(d_N \le M N^{-\alpha}\)), then the error of our Weak-Greedy model decays algebraically at the exact same rate:

\[\sigma_N(\mathcal P) \le C M N^{-\alpha}.\]

The Exponential Case (e.g., Heat/Elliptic): If the Kolmogorov \(N\)-width decays exponentially (\(d_N \le M e^{-a N^\alpha}\)), then the error of our Weak-Greedy model also decays exponentially:

\(\sigma_N(\mathcal P) \le C M e^{-c N^\beta}.\)

💡 The Exam Detail: Why $\beta$ instead of $\alpha$?

Look closely at the exponent for the exponential case. It is no longer \(\alpha\), but rather \(\beta = \frac{\alpha}{\alpha+1}\). It is slightly degraded. Because the Greedy algorithm does not magically foresee the absolute global optimum, but rather pieces the space together step-by-step (“myopically”), it pays a tiny mathematical penalty in the exponent’s power. Despite this slight loss, it remains an exponential decay, which is a phenomenal result!

The takeaway:

• If the physics are “bad” (wave/transport) and the optimal rate is slow, the Greedy algorithm builds a space that is slow, but no worse than the limit.
• If the physics are “good” and analytic (heat/elliptic), the Greedy algorithm is smart enough to find that exponential structure!

The Engine of the Proof: Why does it work?

Why is it called the Weak Greedy algorithm? In a Strong Greedy approach, we would maximize the true error in every step. But we don’t know the true error, so we maximize our computable error estimator \(\Delta_N(\mu)\) instead.

For this to work, the estimator must perfectly “sandwich” the true projection error:

\[\gamma \Delta_N(\mu) \le \|u_\mu^\delta - P_N u_\mu^\delta\|_{\mathcal U} \le \Delta_N(\mu).\]

Here, the constant \(\gamma := \left(\frac{\beta^\circ}{C_\mu}\right)^2 \in (0,1]\). Notice what just happened: the LBB stability constant (\(\beta^\circ\)) and the continuity constant (\(C_\mu\)) have reappeared! The condition of your system matrices literally determines \(\gamma\), which in turn drives the convergence.

Lemma 1.6.3(i) — The Heartbeat of the Proof

\(\gamma \sigma_{n-1} \le |a_{n,n}| \le \sigma_{n-1}.\)

This is arguably the most important formula of the chapter. It translates the mathematical sandwich into a strict guarantee of progress:

• $\lvert a_{n,n} \rvert$: The “new” portion of the solution that the $n$-th basis vector covers.
• $\sigma_{n-1}$: The worst-case error that was still left over from the previous steps.

The formula states: Every single step of the Greedy algorithm is guaranteed to vacuum up at least a \(\gamma\)-fraction of the remaining worst-case error. Because \(\gamma > 0\) (thanks to LBB stability), the error has no choice but to shrink continuously!

The Technical Bridge (Lemma 1.6.4)

How do we turn this “step-by-step shrinking” into a proof about the global convergence rate (Theorems 1.6.1 and 1.6.2)?

That is precisely what Lemma 1.6.4 does. It is a highly technical matrix lemma, but conceptually it says: If the error has not become significantly smaller (by a factor of \(\theta\)) after another \(q m\) Greedy steps, then the current Greedy error \(\sigma_n\) is bounded from above by the Kolmogorov \(N\)-width of dimension \(m\).

The convergence theorems are then proven by contradiction. You assume the Greedy error falls slower than the \(N\)-width, apply Lemma 1.6.4, and arrive at a contradiction with the fact that \(\gamma > 0\) guarantees us relentless forward progress.

🎙️ Elevator Pitch for the Oral Exam

“Chapter 1.6 closes the gap between theory and practice. In Chapter 1.5, we saw that the Kolmogorov \(N\)-width sets the theoretical limit for model reduction. However, since finding that ideal space is NP-hard, we construct our RB space practically using the Weak-Greedy algorithm.

Theorems 1.6.1 and 1.6.2 prove that this heuristic is asymptotically optimal: if the \(N\)-width falls algebraically or exponentially, the Greedy error does too. The mathematical engine driving this proof is our error estimator. As long as the estimator bounds the true error from below by a constant \(\gamma > 0\) — which is guaranteed by the LBB stability of our Petrov-Galerkin method — Lemma 1.6.3 ensures that every Greedy step eliminates a fixed fraction of the remaining error. This relentless progress ultimately forces the Greedy algorithm to match the asymptotic behavior of the Kolmogorov \(N\)-width.”

From Worst-Case to Average: The Proper Orthogonal Decomposition (POD)

In Definition 1.5.1 of the Kolmogorov \(N\)-width, we considered the worst-case error over the parameter domain \(\mathcal P\), meaning we measured the error in \(L_\infty(\mathcal P)\). This is exactly what the Greedy algorithm targets.

But what if we are interested in the mean-squared error, i.e., in \(L_2(\mathcal P)\)? If we have a fixed training set of parameters \(\mathcal P_{\text{train}} = \{\mu_1^{\text{train}}, \dots, \mu_{N_{\text{train}}}^{\text{train}}\}\) and have already computed their corresponding snapshots, the Proper Orthogonal Decomposition (POD) yields the absolute optimal linear subspace of dimension \(N\) to represent this specific data.

Let our snapshot space be \(X^{N, \text{train}} := \text{span}\{u^{\mathcal N}(\mu_i^{\text{train}}) : 1 \le i \le N_{\text{train}}\}\).

Definition 1.7.1 — The POD Space

The POD space \(X_N^{\text{POD}}\) of dimension \(N\) is the subspace that minimizes the mean-squared projection error of the snapshots:

\[X_N^{\text{POD}} := \arg\inf\limits_{\substack{X_N \subset X^{N, \text{train}} \\ \dim X_N = N}} \frac{1}{N_{\text{train}}} \sum\limits_{\mu \in \mathcal P_{\text{train}}} \sigma_X(u^{\mathcal N}(\mu); X_N)^2,\]

where \(\sigma_X(u; Y) := \inf\limits_{v \in Y} \|u - v\|_X\) is the best approximation error.

It turns out that finding this optimal space does not require a complex search. It boils down to solving a standard eigenvalue problem. We define the correlation matrix \(\boldsymbol C^{\text{POD}} \in \mathbb{R}^{N_{\text{train}} \times N_{\text{train}}}\) with entries:

\[C_{i,j}^{\text{POD}} := \frac{1}{N_{\text{train}}} (u^{\mathcal N}(\mu_i^{\text{train}}), u^{\mathcal N}(\mu_j^{\text{train}}))_X.\]

Since this matrix is symmetric and positive (semi-)definite, it has positive eigenvalues \(\lambda_1^{\text{POD}} \ge \dots \ge \lambda_{N_{\text{train}}}^{\text{POD}} > 0\) and corresponding orthonormal eigenvectors \(\boldsymbol e_k^{\text{POD}}\).

But how do we get from these small eigenvectors back to our high-dimensional physical space?

Lemma 1.7.2 — Constructing the POD Basis

The orthonormal basis functions \(\xi_k^{\text{POD}}\) spanning the optimal POD space are given by:

\(\xi_k^{\text{POD}} := \frac{1}{\sqrt{\lambda_k^{\text{POD}}}} \sum\limits_{m=1}^{N_{\text{train}}} (\mathbf{e}_k^{\text{POD}})_m \, u^{\mathcal N}(\mu_m^{\text{train}}).\)

💡 Decoding Lemma 1.7.2: The optimal "recipe" for snapshots

Formulas with multiple indices can be intimidating, but this one is structurally beautiful. Let us completely decode it:

1. Who is who in the formula?
   • \(u^{\mathcal N}(\mu_m^{\text{train}})\): These are your raw, unmodified original snapshots (e.g., the spatial PDE solution for the \(m\)-th parameter). They are gigantic vectors of dimension \(\mathcal N\).
   • \(\mathbf{e}_k^{\text{POD}}\): This is the \(k\)-th eigenvector you computed from the correlation matrix. Since the correlation matrix is of size \(N_{\text{train}} \times N_{\text{train}}\), this eigenvector is a tiny vector of length \(N_{\text{train}}\)!
   • \((\mathbf{e}_k^{\text{POD}})_m\): This is simply the \(m\)-th entry of that \(k\)-th eigenvector. It is just a real number—it is a weight!
   • \(\frac{1}{\sqrt{\lambda_k^{\text{POD}}}}\): This is just a scaling factor. It ensures that the final POD vector \(\xi_k^{\text{POD}}\) has a length (norm) of exactly 1.

2. A Concrete Example (3 Snapshots) Imagine you only computed \(N_{\text{train}} = 3\) parameter points. You have three huge spatial snapshots: \(S_1, S_2, S_3\). You build your tiny \(3 \times 3\) correlation matrix and calculate its largest eigenvalue \(\lambda_1\). The corresponding first eigenvector \(\mathbf{e}_1\) pops out of your solver. Let’s say it is:

   \[\mathbf{e}_1 = \begin{pmatrix} 0.5 \\ 0.8 \\ -0.2 \end{pmatrix}\]

   How do you build your very first, most important spatial POD basis vector \(\xi_1^{\text{POD}}\)? You just take these three numbers and use them to weigh your snapshots:

   \[\xi_1^{\text{POD}} = \frac{1}{\sqrt{\lambda_1}} \Big( 0.5 \cdot S_1 + 0.8 \cdot S_2 - 0.2 \cdot S_3 \Big).\]

3. What does this tell us physically/geometrically? The eigenvector has “discovered” that Snapshot 2 (\(S_2\)) obviously contains the most dominant pattern for this mode (highest weight \(0.8\)). Snapshot 1 is also quite important (\(0.5\)), and Snapshot 3 is subtracted slightly (\(-0.2\)), perhaps to correct for a specific asymmetry or noise in the data.

   The mathematical elegance: Instead of searching blindly for completely new spatial functions, the mathematics prove that the optimal vectors already live in the span of your training data. The eigenvectors \(\mathbf{e}_k^{\text{POD}}\) act as the perfect “blueprints” or “recipes.” They tell you exactly how to mix your existing snapshots together to capture the maximum amount of energy (and therefore leave the smallest possible error).

🎙️ Elevator Pitch for the Oral Exam

“When we solve the eigenvalue problem for the \(N_{\text{train}} \times N_{\text{train}}\) correlation matrix, we obtain the eigenvectors \(\mathbf{e}_k\). These eigenvectors directly contain the optimal linear combination weights. The \(m\)-th entry of the \(k\)-th eigenvector is exactly the weight by which the \(m\)-th snapshot must be multiplied to form the \(k\)-th spatial POD mode. By dividing by the square root of the eigenvalue, we normalize the mode. So, instead of being arbitrary mathematical constructs, the eigenvectors of the correlation matrix are simply the perfect ‘mixing ratios’ for our snapshots to minimize the mean-squared error.”

Every PDE so far in this series has been stationary — solved once, at each parameter \(\mu\), with no notion of time. This post finally lets time in, using the classic parabolic heat equation as the running example, and compares two entirely different philosophies for handling it: chopping time into steps, versus refusing to treat time specially at all.

Warm-up: the thermal block is still the easy case

Before adding time, it’s worth remembering why the stationary thermal block from the first post was so friendly. Its bilinear form,

\[b(u,v;\mu) = \sum_{p=1}^P \mu_p \int_{\Omega_p} \nabla u(x) \cdot \nabla v(x)\, dx, \qquad \mu_p \ge \mu_0 > 0,\]

is strictly coercive, so Lax–Milgram applies directly, and standard Galerkin (\(\mathcal U^\delta = \mathcal V^\delta\)) is automatically LBB-stable with \(\beta_\mu^\circ = \mu_0\) — no supremizers needed. Keep this baseline in mind; it’s exactly the case that breaks once time enters the picture in the way we’re about to see.

The classical route: time-stepping

The parabolic heat equation reads \(u_t + A(\mu)u = f(t;\mu)\). The traditional numerical approach discretizes time first — slice \(I = (0,T)\) into steps of size \(\Delta t\) — and solves a spatial problem at every single step:

\[\mathcal L_{\Delta t}^I u_h^{k+1} = \mathcal L_{\Delta t}^E u_h^k + b_{\Delta t}^k,\]

a general time-stepping scheme that, for implicit Euler, becomes the familiar linear system

\[(\boldsymbol M_h + \Delta t\, \boldsymbol A_h)\, \mathbf u_h^{k+1} = \Delta t\, \mathbf b_h^{k+1} + \boldsymbol M_h\, \mathbf u_h^k,\]

with \(\boldsymbol M_h\) the spatial mass matrix and \(\boldsymbol A_h\) the spatial stiffness matrix.

The trap: an error estimator that compounds over time

Theorem 1.8.2 — the time-stepping error estimator \(\|u_h^k - u_N^k\|_V \;\le\; \sum_{i=1}^{k} \left[\frac{\gamma_{\Delta t}^E}{\alpha_{\Delta t}^I}\right]^{k-i} (\dots)\)

The error at step \(k\) depends on the error at step \(k-1\) (proved by induction on the time steps), and at each step the old error gets multiplied by the ratio of the continuity constant \(\gamma\) to the coercivity constant \(\alpha\). Since \(\gamma > \alpha\), this ratio is strictly greater than \(1\) — and raising it to the power \(k-i\) means the error bound grows exponentially over time. Watch this carefully: it’s the numerical stability of the scheme itself that’s fine, but the certified error bound becomes uselessly pessimistic over long time horizons, regardless of how accurate the actual computed solution is.

Algorithm 1.8.5 — POD-greedy (Haasdonk & Ohlberger)

1. Greedy over the parameter space: pick the \(\mu\) where the current error estimator is largest.
2. Compute the entire time trajectory (all time steps) for that \(\mu\).
3. Run a POD (Proper Orthogonal Decomposition) across the time steps of that one trajectory, extracting the dominant spatial mode to add to the reduced basis.

The name gives away the split: greedy handles the parameter direction, POD handles the time direction, within each single trajectory.

A cleaner picture: what a Bochner space actually is

Before going further, it helps to change how you look at \(u(t,x)\) entirely. Space and time play very different roles in the heat equation — one derivative in time (\(u_t\)), two in space (\(\Delta u\)) — so treating them as two interchangeable coordinates of one big function obscures more than it reveals.

Think of it like the difference between a photo and a film. A function in \(L^2(\Omega)\) assigns each spatial point \(x\) a value — a single photograph. A function in the Bochner space \(L_2(I;V)\) instead assigns each time \(t\) an entire photograph — a whole function \(u(t) \in V\). It’s a film: a sequence of complete spatial states, one per instant. This lets the two regularity requirements be stated separately and cleanly: at each fixed time, \(u(t)\) must be spatially well-behaved (living in the Sobolev space \(V\)); across time, the film only needs to be square-integrable.

(Note on notation: in this section \(V\) denotes the spatial Sobolev space, e.g. \(H^1_0(\Omega)\) — not the “test space” role \(\mathcal V\) played in the well-posedness posts. Same letter, new job, since it now refers to the fixed spatial space every time-slice lives in.)

Spacetime view — u(t,x) as one continuous surface. Time and space are just two coordinates of the same 2D domain. The highlighted line is the current time slice.

Bochner view — the same solution, read as a trajectory u(t) ∈ V. Each glowing line is a single element of the Hilbert space V; the whole stack, accumulating as time moves forward, is the trajectory living in the Bochner space L²(I; V).

t = 0.00

Drag the slider (or hit play) and watch the same underlying solution from two angles: on the left, one continuous spacetime surface; on the right, that same solution read as a trajectory \(u \in L_2(I;V)\) — a path traced through the Hilbert space \(V\), one snapshot at a time, which is exactly the object the POD-greedy algorithm above works with.

The space-time formulation: stop treating time specially

Given the exponentially-growing error bound above, Prof. Urban’s preferred alternative is to abandon time-stepping altogether: treat time as just another spatial dimension, turning the domain \(\Omega\) into a space-time cylinder \(I \times \Omega\), and solve the entire trajectory in one shot as a single variational problem.

The natural solution space bakes the initial condition directly in:

\[\mathcal U := H^1_{(0)}(I; V') \cap L_2(I;V),\]

where the subscript \((0)\) means every \(u \in \mathcal U\) automatically satisfies \(u(0) = 0\) — the initial condition is part of the space itself, not a side constraint.

Where the bilinear form comes from. This is really the same move you’ve seen throughout this whole series: multiply the strong form by a test function \(v(t,x)\) and integrate over everything — both \(\Omega\) and the time interval \(I\):

\[\int_0^T \int_\Omega (u_t + Au)\cdot v \; dx\, dt = \int_0^T \int_\Omega f \cdot v \; dx\, dt.\]

Splitting the integral gives the space-time bilinear form directly:

\[b(u,v;\mu) := \int_I \langle u_t, v\rangle\, dt + \int_I a_\mu(u,v)\, dt,\]

where \(\langle \cdot,\cdot\rangle\) is the \(V' \times V\) duality pairing at each fixed time (since \(u_t(t) \in V'\) in general), and \(a_\mu\) is the familiar spatial bilinear form (e.g. the coercive diffusion form from the first post).

Proposition 1.8.7: the inf-sup constant is exactly 1

For the standard heat equation (\(A = -\Delta\)), this space-time formulation has a strikingly clean stability result.

Proposition 1.8.7 For the space-time formulation of the standard heat equation, the inf-sup constant is exactly \(\beta_\mu = 1\).

🤔 Where does the supremizer $$s_u = A^{-1}u_t + u$$ come from?

This is a genuine “reverse engineering” trick, not a formula that falls from the sky. The goal is to find a test function \(v\) making \(b(u,v) = \|u\|_{\mathcal U}^2\) exactly, so start from how the norm on \(\mathcal U\) is defined:

\[\|u\|_{\mathcal U}^2 = \|u_t\|_{V'}^2 + \|u\|_{L_2(I;V)}^2 + \|u(T)\|_H^2.\]

Now look at the bilinear form: \(b(u,v) = \langle u_t, v\rangle + a(u,v)\). The first term involves \(u_t\), which lives in the “rough” dual space \(V'\) — but \(v\) needs to come from the “smooth” space \(V\) to be a legal test function. The operator \(A\) (e.g. \(-\Delta\)) takes smooth functions in \(V\) and makes them rough, landing in \(V'\); its inverse \(A^{-1}\) does the opposite — it takes something rough and translates it back into something smooth. So building \(v\) out of \(A^{-1}u_t\) is precisely a way to “translate” the rough time derivative back into a legal test function. That’s exactly why the construction is

\[s_u := A^{-1}u_t + u.\]

Proof — $$\beta_\mu = 1$$ by explicit construction

Substitute \(s_u = A^{-1}u_t + u\) into the bilinear form and expand, term by term:

\[b(u, s_u) = \big\langle u_t,\, A^{-1}u_t + u\big\rangle + a\big(u,\, A^{-1}u_t + u\big).\]

• \(\langle u_t, A^{-1}u_t\rangle = \|u_t\|_{V'}^2\), exactly the first piece of the norm.
• \(a(u,u) = \|u\|_{L_2(I;V)}^2\), exactly the second piece.
• The two cross terms, \(\langle u_t, u\rangle\) and \(a(u, A^{-1}u_t)\), combine into \(2\langle u_t, u\rangle\) once integrated over time (using that \(A\) and \(A^{-1}\) are literally inverse to each other).

Lemma 1.8.6 — chain rule in time \(\int_I 2\langle u_t, u\rangle\, dt = \|u(T)\|_H^2 - \|u(0)\|_H^2.\)

Since \(\mathcal U\) enforces \(u(0) = 0\), this cross-term integral collapses to exactly \(\|u(T)\|_H^2\) — the third and final piece of the norm. Summing all three pieces:

\[b(u, s_u) = \|u_t\|_{V'}^2 + \|u\|_{L_2(I;V)}^2 + \|u(T)\|_H^2 = \|u\|_{\mathcal U}^2.\]

Carrying out the same expansion for \(\|s_u\|_{\mathcal V}^2\) produces algebraically the exact same sum, so \(\|s_u\|_{\mathcal V} = \|u\|_{\mathcal U}\). Then

\[\beta_\mu \;\ge\; \frac{b(u,s_u)}{\|u\|_{\mathcal U}\|s_u\|_{\mathcal V}} = \frac{\|u\|_{\mathcal U}^2}{\|u\|_{\mathcal U} \cdot \|u\|_{\mathcal U}} = 1,\]

and combined with the general bound \(\beta_\mu \le 1\), this gives \(\beta_\mu = 1\) exactly.

The problem is, in this precise sense, perfectly stable — no exponentially growing estimator anywhere in sight.

Proposition 1.8.8 — the general coercive case For a general coercive spatial operator \(A_\mu\) (not just the Laplacian), the same construction \(v_w := A_\mu^{-1}\dot w + w\) still works, but \(\beta_\mu\) is no longer exactly \(1\) — it is strictly bounded below by a constant tied to the coercivity constant of \(A_\mu\).

Discretizing the cylinder: tensor products and Kronecker products

Discretizing space-time means building tensor-product spaces, \(\mathcal U^\delta = S_{\Delta t} \otimes V_h\) — piecewise-linear-in-time, finite-element-in-space. The resulting matrix uses the Kronecker product (\(\otimes\)):

\[\boldsymbol B_\mu^\delta = \boldsymbol N_{\Delta t}^{\text{time}} \otimes \boldsymbol M_h^{\text{space}} + \boldsymbol M_{\Delta t}^{\text{time}} \otimes \boldsymbol A_{h;\mu}^{\text{space}},\]

where \(\boldsymbol N_{\Delta t}\) is the time-derivative matrix and \(\boldsymbol M_h\) is the spatial mass matrix.

Multiply this enormous Kronecker system out and reorder it carefully, and something striking happens: you land, algebraically, on exactly the Crank–Nicolson scheme. So why go through all this extra machinery for a result you could have written down directly? Because you now have it as one single all-at-once system, whose LBB stability you’ve rigorously proven via Proposition 1.8.7/1.8.8 — replacing the exponentially-exploding, time-stepping error estimator of Theorem 1.8.2 with one clean, global, provably stable bound.

🔄 Question

Question: You said the space-time system is huge. How big exactly, and how do you solve it without running out of memory?

Answer: With \(N_h\) spatial degrees of freedom and \(N_t\) time steps, the all-at-once system has dimension \(N_h \cdot N_t\) — far beyond what a normal direct solver or dense storage scheme can handle. The rescue is exploiting the Kronecker structure directly (viewing the problem as a Sylvester equation) combined with low-rank approximation methods, such as the Rational Krylov Subspace Method (RKSM): assume the solution is (approximately) a low-rank matrix, and compress space and time on the fly rather than ever forming the full system.

Two philosophies, one problem

Classical time-stepping combined with POD-greedy is simple and familiar, but its certified error bound can grow exponentially over long time horizons, making the certification almost useless exactly when you’d want it most. The space-time formulation trades a conceptually bigger, unified system for a rigorously provable, time-independent stability constant — buying back the certification that time-stepping quietly loses. Which one you reach for in practice depends on whether that certification, over the whole horizon, is worth the size of the system you have to solve for it.

---

This page builds on the lecture Model Reduction by Prof. Urban, Summer Semester 2026.