Relaxation is not stupid.

It is nearsighted.

That is the mental picture I wish I had learned first. A Jacobi or Gauss-Seidel sweep is very good at killing error that alternates from one grid cell to the next. It is bad at killing error that looks like a slow hill across the whole domain. Not because the slow hill is small. Because a local stencil cannot feel much difference between neighboring points on that hill.

Multigrid does not fix this by inventing a heroic smoother.

It changes the grid.

On a coarser grid, the same smooth error has fewer points per wavelength. It is less invisible. The local solver gets another chance to see it.

That is the whole trick, and it is a good one.

The Local Solver Has Bad Eyesight

Use the one-dimensional Poisson problem with zero boundary values:

\[-u''(x)=f(x),\qquad u(0)=u(1)=0.\]

On a uniform grid, the standard second-difference system is

\[\frac{2u_i-u_{i-1}-u_{i+1}}{h^2}=f_i.\]

Weighted Jacobi updates a current guess by looking at the residual through the diagonal:

\[u_i^{new} =u_i+\omega\frac{h^2}{2} \left(f_i-\frac{2u_i-u_{i-1}-u_{i+1}}{h^2}\right).\]

For the error, this is a filter. A sine mode with angle theta is multiplied by

\[1-2\omega\sin^2(\theta/2).\]

High-frequency error has theta near pi, so a weight near 2/3 damps it hard. Low-frequency error has theta near zero, so the multiplier is close to one. The sweep barely notices it.

This is why a residual plot from plain relaxation can feel insulting. The first few sweeps do something real. Then the curve flattens, not because the remaining error is mysterious, but because it is smooth.

Brandt’s 1977 multilevel paper made this scale separation operational for boundary-value problems.1 The tutorial treatment by Briggs, Henson, and McCormick popularized the clean engineering recipe: smooth, compute residual, restrict to a coarse grid, solve a correction problem, prolong the correction, and smooth again.2

The method sounds recursive because the obstruction is recursive. A coarse problem has its own smooth error. So it gets an even coarser problem.

The Residual Is a Complaint About the Error

Suppose the current guess is u and the exact discrete solution is u*. The error is

e = u* - u

The residual is

r = f - A u

Because A u* = f, those two objects satisfy the correction equation

\[A e = r.\]

That equation is the permission slip for multigrid. The solver does not know the error, but it can compute the residual. If the residual is restricted to a coarser grid, a cheaper problem can estimate the missing smooth correction.

One V-cycle looks like this:

smooth on the fine grid
compute residual
restrict residual to the coarse grid
approximately solve A_coarse e_coarse = r_coarse
prolong e_coarse back to the fine grid
correct the fine-grid solution
smooth again

There are many choices hiding in those verbs. Jacobi, Gauss-Seidel, Chebyshev, block smoothers, Schwarz smoothers, aggregation, Galerkin coarse operators, geometric coarsening, algebraic coarsening. Modern algebraic multigrid takes the same idea and asks how to build the hierarchy directly from the matrix, without an obvious physical grid.3

But the base bargain remains:

smoother removes oscillatory error
coarse grid removes smooth error

Local Fourier analysis is the microscope for this bargain. Thompson, Brown, and He describe LFA as a way to investigate and tune multigrid methods, including Jacobi and Chebyshev smoothing choices for two-grid schemes.4 The details get sophisticated quickly, but the first diagnostic is humble: which wavelengths survive one sweep?

A Browser-Sized V-Cycle

The lab below solves a 1D Poisson system where the exact discrete solution is known. The right-hand side is built by applying the discrete operator to that known solution, so a tridiagonal direct solve can audit the whole experiment.

The implementation uses:

  • weighted Jacobi smoothing;
  • full-weighting residual restriction;
  • linear interpolation for prolongation;
  • direct tridiagonal solves on the coarsest grid;
  • a recursive V-cycle.

The default has 127 interior grid points, Jacobi weight 0.67, two pre/post smoothing sweeps, six V-cycles, and a mixed-mode target solution. The audited numbers from the same JavaScript file are:

V-cycle relative residual        4.23e-7
Jacobi relative residual          0.159
fine-grid work units              46.6
high-mode energy after smoothing  1.33%
low-mode energy after one V       0.211%
direct-solve audit error          5.83e-14
36-case audit grid                36 / 36 passed

The comparison is not meant to be a production benchmark. It is a scale experiment. Jacobi gets roughly the same fine-grid work budget as the V-cycles. It smooths. It does not change scales.

Deterministic 1D Poisson solver for the mechanism, not a replacement for a production multigrid package. The executable artifact is assets/js/multigrid-error-lab.js.

The Coarse Grid Is Not a Smaller World

A coarse grid is not just a cheap approximation to the fine grid.

It is a change of coordinates for the error.

That distinction matters because it explains why a tiny coarse solve can be worth so much. The coarse correction does not need to reproduce every fine-scale wiggle. The smoother already handled those. The coarse solve is there for the slow components that were invisible to the local iteration.

This also explains why multigrid can fail. If the smoother does not damp the right high-frequency components, the coarse grid receives garbage. If the restriction and prolongation operators do not represent the slow error well, the correction misses the thing it was hired to fix. If coefficients are anisotropic or discontinuous, a naive geometric hierarchy may stop matching the physics. Then the method becomes an art again.

The recent high-order discontinuous Galerkin study by Saye is a nice reminder that this is still active engineering: using a standard V-cycle template, it tests different block smoothers across Poisson, interface, and Stokes problems, and reports cases where 4 to 8 iterations reduce residuals by 10 orders of magnitude.5 The V-cycle is old. The choice of smoother and hierarchy is still alive.

Why This Feels Like Cheating

Direct solvers eliminate unknowns.

Krylov methods build better and better polynomials on the spectrum.

Multigrid does something slightly more physical. It asks what wavelength the remaining error has, then changes the measuring device.

On the fine grid, a slow error mode is a broad hill. Jacobi sees neighboring points that are almost equal and shrugs. On a coarse grid, that same hill is sampled with fewer points. It bends more sharply. A local correction can touch it.

Then the correction is interpolated back, and the fine grid smooths away the interpolation scars.

That last detail is part of the beauty. Multigrid is not a single clever move. It is a division of labor:

fine grid: remove oscillations
coarse grid: remove hills
finer grid again: clean up the stitching

If the problem really has elliptic, multiscale structure, that division can be almost embarrassingly efficient. The cost of visiting all coarser grids is a geometric series. In one dimension, the total number of unknowns across the hierarchy is less than twice the fine-grid unknown count. In higher dimensions, the same geometric idea survives even more comfortably.

This is why multigrid has the aura of a trick, but it is not magic.

It is scale accounting.

The Useful Lesson

The lesson I want to keep is not “use multigrid everywhere.”

It is sharper:

slow convergence is often a scale mismatch

If the remaining error is smooth, do not keep asking a local fine-grid iteration to discover it from tiny neighbor differences. Move the complaint to a grid where it has the right wavelength. Solve there. Come back. Smooth the seams.

The error changes grids because the solver has eyes at different scales.

  1. Achi Brandt, “Multi-Level Adaptive Solutions to Boundary-Value Problems”, Mathematics of Computation 31(138), 1977. DOI. ↩

  2. William L. Briggs, Van Emden Henson, and Steve F. McCormick, A Multigrid Tutorial, 2nd edition, SIAM, 2000. DOI. ↩

  3. Jinchao Xu and Ludmil T. Zikatanov, “Algebraic Multigrid Methods”, Acta Numerica 26, 2017. arXiv. ↩

  4. Jeremy L. Thompson, Jed Brown, and Yunhui He, “Local Fourier Analysis of P-Multigrid for High-Order Finite Element Operators”, 2021. arXiv. ↩

  5. Robert I. Saye, “A comparative study of efficient multigrid solvers for high-order local discontinuous Galerkin methods: Poisson, elliptic interface, and multiphase Stokes problems”, 2024. arXiv. ↩