Divide by Four and Conquer
Every stationary solver in this suite eventually hits the same wall: smooth error, invisible to local updates, contracting at $1-O(h^2)$. Multigrid walks around the wall — smooth the rough half, ship the smooth half to a grid where it is rough, recurse. This page walks the walk, one arrow-key at a time, and ends on the flattest row of numbers in the suite.
One problem, five grids, two moving parts
The problem is the suite's 2-D workhorse: $-\nabla\!\cdot\!\nabla u = f$ on the unit square with Dirichlet walls, discretized by poisson.poisson_2d at $n = 31$ interior points per side ($h = \tfrac1{32}$, $N = 961$ unknowns, the 5-point stencil: $4/h^2$ on the diagonal, $-1/h^2$ to each of the four neighbors). The right-hand side is the hot/cold-rod pair of reports 11 and 13, mapped from their $n=32$ grid to $n=31$ by relative position: hot rod $+1$ at column 4, rows 3–8; cold rod $-1$ at column 26, rows 22–27 (0-based interior indices). The map lands on the exact point-symmetric image of the originals — $b + \mathrm{rot}_{180}(b) = 0$, and therefore the solution is point-antisymmetric too, a parity fact that will quietly matter in S1. The first guess is seeded noise (standard_normal(seed 0), scaled to $\|x^\star\|$): half smooth, half rough, maximally unhelpful.
Two moving parts do all the work. The smoother is damped Jacobi with $\omega = \tfrac45$ — report 12's high-frequency-optimal damping for the 2-D 5-point stencil (Trottenberg et al., Multigrid), run $\nu_1 = \nu_2 = 2$ sweeps before and after each hand-off. The transfer pair is full weighting down (the $\tfrac1{16}[1\,2\,1;\,2\,4\,2;\,1\,2\,1]$ stencil) and bilinear interpolation up, tied by the classical duality $P = 4R^\top$ — exactly, the constant being $2^d$ in our $A = \mathrm{stencil}/h^2$ scaling: $R$'s rows sum to 1 (an average), $P$'s columns sum to 4 (a unit spread over $\sim$4 fine nodes), and that factor is precisely what makes $RAP$ a $2h$-scale operator. Coarse operators are rediscretized — poisson_2d again at $2h, 4h, \dots$ — with the Galerkin alternative $RAP$ computed alongside and quantified honestly in §12. Grids shrink $31 \to 15 \to 7 \to 3 \to 1$, all odd, standard $h \to 2h$; the $1{\times}1$ operator at the bottom is exactly $[[16]] = 4/h^2$ at $h = \tfrac12$, so the coarsest "solve" is a single division.
Everything in the stepper below is a baked field from results/multigrid.json — no number on this page was computed by the browser's arithmetic alone; the self-test badge re-derives a handful of identities per stage from the baked matrices at every load. Use ←/→ (or the buttons) to step, Home/End to jump; arrow keys always step the stage, never the sliders.
Nine moves, one V
Start honest: take the rod problem, the seeded random first guess, and do the obvious thing — damped-Jacobi sweeps, ω = 4/5. Watch both the error field and its 2-D DST spectrum: the rough corner of the spectrum is dead within a handful of sweeps (rough energy fraction 0.398 → 2.95×10⁻³² by sweep 100), while mode (1,1) still holds 0.6798 of its starting amplitude — exactly μ₁₁¹⁰⁰ with μ₁₁ = 0.996148. A hundred sweeps buy a factor of 4.84 in ‖e‖. The smoother is not weak; it is specialized.
‖e‖ = …
rough energy fraction …
Blame, assigned mode by mode. Damped Jacobi never mixes DST modes; mode (k, l) contracts by exactly μkl = 1 − ω (sin²(kπh/2) + sin²(lπh/2)) per sweep — report 12's per-mode factor 1 − ωλk/(4/h²), now in 2-D. The smoothing factor is the worst |μ| over the high-frequency quadrant only (max(k, l) ≥ 16): the continuous symbol gives sup = max(|1−ω/2|, |1−2ω|) = 3/5 exactly at ω = 4/5, attained at both θ = (π/2, 0) and (π, π) — and 4/5 is the argmin: no other damping does better. The discrete grid sits strictly inside the quadrant, so its worst rough mode, (1, 16), does slightly better: 0.598074, within 0.4% of 3/5 — and a straight-line fit to the measured rough-mode decay (sweeps 10 → 25) returns 0.598074 again.
gold corner = the smooth 15×15 block; outside it, the HF quadrant
hover for (k, l, μ)
slope fit (sweeps 10→25): 0.598074
The bet multigrid makes: after a few sweeps, what remains of the error fits on a grid with a quarter of the points. Test it directly — decimate the field to the 15×15 coarse nodes (keep every second fine node), interpolate straight back, and measure what the round trip lost. The 5-sweep error loses 5.9%; the 20-sweep error 2.7%; the raw sweep-0 error 74%. A pure rough-mode field loses 101% — more than all of it, because the round trip aliases rough into smooth. The representation error falls monotonically with smoothing: 0.7448 > 0.0588 > 0.0266 > 0.0163.
bilinearly interpolated back
rel. error …
So hand the problem down — but hand down the residual, not the iterate: the coarse grid will solve A₂h e = R r for the correction. Full weighting averages each coarse node over its 3×3 fine neighborhood with the 1/16 [1 2 1; 2 4 2; 1 2 1] stencil — the kron square of the 1-D [¼, ½, ¼] — so every row of R sums to 1 (all 225 rows, boundary-adjacent ones included) and constants pass through untouched. Fine 31×31 residual, norm 9.741 → coarse 15×15 right-hand side, norm 1.582; the baked coarse field is R r exactly, at every level of the hierarchy.
‖r‖ = …
‖R r‖ = …
rows of R sum to 1 · columns of P sum to 4
The coarse problem is the same equation — poisson_2d rediscretized at 2h — with a quarter of the unknowns: 961 → 225 → 49 → 9 → 1. So don't solve it; apply the same idea again. Smooth, restrict, recurse. Four hops down, the "grid" is one point, the operator is exactly [[16]] = 4/h² at h = 1/2, and the solve is a division: f = 8.50722×10⁻⁵ → u = f/16 = 5.31701×10⁻⁶. Residual there: exactly 0. The recursion has bottomed out in arithmetic a calculator could do.
Coming back up: bilinear-interpolate the coarse solution to the fine grid (P = 4 Rᵀ, exactly) and add it. Press play to watch the addition happen: the pre-smoothed fine error, norm 0.01258, meets the interpolated coarse correction and drops to 0.003058 — a factor of 4.1 from one addition. This is the moment the coarse grid pays the fine grid back: everything the smoother could not touch, solved at quarter price and shipped home.
‖e‖ = 0.0125765
(what came back up the V)
‖e(t)‖ = …
Interpolation is not free: bilinear P plants a faint high-frequency ripple on the fine grid (look at the spectrum — the rough quadrant lights back up: rough fraction of the error 0.270 right after the correction). But rough error is exactly what the smoother is for. Two post-sweeps cut each rough amplitude by roughly μ² ≈ 0.358 and the rough fraction collapses to 0.0289. The V-cycle's choreography is precisely this: every actor cleans up the mess the other one makes.
‖e‖ = …
rough energy fraction …
Assemble everything: pre-smooth and restrict on the way down (31 → 15 → 7 → 3), divide at 1, correct and post-smooth on the way up. Eighteen recorded stages, one V. Press play and walk it: the marker descends the levels while the error gauge ticks down from 0.0179 to 0.0023 — a one-cycle contraction of 0.1296. And the whole walk is linear in (x₀, b): vcycle(x₀, b) = x₀ + vcycle(0, b − A x₀) to 1.09×10⁻¹⁵ — one V-cycle is one application of M⁻¹, the suite's standard object, ready for CG in S9.
Now iterate the cycle — and change the mesh. The exact spectral radius of the V(2,2) error propagator E = I − M⁻¹A is 0.1872 / 0.1901 / 0.1908 at n = 15 / 31 / 63 (dense eigenvalues at 15 and 31, Arnoldi at 63): a spread of 0.0036 while N grows 18×. The measured per-cycle curves below fit 0.1660 / 0.1752 / 0.1789 over cycles 2–15 — the fit runs slightly below ρ because it averages the transient. Damped Jacobi on the same three problems: ρ = 0.98463 / 0.99615 / 0.99904, marching to 1 like 1 − O(h²). One number is a property of the problem size; the other is a property of the algorithm.
| n | levels | ρV exact | ρV fit (cycles 2–15) | ρJacobi | Jacobi sweeps to 10⁻¹⁰ | Jacobi sweeps ≡ one V-cycle | V-PCG iters | plain CG iters |
|---|---|---|---|---|---|---|---|---|
| 15 | 15→7→3→1 | 0.187201 | 0.165967 | 0.984628 | 1487 | 108 | 8 | 35 |
| 31 | 31→15→7→3→1 | 0.190089 | 0.175199 | 0.996148 | 5966 | 430 | 8 | 74 |
| 63 | 63→31→15→7→3→1 | 0.190818 | 0.178873 | 0.999036 | 23884 | 1718 | 8 | 149 |
The division of labor, run at every scale
The eigenmodes ledger priced every preconditioner by one question: which directions do you leave wrong? Multigrid's answer is structural: none — but no level is right about more than its own octave. The smoother at spacing $h$ owns the modes with wavelength near $2h$ — the 736 high-frequency modes here, contracted at $\le\tfrac35$ per sweep, mesh-independently. The 225 smooth modes it cannot touch are exactly the ones full weighting carries to the $2h$ grid without loss (S3's measured bet) — and on that grid, half of them have become high-frequency: a mode keeps its index $(k,l)$ but the grid's Nyquist has halved, so "smooth" is relative to where you stand. Recurse, and every mode of the original problem is high-frequency on exactly one level of the hierarchy — and gets killed there, at that level's $h$-independent rate. No mode is ever the slowest racer for long, which is why no $\kappa$ appears in the rate.
O(N) per cycle: the geometric series, audited in flops
Quartering the unknowns at every level makes the whole hierarchy cost a geometric series: $N(1 + \tfrac14 + \tfrac1{16} + \cdots) < \tfrac43 N$ cells to visit per cycle. The measured level sums land where the series says: $\sum_\ell N_\ell / N_{\rm fine} = $ 1.2622 / 1.2955 / 1.3137 at $n = 15/31/63$, increasing toward and bounded by $\tfrac43$. In report 11 Part C's flop convention (1 MAC = 2 flops, $A$-matvec $= 2\,\mathrm{nnz}(A)$; the same ledger update views uses), one V(2,2)-cycle costs 19,237 / 87,601 / 373,693 flops at the three sizes — i.e. 85.5 / 91.2 / 94.2 flops per unknown — a 10% total rise across an 18× range of $N$, settling as it grows (the last refinement, $31 \to 63$, moves it under 4%): the definition of $O(N)$, with the constant printed. In matvec currency: one cycle $=$ 9.03 / 9.36 / 9.54 fine-grid matvecs, and the whole-cycle overhead relative to fine-level-only work is 1.236 / 1.281 / 1.306 $< \tfrac43$ — the coarse levels are nearly free.
Put §9 and the series together and the economics close: to $10^{-10}$, the V-cycle spends 14 cycles $\times\ \sim\!9.4$ matvec-equivalents $\approx$ 130 matvecs of work at any size — while damped Jacobi spends 1487 / 5966 / 23,884 sweeps, each about one matvec. In drop-per-cycle terms, one V-cycle is worth 108 / 430 / 1718 Jacobi sweeps at the three sizes, and the exchange rate quadruples with every refinement (430/108 = 3.98, 1718/430 = 3.99 — the mirror image of Jacobi's $1-\rho \propto h^2$ gap shrinking $4\times$ per $h$-halving). There is no crossover to wait for; the V-cycle is ahead at $n=15$ and the lead compounds.
Why report 12 kept measuring 0.357
The two-grid error propagator is $E_{TG} = S^{\nu_2}\,(I - P A_{2h}^{-1} R A)\,S^{\nu_1}$: smooth, correct exactly on the coarse space, smooth again. The classical convergence argument (Trottenberg; transferred to the learned setting as report 06's Theorem 4.1 digest) needs exactly two properties — the smoothing property (S2: rough modes die at a rate independent of $h$) and the approximation property (S3: what survives smoothing is coarse-representable) — and this page has measured both. Here the bound closes into an identity: at $\nu = 1{+}1$ the measured two-grid spectral radius at $n=31$ is $\rho_{TG} = $ 0.3577 — for the Galerkin coarse operator and for the rediscretized one, equal to machine precision — and both equal $\mu_{\rm disc}^2 = 0.598074^2 = $ 0.3576924. The binding mode is the pure high-frequency (1, 16): every coarse correction leaves it untouched, so the rate is one pre-sweep times one post-sweep on the worst rough mode, and the choice of coarse operator cannot matter. That is why report 12's bilinear two-grid measured its mesh-independent house number 0.357 at $N = 1024$: same propagator, same identity.
The Galerkin-vs-rediscretization difference, quantified honestly since the stepper runs rediscretization: $RAP$ is symmetric positive definite and translation-invariant in the interior, but it is the classical 9-point stencil against rediscretization's 5-point one,
$$H^2 A_{RAP}=\begin{pmatrix}-\tfrac14&-\tfrac12&-\tfrac14\\[-2pt]-\tfrac12&3&-\tfrac12\\[-2pt]-\tfrac14&-\tfrac12&-\tfrac14\end{pmatrix} \qquad\text{vs}\qquad H^2 A_{2h}=\begin{pmatrix}0&-1&0\\[-2pt]-1&4&-1\\[-2pt]0&-1&0\end{pmatrix},$$
relative Frobenius difference 0.330, max entry ratio off by 0.25. They are nonetheless the same operator where it matters: both are $O(H^2)$ discretizations, and their Rayleigh quotients on the three lowest coarse modes agree within 2%. Rough-mode convergence identical (the identity above), smooth-mode behavior identical to 2% — which is why the stepper's rediscretized V-cycle and the Galerkin one are interchangeable here, and why the distinction only starts to bite when there is no grid to rediscretize on.
When there is no grid to coarsen
Everything above leaned on geometry: rediscretizing needs a mesh, and full weighting needs to know which nine fine nodes sit under a coarse one. The Galerkin product $R A P$ is the escape hatch — it builds the coarse operator from the matrix alone, and §12 just showed it loses nothing that matters. Algebraic multigrid pushes the same move all the way: choose coarse variables and interpolation weights from the matrix graph itself, no mesh required — the automation whose fragility report 05 §5 flagged as the price of leaving the model problem. One rung further and the transfer operators are learned: report 06's NAMG hard-wires exactly this page's information flow — relax, restrict, coarse solve, prolong, relax — and trains the restriction, its attention weights becoming report 11 §5.3's learned coarse regressors: data-dependent block averages chosen by gradient descent rather than by geometry.
And that is the statistical reading of the whole page, the one reports 11 and 15 built rung by rung: coarse-grid correction is regression on aggregated variables, the smoother is the local stencil regression, and a V-cycle is multiscale regression run at every octave — each level's regressors explain the correlation lengths the level below cannot reach, which is 11 §5's variance split recursed. Report 13 gave the same fact its geometric name: of the decoupling ladder's axes, scale is the one with no finite exact rung — so multigrid recurses it, and S9's flat row of rates is what the bottom of that recursion looks like.
Where the suite has already met the V
Physics
05 §5 promised "the road to multigrid" in one paragraph — smooth, restrict, recurse, $O(N)$ — and called it the end of the conversation for Poisson. This page is that paragraph, run and audited.
12 ran the two-grid cycle as a stationary predictor and kept measuring 0.357; §12 above closes that number as $\mu_{\rm disc}^2$. Update views' two-level ledger row is the same object one language over.
Probability
11 §5 measured block averages absorbing the long-range dependence of the field — the crudest coarse regressors — and §5.3 recursed the idea into multigrid-as-inference-on-aggregates. S3 is its approximation property, measured at $2h$.
15's prediction ladder ends at the same rung: a preconditioner is a statistical model, and the V-cycle is the multiscale one.
Synthesis
The eigenmodes ledger asked which directions to leave wrong; §10 is multigrid's answer — each direction handled where it is cheap. 13's scale axis, recursed; 14's separator hierarchy is the same tree read by a direct solver; 06's NAMG is the same V, learned.