The error smoothing cannot reach

A rate that is known in advance establishes that the stationary methods converge at rates the model problem gives in closed form, and that those rates are bad in a specific and worsening way. Jacobi contracts at cos(πh). On a grid of 24 interior points that is 0.99211, so the residual falls by less than a per cent a sweep, and 3,000 sweeps do not reach 10⁻¹³. Refine the grid and it gets worse, because the rate goes to one like 1 − O(h²).

Optimal SOR improves the exponent of that statement — 128 iterations where Jacobi takes thousands — and does not remove it. Conjugate gradients improve it again, to √κ rather than κ, and do not remove it either. Every method in the field so far needs more work as the grid is refined, and the reason is always the same: the condition number of the discrete Laplacian grows like h⁻², and every one of these methods is governed by it.

This essay is the first step of the argument that removes it, and the step is a reframing rather than an algorithm. The claim is that the smoother is not bad at everything, and that what it is good at and what it is bad at can be separated exactly.

One sweep, one mode, one number #

The measurement behind the figure is worth stating precisely, because it is the cleanest instance of this site’s two routes habit anywhere on it.

Weighted Jacobi is x ← x + ω D⁻¹(b − Ax). Subtract the exact solution from both sides and the error obeys e ← (I − ω D⁻¹A) e, so a sweep is a matrix applied to the error and nothing else. For the one-dimensional discrete Laplacian, D is a multiple of the identity, and the eigenvectors of A are the grid sines sin(kπjh) for k = 1 … n. So every mode is an eigenvector of the sweep, and applying a sweep to a pure mode multiplies it by a scalar.

That scalar is available in closed form. The eigenvalue of A at mode k is (4/h²)sin²(θ/2) with θ = kπh, and the eigenvalue of the iteration matrix is therefore

$$\mu(\theta) = 1 - 2\omega\sin^2(\theta/2)$$

which is the curve in the figure. The dots are what happens when a sweep is actually applied to that mode and the ratio taken. They agree to 5.6·10⁻¹⁶, which is one operation’s worth of rounding, and the agreement is asserted at every position of the slider rather than at the one the essay was written against.

Nothing here is an approximation of anything. The mode is an eigenvector exactly, the eigenvalue is a sine exactly, and the measurement is a ratio of two norms. It is the rare case where a numerical claim is a claim about arithmetic rather than about a method.

The two halves, and which one is the problem #

The curve says two things and they point in opposite directions.

At the left, where θ is small, sin²(θ/2) is about θ²/4 and μ is about 1 − ωθ²/2. The smoothest mode on a grid of 63 points has θ = π/64, so μ = 0.9992: a sweep removes eight parts in ten thousand of it. A hundred sweeps remove 7%. This is the whole of the bad news about stationary methods, written as a property of one mode rather than as a rate, and it gets worse as the grid is refined because θ for the smoothest mode is πh.

At the right, where θ approaches π, sin²(θ/2) approaches 1 and μ approaches 1 − 2ω. At ω = 2/3 that is −1/3. The mode changes sign every sweep and loses two thirds of its size while doing it. Three sweeps take it to 3.7% and six to 0.14%.

So the error after a few sweeps is not small. It is smooth, which is a different and much more useful statement, and it is the one the next essay builds on.

Why the weight, and why two thirds #

Plain Jacobi is ω = 1, and the figure at that weight is the argument for not using it.

At ω = 1 the highest mode has μ = −1. It is not damped at all — it oscillates, sweep after sweep, at full size. Plain Jacobi does reduce the middle of the spectrum well, which is why it converges at all, and it is useless as a smoother because the thing a smoother has to guarantee is that nothing oscillatory survives.

The quantity that guarantee is stated in is the smoothing factor: the worst damping over the oscillatory half θ ∈ [π/2, π]. Because μ is monotone in θ, the worst case is at one end or the other, and the maximum is

$$\mu^* (\omega) = \max\bigl(|1 - \omega|,\ |1 - 2\omega|\bigr)$$

with no search required. The first branch falls as ω rises and the second rises; they cross where 1 − ω = 2ω − 1, which is ω = 2/3, and the value there is exactly one third.

That “exactly” is not an approximation of 0.3333. The two branches are equal at the crossing and both evaluate to 1/3, so the optimum is a rational number and the site checks it as one — against a search over three thousand weights, which is a grid chosen so that 2/3 is a point on it. The first version of that search used two thousand steps, which misses 2/3, and reported the optimum as 0.3335. A discrepancy of two parts in ten thousand reads as a numerical artefact and was entirely a property of the search.

What the coarse grid is for, stated before it exists #

The next essay builds the hierarchy. The reason it can be built at all is visible here, and it is one sentence: the modes a smoother cannot damp are the modes a coarser grid can represent.

A mode with θ near zero varies slowly across the grid. Throw away every second point and it is still there — sampled at half the resolution, and perfectly recognisable. A mode with θ near π alternates from point to point, and throwing away every second point does not sample it, it destroys it.

So the two halves of the spectrum split along exactly the same line twice, for two unrelated reasons. The smoother handles the half that a coarse grid cannot see. The coarse grid can see the half the smoother cannot handle. Nothing arranged that; it is a consequence of both facts being about the same quantity, the number of mesh intervals per oscillation.

The number this replaces #

It is worth being precise about what is being claimed and what is not, because the field so far has been about rates and this essay has not produced one.

A smoothing factor of 1/3 is not a convergence rate. Weighted Jacobi at ω = 2/3 converges more slowly than plain Jacobi does, and considerably more slowly than Gauss–Seidel: its contraction factor on the 63-point grid is 0.998, against Jacobi’s 0.9995 and optimal SOR’s 0.905. Judged as a solver it is the worst method in the field.

Judged as a smoother it is the best of the three, and the two judgements are about different halves of the spectrum. This is the reframing the whole method rests on: stop asking how fast the residual falls, and ask which components of the error are left.

The other smoothers, and why this one is drawn #

Weighted Jacobi is not what a production code smooths with, and it is what this essay draws, so the gap is worth closing.

Gauss–Seidel is a better smoother and a worse figure. Its iteration matrix is not symmetric, its eigenvectors are not the grid sines, and its damping is therefore not a number one can attach to a mode — which removes every closed form above and with them the two-route check that makes the figure evidence rather than an illustration. Measured rather than analysed, its smoothing factor on this problem is about 0.5 at the worst oscillatory mode with no weight at all, so it needs no ω and does better than weighted Jacobi at the same cost.

Red-black Gauss–Seidel recovers the analysis. Colour the grid points alternately, update all the red ones and then all the black ones, and the resulting operator does have the grid modes as eigenvectors — in pairs, since the colouring couples the mode at θ to the mode at π − θ. The smoothing factor comes out at 1/4, better than weighted Jacobi’s 1/3, and the pairing is the reason its analysis is a page rather than a line. It is also the version that parallelises, since every red point can be updated at once, which is why it is the one in the codes.

The choice does not change anything in this essay. All three damp the oscillatory half and none of them touches the smooth end, and the entire argument of the field turns on that shape rather than on which constant sits at the right-hand edge of the plot. Weighted Jacobi is drawn because it is the one whose every claim is a sine — and because a smoother chosen for the clarity of its analysis, with the better options named, is a more honest figure than the best smoother drawn with its analysis omitted.

One consequence is worth carrying forward. Because the smoothing factor is what matters and the convergence rate is not, the usual instinct about relaxation methods inverts: the parameter that makes a stationary method converge fastest is not the parameter that makes it smooth best. Optimal SOR is tuned to minimise the spectral radius of the whole iteration matrix, which is dominated by the smooth modes it will never fix, and as a smoother it is unremarkable. Weighted Jacobi at ω = 2/3 is tuned to minimise the maximum over the oscillatory half alone, and is the worse solver of the two by a wide margin. Choosing between them requires knowing which of the two jobs the method is being asked to do, and the field before this essay had only one job to offer.

What is asserted, and what refuses #

Every mode’s damping, against its closed form. Twelve modes on a 63-point grid, at every position of the weight slider, agreeing to 5.6·10⁻¹⁶. The two routes share nothing: one is a norm ratio after a sweep, the other is a sine.

The optimum is exactly one third, checked against a search over the weight, and the closed form for the smoothing factor is checked against a scan over the oscillatory frequencies it maximises over. Neither is trusted from the algebra.

The smoothest mode survives whatever the weight is. This is asserted at every slider position, and it is the assertion that would catch a weight that appeared to solve the problem. No ω damps the left-hand end, and a figure implying one did would be drawing something false.

The refusal: the claim that weighting Jacobi makes no difference to the smoothing factor is fed to the assertion and must be rejected. µ*(1) = 1 and µ*(2/3) = 1/3, so a build in which those are reported as equal fails.

And one more, which belongs to the next essay and is worth stating here because it is what stops this one from being a story about a clever weight: a cycle with the coarse-grid correction removed is ν₁ + ν₂ Jacobi sweeps in a costume, and its convergence factor on the 63-point grid is 0.968. The assertion that it converges at the multigrid rate is fed that number and refuses it.

Where this leaves the field #

Three of the site’s iterative essays end with a version of the same sentence: the method is good, the rate is provable, and the rate still degrades as the problem grows.

This one ends differently. The smoother’s rate degrades and its smoothing factor does not — 1/3 at every grid size, because the closed form for it has no h in it at all. That is a quantity which does not care how fine the mesh is, and it is the first such quantity in the field.

What remains is to do something with the half of the error that is left, and the whole of the next essay is that: restrict the residual to a coarser grid, solve the same problem there, and interpolate the correction back. On the coarse grid the modes that were smooth are no longer smooth — halving the points doubles every frequency relative to the mesh — so the smoother that could not see them there can. Recurse, and every frequency in the problem is somebody’s high frequency.

A note on the model problem, and on what generalises #

Every closed form in this essay belongs to the one-dimensional discrete Laplacian, and a reader is entitled to ask how much of the argument survives leaving it.

The exact statements do not. On a general matrix the eigenvectors are not sines, they are not known, and the damping factor of a mode is not a number anybody can write down. The figure above could not be drawn.

The structure does, and it is what the field is built on. Local relaxation reduces the components of the error that vary rapidly between neighbouring unknowns, because it only ever looks at neighbouring unknowns; it cannot reduce components that are nearly constant across a neighbourhood, because to such a component the local equation is nearly satisfied already. That argument mentions no eigenvector and no sine. It is why multigrid works on unstructured problems, and why algebraic multigrid — which builds its coarse levels from the matrix rather than from a grid — is possible at all.

The model problem is here for the reason it is here everywhere else on this site: it is the case where the claim can be checked against something that was not measured. A rate observed on a general matrix is an observation. A rate checked against a sine is evidence that the machinery is doing what it is described as doing, and the machinery is then run on cases where no such check exists.

That is the same trade the rate the condition number predicts makes with the conjugate-gradient bound, and the same one a rate that is known in advance makes with the three stationary rates. The pattern is worth naming, since this phase leans on it three more times: build the machinery where the answer is known, assert against the answer, and only then take it somewhere the answer is not known.