Multigrid Methods: Killing Error at Every Scale

The precise claim: mesh-independent, O(N) convergence

Fix a bounded domain and discretize an elliptic PDE (the model case is −Δu = f) on a nested hierarchy of grids Ω₁ ⊂ Ω₂ ⊂ ⋯ ⊂ Ω_L with spacings h₁ > h₂ > ⋯ > h_L, giving symmetric positive-definite (SPD) stiffness matrices A_k. The multigrid V-cycle is a linear stationary iteration x ↦ x + M(b − A_L x) whose error-propagation matrix E_MG = I − M·A_L satisfies

• ‖E_MG‖_{A} ≤ ρ₀ < 1, where ‖·‖_A is the energy norm and the constant ρ₀ is independent of the mesh size h and of the number of levels L.

Because each V-cycle costs O(N) arithmetic (the grid sizes form a geometric series that sums to a constant times the finest), reducing the error by any fixed factor takes O(1) cycles, so total work is O(N). Full multigrid (FMG) — nesting V-cycles from coarse to fine — sharpens this to solving to the level of discretization error in O(N) work. This optimality, uniform in h, is the whole point.

The picture: high frequencies die locally, low frequencies die on coarse grids

Expand the error in the eigenvectors of A. Classical relaxations — Jacobi, Gauss–Seidel — are excellent at damping high-frequency (oscillatory) error: a weighted-Jacobi sweep multiplies the amplitude of the mode with frequency θ by roughly |1 − ω(1−cos θ)|, which is tiny for θ near π but ≈ 1 for θ near 0. So after a few sweeps the error is smooth — that is the smoothing property.

But smooth error is exactly what relaxation cannot kill, because a smooth function looks locally like a constant and local averaging barely changes it. The trick: a function that is smooth on a fine grid is oscillatory relative to a coarse grid. Transfer the smooth residual down to a grid with twice the spacing, solve the residual equation there (recursively), and interpolate the correction back up. On the coarse grid the once-smooth error is now high-frequency and again cheap to relax. Recursing through all levels, every frequency band meets a grid on which it is oscillatory — hence 'killing error at every scale.'

The mechanism of the proof: smoothing property × approximation property

The rigorous convergence proof (Bank–Dupont, Braess–Hackbusch, Hackbusch 1985) factors the two-grid error operator as a product and bounds the factors separately. Two-grid error propagation is

E_TG = S^ν · ( I − P A_H⁻¹ R A_h ) ,

where S is the smoother, ν the number of smoothing steps, P interpolation (prolongation), R restriction (R = Pᵀ up to scaling), and A_H the coarse operator. One establishes two estimates in the energy inner product:

• Smoothing property: ‖A_h S^ν v‖ ≤ η(ν) h⁻² ‖v‖ with η(ν) → 0 as ν → ∞ — many smoothing steps flatten the residual, and the h⁻² is exactly the norm of A_h.
• Approximation property: ‖(A_h⁻¹ − P A_H⁻¹ R) v‖ ≤ C h² ‖v‖ — the coarse space approximates A_h⁻¹ to order h², from elliptic regularity and finite-element approximation theory.

The h⁻² and h² cancel, giving ‖E_TG‖ ≤ C η(ν) < 1 for ν large enough, uniformly in h. Extending from two grids to the full V-cycle uses a perturbation/energy argument (or Xu–Zikatanov's sharp subspace-correction identity, 2002) to keep the bound level-independent.

Worked example: 1D Poisson and the smoothing factor

Take −u″ = f on (0,1) with u(0)=u(1)=0, discretized by second differences on a uniform grid of spacing h = 1/n. The matrix is the tridiagonal A = h⁻²·tridiag(−1, 2, −1), with eigenvectors the sine modes v_k(j) = sin(k π j h) and eigenvalues λ_k = 4h⁻² sin²(k π h / 2), k = 1,…,n−1.

Apply weighted Jacobi x ← x + ω D⁻¹ r with D = 2h⁻² I. Mode k is multiplied each sweep by g_k = 1 − ω·2 sin²(kπh/2). Choose ω = 2/3. Then for the oscillatory half of the spectrum (n/2 ≤ k ≤ n−1) one checks |g_k| ≤ 1/3 — the smoothing factor. High modes lose at least ⅔ of their amplitude per sweep, so a few sweeps leave a genuinely smooth error.

• The smooth modes (k small) have g_k ≈ 1 − O(k²h²) ≈ 1: Jacobi cannot touch them.
• Restricting to the grid with spacing 2h, a mode k < n/2 has wavenumber that is now a large fraction of the coarse Nyquist limit — oscillatory again — so coarse relaxation crushes it.

This 1/3 smoothing factor, combined with the O(h²) coarse approximation, is precisely what yields the h-independent V-cycle rate.

Why the hypotheses matter — and what breaks

The clean O(N), h-independent result rests on structure that can fail:

• Ellipticity / SPD-ness and H² regularity. The approximation property ‖(A_h⁻¹ − P A_H⁻¹ R)‖ ≤ C h² uses elliptic regularity ‖u‖_{H²} ≤ C‖f‖_{L²}. On domains with re-entrant corners regularity drops to H^{1+s}, s<1, and standard analysis gives only a suboptimal, though still bounded, rate.
• Smoother must actually smooth. For strongly anisotropic operators (−ε u_xx − u_yy) or convection-dominated problems, point relaxation smooths in only one direction; error stays oscillatory across the strong coupling, and pointwise multigrid stalls. Cures: line/plane smoothers, semi-coarsening, or matrix-dependent transfers.
• Coarse operator must represent the fine problem. Using naive rediscretization for a jumping-coefficient problem destroys the approximation property; Galerkin coarsening A_H = P ᵀA_h P and operator-dependent interpolation are essential — the origin of algebraic multigrid (AMG, Ruge–Stüben 1987), which builds the hierarchy from A alone with no geometric grid.

Drop any one hypothesis blindly and you lose either O(N) cost or h-independence — usually both.

Why it matters: what multigrid unlocks

Multigrid is the reason large-scale elliptic and parabolic simulation is feasible at all. Its reach:

• Optimal PDE solvers. Pressure-Poisson solves in incompressible CFD, structural mechanics, electromagnetics, and reservoir simulation all lean on multigrid or multigrid-preconditioned Krylov methods for O(N) scalability.
• Preconditioning. A single V-cycle is a superb preconditioner for conjugate gradients / GMRES; MG–CG is robust even when pure MG's constant degrades, and inherits CG's Krylov acceleration.
• Algebraic multigrid. AMG needs only the matrix, so it applies to unstructured meshes and problems with no geometric grid — the workhorse inside solvers like hypre/BoomerAMG and PETSc's GAMG.
• Beyond linear. The full approximation scheme (FAS, Brandt) extends the idea to nonlinear systems, and the same multiscale philosophy seeds multilevel Monte Carlo, wavelet methods, and multigrid-in-time (parareal, MGRIT).

The unifying theory — Xu's subspace correction framework and the Xu–Zikatanov identity — recasts multigrid, domain decomposition, and Schwarz methods as one family, giving sharp convergence identities for all of them.

How multigrid compares to other iterative solvers for the 2D/3D Poisson model problem (N unknowns, spacing h ~ N^(−1/d))

Method	Iterations to fixed accuracy	Work per iteration	Total work	h-dependence
Jacobi / Gauss–Seidel	O(h⁻²) = O(N^(2/d))	O(N)	O(N^(1+2/d))	Deteriorates badly as h→0
Conjugate gradient (CG)	O(h⁻¹) = O(N^(1/d))	O(N)	O(N^(1+1/d))	Deteriorates as h→0
CG + optimal preconditioner	O(1)	O(N)	O(N)	Depends on preconditioner
Multigrid V-cycle	O(1) (ρ ≤ ρ₀ < 1)	O(N)	O(N)	None — mesh-independent
Full multigrid (FMG)	O(1)	O(N)	O(N) to discretization error	None