The same problem on a coarser grid

The error smoothing cannot reach ends with an error that is smooth and not small. Three weighted Jacobi sweeps take the oscillatory part of it from 82% of the energy to 2%, and leave the smooth part at 99% of the size it started at. Running more sweeps is pointless — the smooth part is being multiplied by 0.9992 each time — and the question is what to do with what is left.

The answer is to solve for it somewhere else. A smooth function on 63 points is a smooth function on 31 points that has been sampled twice as finely; nothing about it needs the fine grid. So restrict the residual to the coarse grid, solve the same equation there, interpolate the correction back, and add it on.

Every clause of that sentence is a matrix, and this essay is about which matrices and why those.

Interpolation, which is three numbers #

Take a fine grid of n = 2ᴸ − 1 interior points and a coarse grid of (n − 1)/2. The sizes are chosen so that every coarse point sits exactly on a fine point — the coarse grid is the fine grid with every second interior point removed — and that choice is what makes the transfer operators trivial rather than a projection problem.

Interpolation from coarse to fine, written P, does the obvious thing. A fine point that coincides with a coarse point takes its value. A fine point between two coarse points takes their average. As a matrix, every column of P has the entries ½, 1, ½ in three consecutive rows, and that is the whole operator.

Restriction from fine to coarse is not the obvious thing, and the difference matters. The obvious choice is injection — take the value at the coincident fine point and discard the two neighbours. The choice made here is full weighting: a quarter, a half, a quarter, which is P transposed and halved.

The halving is a normalisation, so that a constant restricts to a constant. The transpose is the part worth arguing for, and there are two arguments.

The first is structural. If R = ½Pᵀ then the coarse operator built below comes out symmetric whenever the fine one is, which means the coarse problem can be handed to the same machinery as the fine one — the same smoother, the same recursion, eventually the same conjugate-gradient wrapper. Injection does not give that.

The second is about what restriction is being asked to do. The residual on the fine grid still has an oscillatory component; smoothing reduced it to 2% rather than to nothing. Injection samples that component and reports it as though it were smooth — the classic aliasing failure, where a high-frequency signal read at half the sampling rate looks like a low-frequency one. Full weighting averages first, which attenuates the oscillatory part before it can be misread. It is a low-pass filter, applied for exactly the reason low-pass filters are applied before decimation anywhere else.

The coarse operator, and the identity that makes it legitimate #

Something has to play the part of A on the coarse grid, and there are two candidates.

Rediscretise. Write down the same differential problem on the coarse mesh and discretise it there. This is what a person would do and it needs a grid to do it.

Galerkin. Take the product R A P. This needs no grid at all — it is three matrices multiplied — and it is what algebraic multigrid must use, since there is no mesh to rediscretise on.

On this problem they are the same matrix. Not close: equal, entry for entry, to 15 digits, at every level of a six-level hierarchy. The site computes both and prints the relative difference on the figure, and the difference is exactly zero at every level.

That is the identity the recursion rests on. If R A P were merely an approximation to the coarse discretisation, then the coarse problem would be a different problem and the correction computed there would be a correction to something else. Because they are equal, the coarse grid is solving the same equation — and the recursion can be applied again, and again, with the same justification each time.

The scaling is where this gets fussy, and the fussiness is worth one paragraph because it is a real trap. iterative.js uses the unscaled Laplacian, tridiagonal(−1, 2, −1), because every claim it makes is invariant to a scalar multiple: a stationary method’s iteration matrix does not change if the matrix is multiplied by two. The Galerkin identity is not invariant. With the unscaled operator, R A P comes out as one quarter of the coarse tridiagonal, and the quarter is the h² that was divided out. So the operators in this field carry their h⁻², and laplace1d(n) is (n+1)² times poisson1d(n). Every rate in these three essays is unchanged by that and the coarse operator is only the coarse operator with it.

The cycle #

Put the pieces together and a two-grid correction is six lines.

Smooth ν₁ times on the fine grid. Compute the residual r = b − Ax. Restrict it: r_c = R r. Solve A_c e_c = r_c on the coarse grid. Interpolate: e = P e_c. Correct: x ← x + e. Smooth ν₂ times.

The step that is cheating is “solve on the coarse grid”, and the whole of multigrid is the observation that it does not have to be a solve. It can be another two-grid correction, and the recursion bottoms out when the grid is one point and the solve is a division. That is the V-cycle, and it is the subject of a rate that does not notice the size.

What is worth noticing here is what the recursion costs. Each level has half the points of the one above, so the total number of unknowns in the whole hierarchy is under twice the finest grid’s, however many levels there are — the readout on the figure above prints the ratio and it converges to 2 from below. One V-cycle therefore costs a bounded multiple of one fine-grid sweep no matter how deep it goes, which is a geometric series doing all the work and is the reason the recursion is free.

Why the smooth error is the one that transfers #

There is a symmetry in the argument that is easy to state and easy to get backwards, so it is worth saying twice.

Interpolation is exact on the modes the coarse grid can represent and wrong on the others. A smooth mode interpolated from the coarse grid comes back as very nearly the same smooth mode. An oscillatory mode does not exist on the coarse grid to be interpolated from, and anything the coarse grid sends up in its place is smooth.

So the correction P e_c is always smooth, whatever the coarse solve returns. That is a limitation and it is the reason the smoother is not optional: whatever oscillatory error is present at the point the correction is applied, the correction cannot remove it, and there is no arrangement of the coarse grid that changes this.

The measured version of that statement is one of this field’s two refusals, and it came out sharper than expected. Take the most oscillatory vector the grid holds and apply a coarse-grid correction to it with no smoothing at all. The expectation is that nothing happens. What happens is that 134% of the error is left — the correction makes it worse. Full weighting nearly annihilates that mode, so the residual restricted to the coarse grid is tiny; the coarse solve then divides that tiny residual by a small coarse eigenvalue and interpolates the result back as a smooth vector that was not there before. The build asserts the increase rather than the expectation, because the expectation was wrong.

The variational property, which is what the transposes buy #

The choice R = ½Pᵀ has a consequence that is worth naming because it is the reason the same three matrices keep reappearing in codes that share nothing else.

When the fine operator is symmetric positive definite and the restriction is the transpose of the interpolation, the coarse-grid correction is an orthogonal projection in the A-norm. That is, among all corrections that could be built by interpolating something from the coarse grid, the one the coarse solve produces is the one that minimises the A-norm of the remaining error. It is not merely a good correction; it is the best available correction of that shape.

Three things follow, and the third is the practical one.

The correction cannot make the A-norm of the error larger. A projection does not increase the norm it projects in. This is why the 134% increase measured above is not a contradiction: that measurement is in the ordinary Euclidean norm, on a vector chosen to be the worst case, and the A-norm of the same error does not increase. Two norms disagreeing about whether something got worse is exactly the trap the rate the condition number predicts records about the conjugate-gradient bound, where a bound in the A-norm quoted against the 2-norm is being compared against a quantity it was never about.

The coarse operator inherits definiteness. RAP with R = ½Pᵀ is ½PᵀAP, and PᵀAP is positive definite whenever A is and P has full column rank. So the coarse problem is a problem of the same kind, which is what lets it be handed to the same smoother — and would not be true of an arbitrary pair of transfer operators.

And the whole cycle can be used as a preconditioner. A symmetric V-cycle — the same number of smoothing sweeps before and after, with a symmetric smoother — is a symmetric positive definite operator, so conjugate gradients can be wrapped around it. That combination is what a great many production codes actually run, and it is a strictly better object than either half: the V-cycle supplies a rate that does not depend on the grid, and the Krylov wrapper supplies robustness when the V-cycle’s assumptions are only approximately true. It is the same move changing the condition number on purpose makes with an incomplete factorisation, applied to a preconditioner whose quality does not degrade with size.

The variational property is also where the field’s aesthetic argument comes from, and it should be resisted rather than repeated. Multigrid is often introduced as elegant, and the elegance is real — three matrices, nine entries, a projection that is optimal in the natural norm. Elegance is not evidence. What makes the method worth the field’s attention is a measured convergence factor that does not move when the grid is refined by a factor of eight, and that number is in the next essay rather than this one.

What is asserted here #

The Galerkin operator is the coarse discretisation, at every level of a six-level hierarchy, to a relative difference of exactly zero. Two routes, one of which multiplies three matrices and the other of which writes down a stencil.

Restriction is half the transpose of interpolation, checked as a matrix identity, and the coarse operator is therefore symmetric — also checked, since the symmetry is what lets the recursion reuse the machinery.

The hierarchy costs under twice the finest grid. Asserted at every depth on the slider, because it is the claim that makes the recursion free and it would be quietly false if the coarsening ratio were wrong.

And the coarse grid cannot see the finest mode. Fed the most oscillatory vector on the grid, the correction leaves more error than it found; the assertion that it removes half of it is refused.

The pieces, and what each one costs #

A final accounting, since the field’s other essays are careful about cost and this one has introduced three new operators.

Interpolation and restriction are three multiply-adds per point, on grids that halve, so their total across a whole cycle is a small multiple of one fine-grid sweep. The coarse operators are formed once and reused. The only per-cycle work that is not a sweep is the transfers, and the figure’s readout counts them separately from the smoothing rather than folding them into a single number — a cost model that hides its own overhead is the kind of accounting this site exists to distrust.

The one genuinely awkward cost is the one nobody draws: setup. On a structured grid the hierarchy is free, because the sizes and the stencils are known. On an unstructured problem, or in algebraic multigrid, the coarse levels have to be chosen — which unknowns survive coarsening, and with what weights — and that choice costs real time before a single cycle runs. It is the reason algebraic multigrid is sometimes slower than a simpler method on a problem solved once, and much faster on a problem solved a thousand times with the same structure. The same distinction the order decides the memory draws about symbolic factorisation applies here for the same reason: a cost paid once amortises over a sequence of solves and does not amortise over one.

A remark on what has and has not been shown #

Nothing in this essay is a convergence proof. Two-grid convergence can be proved for this problem, by splitting the error into the modes the coarse grid represents and the modes it does not and bounding each separately, and the proof is a page of Fourier analysis that this site does not run.

What has been shown is that each piece is what it claims to be: the transfers are adjoint, the coarse operator is the coarse discretisation exactly, the correction is smooth by construction, and the smoother handles what the correction cannot. The measurement that the whole thing converges — and at what rate, and whether that rate cares about the grid — is the next essay’s, and it is a measurement rather than a proof there too.

That division is deliberate and it is this site’s usual one. The pieces are checked against closed forms; the assembly is measured. A proof of two-grid convergence quoted without either would be a third thing, and the least useful of the three.