Iterating, instead of factorising

A rate that does not notice the size

The V-cycle reduces the residual by a factor of ten a cycle at fifteen points and at a hundred and twenty-seven. Jacobi on the same four problems goes from 0.981 to 0.9978, climbing towards one. One of those is a constant and the other is an exponent, and that is the whole distinction the field turns on.

Worth reading first: The error smoothing cannot reach · The rate the condition number predicts.

Every convergence claim in this field so far has the same shape. A method contracts at some factor ρ, the factor has a closed form or a bound, and the closed form or bound contains the condition number — which on the model problem is 4/(π²h²) and grows without limit as the mesh is refined.

Jacobi: ρ = cos(πh), which is 1 − π²h²/2 + O(h⁴). Gauss–Seidel: the square of that. Optimal SOR: a genuinely better exponent, and still 1 − O(h). Conjugate gradients: a bound in √κ, so 1 − O(h) again. Four methods, four rates, and all four go to one as h goes to zero — which is the rate the condition number predicts, holding exactly as it is supposed to.

This essay is the measurement that says the V-cycle does not, and the number it turns on is a constant of about a tenth that does not move across a factor of eight in grid size.

Convergence factor against grid size, with 2 pre-smoothing sweepsThree curves against grid size on a logarithmic horizontal axis. Two rise towards one; the third is a flat line near a tenth.10¹10¹.³10¹.⁶10¹.⁹⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰¹10².²00.250.50.751grid points nresidual reduction per stepJacobiGauss–SeidelV-cycleV-cycle spread, 8× in size0.003Jacobi at n = 1270.99work exponent, fitted0.079the dashed curve is cos(πh), Jacobi's closed formthe flat line is the whole method
Fig. 1 Convergence factor against grid size. Two of the three curves climb towards one as the grid is refined and the dashed prediction beside them is Jacobi’s closed form; the flat line is the V-cycle. Drag the number of pre-smoothing sweeps and watch the flat line fall without tilting.

The measurement

Four grid sizes — 15, 31, 63 and 127 interior points — and on each one the same right-hand side, built from a random vector so that every frequency in the problem is present and nothing converges by luck.

For each, a V-cycle is run to a relative residual of 10⁻¹², and the asymptotic convergence factor is taken from the geometric mean of the last few residual ratios. The factor rather than the iteration count, because the count depends on the starting guess and the factor does not: it is the property of the method, and it is the quantity that has to be bounded away from one and stay there.

The four numbers are 0.1022, 0.1003, 0.1006 and 0.1001. The spread across the whole range is 0.0021, which the build asserts is under 0.08 — a bound loose enough that the assertion is about the shape rather than about the particular arithmetic.

The same four problems given to weighted Jacobi produce 0.9808, 0.9952, 0.9953 and 0.9978, and each is checked against cos(πh) to within 0.02. So the comparison is not between a measured rate and a bound; both curves are measured the same way from the same kind of history, and one of them has a closed form that says where it is going.

Ten cycles, forever

The consequence is worth stating in the units a person planning a computation cares about.

A factor of 0.1 a cycle means ten cycles for ten decades of residual reduction, which is more accuracy than the discretisation error of any of these problems justifies. Ten cycles at fifteen points, ten cycles at a hundred and twenty-seven, ten cycles at a million — the count does not depend on the size, because the factor does not.

Jacobi at 0.9978 needs log(10⁻¹⁰)/log(0.9978) ≈ 10,500 sweeps at 127 points, and about four times that at 255. That is the difference between an algorithm and a table of numbers.

The honest comparison is not in iterations, though, because a V-cycle is much more expensive than a Jacobi sweep. Which is the next section.

A 6-level hierarchy on 63 points, and the coarse operator it impliesRows of dots, each row half the length of the one above it, with lines joining every coarse point to the three fine points it interpolates to.633115731pointsfinest gridone unknown — the recursion bottoms out in a divisionthe coarse operator, two ways‖RAₕP − A₂ₕ‖/‖A₂ₕ‖10⁻¹⁸unknowns / finest grid1.9cycle cost, in fine sweeps14each coarse point reaches three fine ones½, 1, ½ — and the restriction is its transpose
Fig. 2 What a cycle costs. Six levels, and the total unknowns across all of them is under twice the finest grid’s — the geometric series that makes the recursion free, and the reason a V-cycle’s cost is a small multiple of one fine-grid sweep however deep it goes.

Work units, and the exponent

The unit this field uses is the work unit: one smoothing sweep on the finest grid. A V-cycle with ν₁ = 2 and ν₂ = 1 costs three sweeps on the finest level, three on a grid of half the size, three on a quarter, and so on — a geometric series summing to under six work units in one dimension — plus the transfers, which the site counts separately rather than folding in, because a cost model that hides its own overhead is the kind of accounting this site distrusts. The figure’s badge prints the total, and it is about 8.4 for a six-level hierarchy.

So the real claim is about total arithmetic against problem size, and it is fitted rather than quoted. Across the four grid sizes, the total work to reach 10⁻¹⁰ — cycles multiplied by work units — comes to 122, 133, 141 and 145 units. Fitted as a power of n, the exponent is 0.079.

Zero would be perfect: the same number of fine-grid sweeps whatever the size, and therefore a total arithmetic cost proportional to the number of unknowns and to nothing else. The measured 0.079 is not zero because the hierarchy gets one level deeper each time the grid doubles, so a cycle costs marginally more; it is asserted to be under 0.25. The section after next is about what that number is and is not.

Conjugate gradients on the same four problems take exactly 15, 31, 63 and 127 iterations, which is a fitted exponent of exactly 1.00. That number is worth a second look, because it is not the √κ bound doing the work — it is finite termination. The one-dimensional model problem is small enough that CG reaches the exact answer at step n, which is the property an orthogonalisation nobody calls one shows is lost in finite precision at larger condition numbers. Here it holds, and the exponent it produces is a clean 1.

0.079 against 1.00. That is the sentence the field was built to be able to write.

The exponent is an approach to a bound

0.079 is small and the temptation is to leave it there, since it is comfortably inside the assertion. It is worth asking what it is, because the answer changes what the measurement claims.

First, where it comes from. The total work is the cycle count times the work per cycle, and the cycle count is exactly 10 at all four sizes — the factor is near a tenth everywhere, so ten cycles is ten decades everywhere. Every part of the 0.079 is the work per cycle, which is what the paragraph above says and it is worth confirming rather than assuming.

Second, whether the work per cycle grows. The essay’s reason is that the hierarchy gains a level, and that is true; what it leaves out is that the new level is at the bottom, where a grid is half the size of the one above it. workUnits sums (ν₁ + ν₂) times the size of each level over the size of the finest, plus three multiply-adds a point for each restriction and prolongation. Both are geometric series in ½, so both converge. Extending the sweep past the essay’s four sizes:

n      cycles   work per cycle   total    increment
15       10         12.200        122.0      —
31       10         13.355        133.5     1.155
63       10         14.048        140.5     0.693
127      10         14.457        144.6     0.409
255       9         14.694        132.2     0.237
511       9         14.830        133.5     0.136

The increments halve. The limit is 3·2 + 3·3 = 15 work units, which the sixth row is within 1.1% of, and the per-cycle cost cannot pass it however deep the hierarchy goes.

So the fitted exponent falls as the range extends: 0.0788 over the essay’s four sizes and 0.0185 over all six. That is the signature of a quantity approaching a bound rather than growing — a power law fitted to a saturating series returns an exponent that goes to zero as the fit is widened, and the number is a statement about where the range starts rather than about the method.

The right way to say it is therefore not “the exponent is 0.079 and that is nearly zero”. It is that the work per cycle is bounded above by fifteen and the cycle count does not grow, so the total work is bounded, full stop. The exponent is what a log-log fit reports while the bound is being approached, and quoting it as the cost of the extra level makes the result weaker than it is.

The last two rows carry a third thing that was not expected. The convergence factor improves with size on this problem — 0.1014 at fifteen points and 0.0872 at five hundred and eleven — so the cycle count falls from ten to nine, and the fitted exponent of the cycle count is −0.034. That is a small effect and it is in the direction nothing predicted: the standard statement is that the rate is bounded away from one independently of h, and here it is not merely bounded but slightly better on the finer grid, because a deeper hierarchy has a coarsest level that is a smaller fraction of the problem and the coarse solve is exact.

None of that changes the essay’s claim; it strengthens it and removes a number that was doing more work in the reader’s head than the measurement supports.

More smoothing buys a constant

The slider is the number of pre-smoothing sweeps, and what it does is the clearest available statement of what has changed.

At ν₁ = 1 the line sits between 0.167 and 0.181; at ν₁ = 4 between 0.062 and 0.065. More smoothing gives a better factor, as it should — the smoother is removing more of the oscillatory error before the coarse grid is asked for the rest.

Convergence factor against grid size, with 1 pre-smoothing sweepThree curves against grid size on a logarithmic horizontal axis. Two rise towards one; the third is a flat line near a tenth.10¹10¹.³10¹.⁶10¹.⁹⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰¹10².²00.250.50.751grid points nresidual reduction per stepJacobiGauss–SeidelV-cycleV-cycle spread, 8× in size0.014Jacobi at n = 1270.99work exponent, fitted0.08the dashed curve is cos(πh), Jacobi's closed formthe flat line is the whole method
Fig. 3 One pre-smoothing sweep, the cheapest cycle the figure draws. The four factors are 0.172, 0.167, 0.181 and 0.181 across 15, 31, 63 and 127 points.
Convergence factor against grid size, with 3 pre-smoothing sweepsThree curves against grid size on a logarithmic horizontal axis. Two rise towards one; the third is a flat line near a tenth.10¹10¹.³10¹.⁶10¹.⁹⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰¹10².²00.250.50.751grid points nresidual reduction per stepJacobiGauss–SeidelV-cycleV-cycle spread, 8× in size0.0046Jacobi at n = 1270.99work exponent, fitted0.078the dashed curve is cos(πh), Jacobi's closed formthe flat line is the whole method
Fig. 4 Three: 0.083, 0.080, 0.079 and 0.080, a spread of 0.004 — and 0.0010 across the last three.

What it does not do is tilt the line, and how flat the line is turns out to depend on the slider in a way the summary above does not carry. The spread across the four grid sizes reads 0.014, 0.0021, 0.0040 and 0.0028 at ν₁ = 1, 2, 3 and 4 — so at one sweep it is seven times what it is at two, and it is the only setting where the largest factor is at the finest grid: 0.167 at 31 points against 0.181 at 63 and 127. One sweep of weighted Jacobi does not quite remove the oscillatory error the coarse grid cannot carry, and what is left of it is a little worse on a finer grid, which is the field’s whole failure mode appearing in miniature at the setting nobody uses.

At two sweeps and above it is gone: 0.0021, 0.0040 and 0.0028 across a factor of eight in size, with no ordering by grid. So more smoothing improves the constant and — up to the first extra sweep — buys the flatness as well. The field’s difficulty was a rate that worsened with size, and adding sweeps to Jacobi does not touch that any more than adding a fifth sweep here tilts a line that is already flat.

And the sweeps buy less each time. The factor at the finest grid reads 0.181, 0.099, 0.080 and 0.063 as ν₁ goes 1, 2, 3, 4 — ratios of 1.83, 1.24 and 1.27 for one extra work unit each. The first extra sweep is worth nearly a factor of two and the third is worth a quarter, on a cost that rises by the same amount every time, which is the whole of why the answer is one or two.

The economics point the same way. Each extra sweep costs a work unit and buys a better factor, so there is an optimum in total work rather than in factor — and it is usually one or two, because the factor improves geometrically and the cost improves linearly. The figure lets the reader find it; the site does not assert where it is, because that is a property of the smoother, the problem and the machine rather than of the mathematics.

What “optimal” means, and what it does not

The word attached to methods of this kind is optimal complexity, and it means one thing specifically: the work to solve the problem to a fixed accuracy is proportional to the number of unknowns. Not to n log n, not to n^1.5 — to n, which is the cost of writing the answer down.

Three qualifications, and each of them is a place where a reader could be misled.

It is not the fastest method on a small problem. A V-cycle has real overhead, and on a few thousand unknowns a well-tuned direct solver beats it comfortably. Optimal complexity is an asymptotic statement, and the order decides the memory records what happens to asymptotic statements at the sizes anybody draws: nested dissection has the better exponent and loses to minimum degree at every size on that essay’s slider. A complexity class is a statement about a limit.

It is not robust by default. The rate measured here is on the Laplacian, which is the problem the components were designed for. Change the operator — add strong convection, or anisotropy, or jumping coefficients — and the smoother stops smoothing in the direction the coupling is strong, or the coarse grid stops representing the modes the smoother leaves, and the flat line tilts. Everything in the field beyond this point is about repairing one or other of those: smoothing a whole line at once, semi-coarsening — which is the repair that is exact along one axis and blind to the other — matrix-dependent transfers, and algebraic multigrid, which gives up on grids entirely and builds its levels from the matrix.

And it is a claim about the solve, not about the setup. On a structured grid the hierarchy is free. On an unstructured one it has to be built — the same asymmetry a tree preconditioner has, where the object is cheap to apply and the construction is the whole cost — and building it is not proportional to anything as tidy as n in practice.

Starting at the bottom, which removes the last logarithm

There is one more idea in the field and it is cheap to describe, because every piece of it has already been built.

A V-cycle starts on the finest grid with whatever initial guess is available, usually zero. That guess is bad, so the first few cycles are spent getting from “nothing” to “roughly right”, and only after that does the tenth-per-cycle rate start doing useful work. The count above — ten cycles for ten decades — includes that opening.

Full multigrid removes it. Solve the same problem on a coarser grid first, where at the bottom it is a division. Interpolate that answer up one level as the initial guess there, and run one V-cycle. Interpolate the result up again, run one V-cycle, and so on to the finest grid. Each level starts from an answer that is already accurate to the discretisation error of the level below it, so one cycle a level is enough to bring it to the discretisation error of its own.

The accounting is the interesting part. One V-cycle per level, on grids that halve, is another geometric series — so the total is a small multiple of one V-cycle on the finest grid, and the result is an answer accurate to the discretisation error of the problem in O(n) work with no logarithm at all. That is as good as this kind of statement gets: the cost of solving the equations is a constant multiple of the cost of writing down the answer.

Two things are worth being careful about, since the claim sounds like more than it is.

“Accurate to the discretisation error” is not “accurate”. The answer full multigrid returns is as close to the true solution of the differential problem as the discrete problem’s own answer is, and no closer. Solving the discrete system to fifteen digits when the discretisation is accurate to four is work spent on digits that mean nothing — which is a general point about iterative methods and is the reason a fixed tolerance is often the wrong stopping rule. It is also a distinction this site usually declines to make, because every other essay here is about the arithmetic rather than about the model it is discretising.

And the interpolation between levels is not the cycle’s interpolation. The transfer used to move an answer up a level should be more accurate than the one used to move a correction, because a correction is only ever required to be smooth and an answer is required to be right. Codes typically use cubic interpolation for the first and the linear ½, 1, ½ for the second, which looks like an inconsistency and is a considered one.

Neither of those is measured here. Full multigrid is named rather than built, because the claim this phase set out to establish is the one the flat line makes — that there is a method whose rate does not depend on the size — and the remaining factor is a constant on top of a constant.

Convergence factor against grid size, with 4 pre-smoothing sweepsThree curves against grid size on a logarithmic horizontal axis. Two rise towards one; the third is a flat line near a tenth.10¹10¹.³10¹.⁶10¹.⁹⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰¹10².²00.250.50.751grid points nresidual reduction per stepJacobiGauss–SeidelV-cycleV-cycle spread, 8× in size0.0038Jacobi at n = 1270.99work exponent, fitted0.077the dashed curve is cos(πh), Jacobi's closed formthe flat line is the whole method
Fig. 5 Four pre-smoothing sweeps rather than two. The flat line sits lower and is still flat, which is the distinction this essay is about: every knob in the field moves the constant, and only the coarse grid touches what happens as the problem grows.

What is asserted, and what refuses

The V-cycle converges on all four grids, and its factor is under a quarter on each, at every position of the sweep slider.

The spread across a factor of eight in size is under 0.08. This is the claim the whole field exists for, and it is stated as a bound on a spread rather than as an equality, because the factor does move slightly and describing it as constant would be describing something false.

The smoother’s factor climbs, and each of its four values is checked against cos(πh) — so the curve the V-cycle is being compared against is itself checked against a closed form rather than merely observed.

The work exponent is under 0.25, fitted from the measured work at four sizes.

And the two refusals are the two halves of the method. A cycle with the coarse-grid correction removed is ν₁ + ν₂ Jacobi sweeps, and its factor is 0.968 — the assertion that it converges like a V-cycle is fed that and refuses it. A coarse-grid correction with no smoothing at all, applied to the most oscillatory vector the grid holds, leaves 134% of the error it found — and the assertion that it removes half of it is refused too. Neither half is a method. Both are asserted to fail, in the direction each one fails in.

What the phase leaves open

Two things, and both are honest gaps rather than omissions.

The proof. Two-grid convergence for this problem can be proved by Fourier analysis, and the V-cycle’s by a rather more involved argument. Neither is run here. What is run is the measurement, checked against closed forms wherever a closed form exists — which is the same standing this site gives the conjugate-gradient bound, where the bound is proved elsewhere and checked here at every step.

The general case. Everything measured in this essay’s own tables is one-dimensional, and the figures at the end reach into two — the anisotropic V-cycle among them — which is where the field’s real difficulty lives. The two-dimensional model problem is the one where multigrid’s advantage is overwhelming rather than merely large, since the direct alternative’s fill grows as n^1.5 there and the iteration counts are worse — and it is also where the figures become expensive to generate at build time, since every drag position runs the whole solve. The one-dimensional case carries every claim in this essay and none of the claims about anisotropy, which is where the field’s real difficulty lives and where a later phase would have to go.

What is settled is the thing the field was missing. Every other method here improves a constant in front of a growth. This one removes the growth, and the measurement that says so is a flat line across a factor of eight and a fitted exponent of 0.079.

What links here

Computed from the collection, not written here: the essays that point at this one.

Named objects

A flat tag is an object no other essay names yet.

Convergence factorJacobi iterationOptimal-complexityThe V-cycleWork units