Data Assimilation and PBDW: When Physics Alone Isn't Enough

Everything in the previous posts assumed the physical model itself was correct: given the right parameter $\mu$, the PDE tells you the truth, and the only errors left to worry about were discretization and reduction error. In practice that’s rarely the whole story. Every model rests on simplifying assumptions, and those assumptions introduce a modeling error that no amount of mesh refinement or better reduced bases can fix, because it was never a numerical error to begin with. The natural next step is to stop relying on the model alone and let real, measured data correct it.

The trap: data assimilation is not calibration

It’s easy to conflate two very different ways of using data, and an examiner is likely to probe exactly this distinction.

Calibration vs. data assimilation

Calibration uses input–output data to find the best-fit parameters of a physical model. The model’s structure never changes; only $\mu$ does. You’re still fully trusting that the model family is correct, and just searching within it.

Data assimilation goes further: it accepts that the model itself has errors (from simplifying assumptions, unmodeled effects, or anything else), and uses data to directly correct the predicted state, not just tune a parameter. The focus of this chapter, the PBDW (Parametrized-Background Data-Weak) method, sits firmly in this second camp.

The difference matters because it changes what you’re allowed to conclude from a good fit. Calibration tells you “this is the best $\mu$ within my model.” Data assimilation tells you “here is a corrected state, even if my model was structurally wrong to begin with.”

The two ingredients

Data assimilation needs two things on the table before anything else:

1. A best-knowledge (bk) model. This is the physical model you already have, typically living in a reduced background trial space built exactly as in the previous posts, \(\mathcal U_N = \mathrm{span}\{\xi_1, \dots, \xi_n\} \subset \mathcal U.\) It’s “best-knowledge” rather than “truth” precisely because you’re about to admit it might be wrong.

2. Observations. You have $M$ real measurements, modeled as linear observational functionals $\ell_m^o \in \mathcal W’$, where $\mathcal W$ is a Hilbert space large enough to contain the state space, $\mathcal U \subseteq \mathcal W$ (the “update space”). Measurements are noisy: \(\ell_m^{obs} = \ell_m^o(u_{true}) + \varepsilon_m.\) A sensor reading a temperature at one point, or a strain gauge on a structure, are typical examples of such functionals: each one extracts a single real number from the (unknown) true state.

Building up to PBDW

Step A — no parameters yet

Strip away parameters entirely and suppose you just have one fixed background solution $u_{bk}$. You want a corrected state $w$ that stays close to $u_{bk}$ while also explaining the measured data as well as possible:

\[u_\kappa^* := \arg\min_{w \in \mathcal W} \left\{ \kappa \|w - u_{bk}\|_{\mathcal W}^2 + \frac{1}{M} \sum_{m=1}^M \big(\ell_m^o(w) - \ell_m^{obs}\big)^2 \right\}.\]

This is a regularized least-squares problem, and the regularization weight $\kappa > 0$ is doing all the interesting work:

🤔 What does $\kappa$ actually control?

• As $\kappa \to 0^+$: the penalty on straying from $u_{bk}$ vanishes, and the optimizer is free to chase the data exactly. If the data is noiseless this is fine; if it’s noisy (as it almost always is), you end up exactly interpolating the noise itself, which is usually a bad idea.
• As $\kappa \to \infty$: the data term becomes negligible next to the penalty, and the optimizer collapses back onto $w = u_{bk}$. You’re trusting the physical model completely and ignoring the measurements.

$\kappa$ is therefore a dial between two failure modes: overfitting noise on one end, and ignoring real information on the other. Its role is directly analogous to a regularization parameter in any inverse problem, and choosing it well is its own (separate) topic.

Step B — bringing parameters back: the PBDW statement

Now reintroduce the parameter dependence. Instead of searching directly over one corrected state $w$, split the unknown into two pieces: a background contribution $u \in \mathcal U_N$ (cheap, reduced, physics-based) and an update $w \in \mathcal W$ (the correction, informed by data). The PBDW formulation reads

\[(u_\kappa^*, w_\kappa^*) := \arg\min_{(u,w) \in \mathcal U_N \times \mathcal W} \underbrace{\left\{ \kappa \|w\|_{\mathcal W}^2 + \frac{1}{M} \sum_{m=1}^M \big(\ell_m^o(u+w) - \ell_m^{obs}\big)^2 \right\}}_{=: J_\kappa(u,w)}.\]

In plain terms: find the best combination of a background state $u$, drawn from the cheap reduced model space, and a correction $w$, drawn from the richer update space, such that the sum $u + w$ explains the real measurements as well as possible, while keeping the size of the correction itself under control via $\kappa$.

Making the update space finite: Riesz representers

$\mathcal U_N$ is already finite-dimensional by construction, spanned by $N$ snapshots. But $\mathcal W$, the update space, is an infinite-dimensional Hilbert space, and you can’t optimize over that directly on a computer. Something has to turn it into a finite-dimensional search too.

The trick is to stop thinking of a sensor as a “measuring machine” and start thinking of it as a function. A measurement functional $\ell_m^o$ eats a function and returns a number; on its own, that’s not something you can span a space with. The Riesz representation theorem fixes this: for every bounded linear functional $\ell_m^o$ on a Hilbert space $\mathcal W$, there is a unique function $q_m \in \mathcal W$ such that

\[q_m := R_{\mathcal W}\, \ell_m^o, \qquad \ell_m^o(w) = (q_m, w)_{\mathcal W} \quad \text{for all } w \in \mathcal W.\]

$q_m$ is the sensor’s “footprint” inside the function space: applying the measurement is literally the same as taking an inner product against $q_m$. A sensor that reads a value at a single point, for instance, has a representer that’s sharply peaked exactly at that point (the same kind of kernel function that reappears later with radial basis functions). Because $q_m$ is an honest function rather than an abstract functional, it can actually be used to build a subspace:

\[\mathcal W_M := \mathrm{span}\{q_1, \dots, q_M\} \subset \mathcal W.\]

This turns the infinite-dimensional search over $\mathcal W$ into a search over an $M$-dimensional space spanned by the sensor footprints themselves, exactly the same finite basis trick used for $\mathcal U_N$, just built from data instead of snapshots.

q_m (Riesz representer)

ker(ℓ_m)

projection of x

x (test function)

a (weight on x₁) 1.20

b (weight on x₂) 0.40

x₁ (test point) 1.40

x₂ (test point) 0.90

q_m = (a, b) = (1.20, 0.40)
ℓ_m(x) = a·x₁ + b·x₂ = 2.04
(x, q_m) = 2.04

---

residual · q_m (should be ≈ 0) = 0.00

Where the optimal update actually lives: Proposition 2.2.1

Proposition 2.2.1

Every solution \((u_\kappa^*, w_\kappa^*)\) of the PBDW problem satisfies

\(u_\kappa^* \in \mathcal U_N, \qquad w_\kappa^* \in \mathcal W_M \cap (\mathcal U_N)^\perp.\)

In words: the optimal update is not just some element of the infinite-dimensional $\mathcal W$. It is built entirely out of sensor footprints, and it is orthogonal to everything the physical background model could already represent.

Proof idea — why the optimal update collapses onto $\mathcal W_M \cap (\mathcal U_N)^\perp$

The argument has two independent parts, and both work by showing that any leftover piece of $w$ outside the claimed space is dead weight: it does nothing useful and only adds cost.

Part 1 — why $w \in \mathcal W_M$. Split any candidate $w \in \mathcal W$ orthogonally with respect to the sensor-footprint space, \(w = P_{\mathcal W_M} w + P_{(\mathcal W_M)^\perp} w.\) Apply a measurement functional to this decomposition. Since \(\ell_m^o(w) = (q_m, w)_{\mathcal W}\) and the inner product is linear, \(\ell_m^o(w) = (q_m, P_{\mathcal W_M} w)_{\mathcal W} + \underbrace{(q_m, P_{(\mathcal W_M)^\perp} w)_{\mathcal W}}_{=\,0},\) where the second term vanishes because $q_m \in \mathcal W_M$ is orthogonal to anything in $(\mathcal W_M)^\perp$. So the sensors are completely blind to $P_{(\mathcal W_M)^\perp} w$: it changes nothing about how well the data is fit. But it does change the cost, since by Pythagoras \(\|w\|_{\mathcal W}^2 = \|P_{\mathcal W_M} w\|_{\mathcal W}^2 + \|P_{(\mathcal W_M)^\perp} w\|_{\mathcal W}^2\), and the regularization term \(\kappa\|w\|_{\mathcal W}^2\) grows with it. A component that helps nowhere and costs something can never belong to a minimizer, so \(P_{(\mathcal W_M)^\perp} w^* = 0\) at the optimum: the update lives entirely in \(\mathcal W_M\).

Part 2 — why $w \in (\mathcal U_N)^\perp$. Now decompose $w \in \mathcal W_M$ with respect to the background space instead, $w = P_{\mathcal U_N} w + P_{(\mathcal U_N)^\perp} w$. The combined solution becomes \(u + w = \underbrace{u + P_{\mathcal U_N} w}_{\text{both terms live in } \mathcal U_N} + P_{(\mathcal U_N)^\perp} w.\) Since $u$ is completely unpenalized in the cost functional, any part of $w$ that already lies in $\mathcal U_N$ can simply be absorbed into $u$ instead, without changing the data-fit term at all, while strictly lowering the regularization cost (because $u$ pays nothing for it, whereas $w$ pays $\kappa$). Concretely, \(J_\kappa\big(u + P_{\mathcal U_N} w,\; P_{(\mathcal U_N)^\perp} w\big) + \kappa\|P_{\mathcal U_N} w\|^2 = J_\kappa(u, w),\) and since $\kappa|P_{\mathcal U_N} w|^2 \ge 0$, the left-hand side without that extra term is never larger. So the minimum can only be attained once $P_{\mathcal U_N} w^* = 0$: the update carries no component that the background model could already have supplied for free.

Putting both parts together: the optimal update must be built purely from sensor footprints (Part 1) and must be orthogonal to whatever the reduced physical model already spans (Part 2). Its whole purpose is to supply exactly what the background model structurally cannot represent, nothing the sensors can’t see, and nothing the model could have given you anyway.

Practice question

Test yourself

Question: In the PBDW statement, why does $\kappa$ only penalize $|w|_{\mathcal W}$, the update, and not $|u|_{\mathcal U_N}$, the background state?

Answer: The background state $u$ is not a free variable in the same sense as $w$ — it’s constrained to lie in the reduced space $\mathcal U_N$, which was already built to represent physically meaningful solutions of the best-knowledge model. There’s no need to additionally penalize its size, because staying inside $\mathcal U_N$ already keeps it “close to the physics” by construction.

The update $w$ is different: it lives in the much larger, essentially unconstrained space $\mathcal W$, precisely because its whole job is to absorb whatever the background model got wrong. Without a penalty, $w$ could grow arbitrarily large and simply overwrite the entire background contribution, at which point the method would degenerate into pure data interpolation, exactly the $\kappa \to 0^+$ failure mode from Step A. Penalizing $|w|$ keeps the correction to the minimal adjustment needed to fit the data, so that the background model still does as much of the work as it can, and the data is used to correct it rather than replace it.

From infinite spaces to linear algebra: The Saddle Point Problem

Knowing that $w_\kappa^* \in \mathcal W_M$ is conceptually great, but a computer needs matrices. By representing the background state as $\mathbf{u} \in \mathbb{R}^N$ (coefficients of $\xi_n$) and the update as $\mathbf{w} \in \mathbb{R}^M$ (coefficients of the Riesz representers $q_m$), the minimization of $J_\kappa$ boils down to setting its derivatives to zero.

This yields a block-linear system known as a Saddle Point Problem (see Proposition 2.2.3):

\[\begin{bmatrix} \kappa M \mathbf{I} + \mathbf{K} & \mathbf{L} \\ \mathbf{L}^T & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{w} \\ \mathbf{u} \end{bmatrix} = \begin{bmatrix} \boldsymbol{\ell}_M^{obs} \\ \mathbf{0} \end{bmatrix}\]

The structure of these matrices is beautiful and tells a physical story:

• $\mathbf{K} \in \mathbb{R}^{M \times M}$ (Sensor overlap): $K_{m,m’} = (q_m, q_{m’})_{\mathcal{W}}$. This Gramian matrix measures how much the “footprints” of the different sensors overlap with each other.
• $\mathbf{L} \in \mathbb{R}^{M \times N}$ (Physics meets Data): $L_{m,n} = \ell_m^o(\xi_n)$. This matrix evaluates what the $m$-th sensor would measure if the true state was exactly the $n$-th background basis function.

Proposition 2.2.5 (The Rank Condition)

The PBDW saddle point problem admits a unique solution if and only if the matrix $\mathbf{L}$ has full column rank, i.e., $\mathrm{rank}(\mathbf{L}) = N$.

Test yourself: Physical meaning of the Rank Condition

Question: What simple inequality must hold for $M$ (number of sensors) and $N$ (number of background modes) for Proposition 2.2.5 to be satisfied, and what does this mean physically?

Answer: Since $\mathbf{L}$ is an $M \times N$ matrix, it can only have full column rank $N$ if $M \ge N$.

Physically, this means we need at least as many sensors as we have physical background modes. If $M < N$, the problem is underdetermined. There would be infinitely many ways to mix our $N$ physical modes to perfectly match a tiny number of sensor readings. The algorithm wouldn’t know which physical reality is the correct one. Having $M \ge N$ ensures we have enough data to “nail down” the coefficients of our reduced background model unambiguously.

---

Making point measurements rigorous: RKHS

Up to this point, we assumed the measurement functionals $\ell_m^o \in \mathcal W’$ were just “some” bounded linear functionals. But in reality, sensors are usually placed at specific coordinates $x_m$. We want to measure the state at a point: $\ell_m^o(w) = w(x_m)$.

The problem: In standard $L_2$ spaces, functions are only defined “almost everywhere”. Evaluating an $L_2$ function at a single, infinitely small point $x_m$ is mathematically meaningless (the Dirac delta is not a function in $L_2$). Consequently, we cannot find a standard Riesz representer for a point measurement in these spaces.

The solution: We must choose an update space $\mathcal W$ in which point evaluation is a bounded, continuous operation. This leads us to Reproducing Kernel Hilbert Spaces (RKHS).

Definition 2.3.1 (Reproducing Kernel)

A function $K: \Omega \times \Omega \to \mathbb{R}$ is called a reproducing kernel of a Hilbert space $\mathcal H$ if:

1. $K(\cdot, x) \in \mathcal H$ for all $x \in \Omega$.
2. Reproducing property: $(f, K(\cdot, x))_{\mathcal H} = f(x)$ for all $f \in \mathcal H$.

This property is pure magic for data assimilation. By using an RKHS, the abstract Riesz representer $q_m$ of a point measurement at $x_m$ simply becomes the kernel function centered at that point: $q_m = K(\cdot, x_m)$.

Instead of an infinitely sharp Dirac delta needle, the kernel acts as a smooth “bell curve” around the sensor. The reproducing property guarantees that taking the inner product of any state $w$ with this bell curve returns the exact value of $w$ at the sensor’s location.

This also makes assembling our matrices incredibly easy:

• $\mathbf{K}{m,m’} = (q_m, q{m’}){\mathcal W} = (K(\cdot, x{m’}), K(\cdot, x_m)){\mathcal W} = \mathbf{K(x_m, x{m’})}$
• $\mathbf{L}_{m,n} = \ell_m^o(\xi_n) = \mathbf{\xi_n(x_m)}$

We just plug the sensor coordinates into the kernel function and the background modes!

Building the Update Space: Wendland Functions

What kind of kernel function should we choose? We want a function that is easy to compute, smooth, and—most importantly—zero if you move far enough away from the sensor.

Definition 2.3.3 (Wendland functions / csRBF)

We typically choose compactly supported Radial Basis Functions (csRBF). A popular choice are Wendland functions $\Phi_{d,k}$, where $d$ is the spatial dimension and $k$ controls the smoothness. The kernel is defined via the Euclidean distance: \(K(x,y) := \Phi_{d,k}(\|x-y\|_2)\)

Because the Wendland functions have a compact support (they drop exactly to zero beyond a certain radius), sensors that are far apart don’t overlap. This means $(q_m, q_{m’})_{\mathcal W} = 0$. Our matrix $\mathbf{K}$ becomes a sparse matrix, which makes solving the saddle point problem computationally blazing fast, even for thousands of sensors!

Theorem 2.3.4 (The Native Space)

The Hilbert space generated by the Wendland kernel $K_{d,k}$ is a Sobolev space of specific smoothness: \(\mathcal{H}_K = H^{k + \frac{d+1}{2}}(\mathbb{R}^d)\)

Why is this theorem so crucial for exams? Remember the fundamental requirement for PBDW from Chapter 2.1: $\mathcal U_N \subset \mathcal W$. The physical background space must be a subspace of the update space. Theorem 2.3.4 gives us the tool to guarantee this! If our PDE solutions are very smooth (e.g., they live in $H^2$ or $H^3$), we simply increase the parameter $k$ of our Wendland function. This creates a smoother “Native Space” $\mathcal W$ that is guaranteed to contain our physical space $\mathcal U_N$.

Where to put the sensors? The SGreedy-plus Algorithm

If you have $M$ sensors to place in a domain $\Omega$, you shouldn’t just scatter them randomly. The SGreedy-plus algorithm (Algorithm 2.3.5) optimizes sensor placement by balancing two competing goals:

1. Maximize Stability (The “SGreedy” loop): We want the saddle point problem to be well-posed, meaning the inf-sup constant $\beta_{N,m}$ must be large. The algorithm iteratively places sensors where the background modes $\xi_n$ have their maximum variance. Intuition: Place sensors where the physics models predict the most “action”. A sensor placed where all basis functions are always zero provides zero information about the system state.

2. Minimize Fill Distance (The “plus” loop): Once stability is reached ($\beta_{N,m} > tol$), the algorithm switches to pure geometry. It places the remaining sensors to minimize the “fill distance” (the radius of the largest empty sphere in the domain without a sensor). Intuition: Do not leave massive “blind spots” in the room. Even if the physics model says nothing interesting happens in the corner, a window might be left open in reality. We need geometric coverage to let the update $w$ catch local anomalies that the background model doesn’t know about.

---

Chapter 3: Nonlinear Problems

Up to this point, all parameterized problems were linear, resulting in a bilinear form $b_\mu(u,v)$. This allowed us to rely on the Lax-Milgram lemma or the Banach-Nečas theorem to guarantee a unique solution. However, for many real-world applications—especially in fluid dynamics (like the Burgers equation or Navier-Stokes equations)—this linear framework completely collapses. Convection terms like $u \cdot \nabla u$ introduce nonlinearities. Lax-Milgram is officially dead.

3.1 Quadratically Nonlinear Problems

In fluid mechanics, the nonlinearity is almost always quadratic (a state multiplied by its own spatial derivative). We can decompose such problems into a linear part and a quadratic part.

Definition 3.1.2 (Quadratically nonlinear form)

A quadratically nonlinear form $b_\mu : \mathcal{U} \times \mathcal{V} \times \mathcal{P} \to \mathbb{R}$ is decomposed as: \(b_\mu(u, v) = a_0(u, v; \mu) + a_1(u, u, v; \mu)\)

• $a_0(u, v; \mu)$: A standard bilinear form (e.g., representing linear diffusion/viscosity: $-\mu_1 \Delta u$).
• $a_1(u, w, v; \mu)$: A trilinear form. It is linear in each of its three arguments. When we plug the state $u$ into the first two slots, it generates the quadratic nonlinearity (e.g., convection $u \cdot \nabla u$).

Note: The form $b_\mu(u,v)$ is still perfectly linear with respect to the test function $v$! The nonlinearity only applies to the trial function $u$.

3.2 Well-Posedness: Fixed-Point Theorems

Since we cannot use Lax-Milgram, we must prove well-posedness using Fixed-Point Theorems. In nonlinear PDE theory, we must strictly separate the proof of existence from the proof of uniqueness.

1. Existence (Brouwer’s Fixed Point Theorem)

To prove that at least one solution exists, we rely on the Brouwer Fixed-Point Theorem (Thm 3.2.1 / Cor 3.2.2). Because Brouwer only works in finite-dimensional spaces, the proof (Thm 3.2.3) uses a classic Galerkin argument:

1. Restrict the problem to a finite-dimensional subspace $\mathcal{U}_m$ and use Brouwer to prove a discrete solution $u_m$ exists.
2. Let the dimension $m \to \infty$ and use weak convergence in the Hilbert space to show the limit $u$ is a valid solution in the infinite-dimensional space $\mathcal{U}$. Drawback: This only guarantees that at least one solution exists. It tells us nothing about uniqueness.

2. Uniqueness (Banach’s Fixed Point Theorem)

To prove that exactly one unique solution exists, we use the Banach Fixed-Point Theorem (Contraction Mapping Theorem). We rewrite the PDE as a fixed-point iteration $u = T(u)f$ and must prove that $T(u)f$ is a contraction.

Theorem 3.2.4 (Uniqueness & Small Data Condition)

The quadratically nonlinear PDE admits exactly one unique solution if:

1. The form is uniformly elliptic in the 2nd and 3rd argument: $b(w,v,v) \ge \alpha |v|^2_{\mathcal{U}}$.
2. The mapping is locally Lipschitz continuous, bounded by a monotonically increasing function $L(\xi)$.
3. The Small Data Condition holds: \(\frac{\|f\|_{\mathcal{U}'}}{\alpha^2} L\left( \frac{\|f\|_{\mathcal{U}'}}{\alpha} \right) < 1\)

💡 Exam Gold: The Physical meaning of the Small Data Condition

Why does the Lipschitz constant $L$ depend on a radius $\xi = |f|/\alpha$? And why does mathematics suddenly require “small forces” ($f$) or “high viscosity” ($\alpha$) to guarantee a unique solution? Does this reflect reality?

Yes! The Small Data Condition is essentially the mathematical equivalent of the Reynolds number ($Re$).

• In fluid dynamics, the coercivity constant $\alpha$ corresponds to the kinematic viscosity (e.g., honey has a large $\alpha$, water a small one). The right-hand side $f$ corresponds to the driving forces.
• If forces are small and viscosity is high (small $Re$, left side of the inequality is $< 1$), the flow is laminar. Physics dictates there is exactly one stable, unique flow state. Mathematics perfectly reflects this by granting uniqueness via Banach.
• If forces are huge or viscosity is tiny (large $Re$, inequality is violated), the flow becomes turbulent or bifurcates (e.g., Kármán vortex street). Physics allows for multiple valid states (does the vortex spin left or right?). Because multiple states physically exist, the mathematical guarantee of uniqueness must break down. The theorem failing at this exact threshold is a beautiful mirror of reality!

Preparing for the Computer: The Fréchet Derivative

To actually solve this nonlinear problem numerically, we will use Newton’s method (Section 3.3). Newton’s method requires us to “take the derivative” of our PDE operator. For a nonlinear form $b_\mu(u,v)$, the directional derivative with respect to $u$ in the direction of a perturbation $z$ is the Fréchet derivative:

\[db_\mu(u,v)[z] = a_0(z,v) + a_1(z,u,v) + a_1(u,z,v)\]

Exam Intuition: This is nothing but the standard Product Rule! Just like the derivative of $f(x) = c x + d x^2$ is $f’(x) = c + 2dx$, the derivative of our form treats the linear part $a_0$ normally (replace $u$ with $z$), and applies the product rule to the quadratic part $a_1$ (differentiate the first slot, keep the second + keep the first slot, differentiate the second).

3.3 Newton’s Method for Nonlinear PDEs

Because the operator is nonlinear, we cannot invert it in a single step using Lax-Milgram. Instead, we must solve the problem iteratively using Newton’s Method in Hilbert spaces (Algorithm 3.3.2).

Just like classical Newton ($x_{k+1} = x_k - f(x_k)/f’(x_k)$), we find an update step $z^{(k)}$ by setting the linear derivative equal to the negative residual of the current guess:

1. The Residual: $-F_\mu(u^{(k)}) = f(v) - b_\mu(u^{(k)}, u^{(k)}, v)$
2. The Derivative (Fréchet Derivative): $db_\mu(u^{(k)}, \cdot)[z]$
3. The Newton Step: Find $z^{(k)}$ by solving the linear system: \(db_\mu(u^{(k)}, v; \mu)[z^{(k)}] = f(v) - b_\mu(u^{(k)}, u^{(k)}, v) \quad \forall v \in \mathcal{V}\)

The Product Rule Intuition: The Fréchet derivative of our form $b_\mu = a_0 + a_1$ in the direction of $z$ is simply the product rule applied to the nonlinear slots: \(db_\mu(u,v)[z] = a_0(z,v) + a_1(z,u,v) + a_1(u,z,v)\) If the trilinear form $a_1$ is symmetric (Def 3.3.1: $a_1(u,w,v) = a_1(w,u,v)$), this beautifully simplifies to $2a_1(u, z, v) + a_0(z, v)$.

The Stability Trap in Newton’s Method

In linear PDEs, we only have to prove LBB stability ($\beta > 0$) once. However, in Newton’s method, the linear operator we invert is the Fréchet derivative, which depends on the current state estimate $u^{(k)}$. This means the inf-sup constant changes in every single iteration step! The problem must remain well-posed along the entire solution trajectory. Mathematical frameworks like Brezzi-Rappaz-Raviart (BRR) theory are required to provide a-posteriori stability bounds for this moving target.

3.4 Model Reduction: The Tensor Problem

To make this efficient, we project the Newton iteration onto the reduced basis space $U^N$. However, the offline-online decomposition suffers from a severe dimensional explosion.

For linear problems, projecting a bilinear form $a(u,v)$ onto an $N$-dimensional space creates an $N \times N$ matrix. But our form $a_1(u,w,v)$ is trilinear (it takes three arguments). Projecting it onto an $N$-dimensional basis generates a 3rd-order tensor $\mathbf{A}_N^{1,q} \in \mathbb{R}^{N \times N \times N}$ (Equation 3.4.3).

The role of $Q^{a_1}$: Just like we did for linear problems, we must affinely decompose the trilinear form to separate parameters from spatial variables: \(a_1(u, w, v; \mu) = \sum_{q=1}^{Q^{a_1}} \vartheta_q^{a_1}(\mu) \, a_{1,q}(u, w, v)\) Here, $Q^{a_1}$ is the number of affine terms. This means in the offline phase, we don’t just store one tensor, but $Q^{a_1}$ different $N \times N \times N$ tensors!

• Storage & Online Cost: $O(Q^{a_1} \cdot N^3)$. If $N$ grows beyond $\approx 50$, the $N^3$ scaling causes the “fast” online phase to become incredibly slow and memory-intensive. (This is the primary motivation for hyper-reduction methods like EIM/DEIM, which approximate the nonlinearity to bypass the tensor entirely).

3.4.1 A Posteriori Error Estimation

To certify the reduced Newton solution $u_N^{(k)}$, we need to guarantee how far it deviates from the high-fidelity (truth) Newton iterate $u^{(k)}$.

Theorem 3.4.1 (A Posteriori Error Bound)

This theorem provides a rigorous, computable a-posteriori error bound $\Delta_N^{(k)}(\mu)$. It mathematically guarantees that at any iteration step $k$, the true error is bounded by: \(\|u^{(k)} - u_N^{(k)}\|_{\mathcal{U}} \le \Delta_N^{(k)}(\mu)\)

Unlike linear problems (where Error $\le$ Residual / $\beta$), the nonlinear error bound (Equation 3.4.4) looks intimidating because it contains a sum over a product.

The Physical Intuition (The “Ship in a Storm” Analogy): In a linear system, an error made in step 1 stays proportional. In a nonlinear system, errors have a memory and can compound exponentially. Imagine steering a ship in a turbulent, nonlinear current: a tiny steering error in step 1 puts you in a slightly wrong position for step 2, where suddenly completely different current forces act on you, magnifying the initial mistake.

The formula captures exactly this:

1. The Sum ($\sum_{i=1}^k$): It accumulates the dual norms of the residuals $r_N^{(i)}$ from all past Newton steps. The error has a memory.
2. The Product ($\prod$): This is the amplification factor. It determines how much a past error is “blown up” by the nonlinearity before reaching the current step $k$. It depends heavily on the Lipschitz constants $L_F, L_{dF}$ (how aggressive the nonlinearity is) and the stability $\beta_\mu$.

Exam Question: Why does the error estimator NOT converge to zero as $k \to \infty$?

Question: If we let Newton’s method run for infinitely many iterations ($k \to \infty$) and it converges perfectly, why does the error bound $\Delta_N^{(k)}$ in Remark 3.4.2 not go to zero?

Answer: Because Newton convergence and RB approximation are two different things! Even if the RB Newton method converges perfectly to a final state $u_N^{(\infty)}$, this state is still confined to the $N$-dimensional reduced space $U^N$. The truth solution $u^{(\infty)}$ lives in the massive space $U^\delta$. There is an inherent, unresolvable projection error (the RB error $\varepsilon$) between these two spaces.

Since the error bound $\Delta_N^{(k)}$ measures the total distance between the truth iterate and the RB iterate, it can never drop below this fundamental RB projection error, regardless of how many Newton steps you take. In fact, because the formula accumulates past errors, $\Delta_N^{(k)}$ grows monotonically with $k$, often leading to a very pessimistic (poor effectivity) bound.

This page builds on the lecture Model Reduction by Prof. Urban, Summer Semester 2026 — written in my own words rather than reproducing lecture material directly.