The coarse problem is a different problem

A rate that does not notice the size ends with a list of what had not been measured, and the first item on it is this one. Every number in the multigrid field — the convergence factor of 0.10, the smoothing factor of exactly ⅓ at ω = 2/3, the fitted work exponent of 0.079 — belongs to the Laplacian on a line. One dimension is where the components were designed, where every closed form is available, and where the arithmetic is cheap enough to run at four grid sizes inside a figure.

It is also the case where multigrid’s advantage is smallest. In one dimension the alternative is a tridiagonal solve, which costs O(n) and is exact; a method whose cost is O(n) and is not exact is not obviously worth building. In two dimensions the direct alternative fills in — the essays on fill and ordering count how much — and the best known elimination order leaves a factor with O(N log N) entries and costs O(N^1.5) to compute. A solver whose cost is proportional to N is then not a constant factor better. It is a different complexity class, and the gap widens with every refinement.

So the second dimension is where the argument actually pays, and this essay measures whether the argument survives getting there. Two of the three headline numbers do not, and the method works anyway — which is the more interesting result and the one that took the phase by surprise.

Three quantities and one knob, and the knob is the anisotropy rather than anything about the method. Running it from one end to the other is what turns three curves into a recommendation:

Both repairs work at the anisotropic end and only one of them keeps working — which is not a small distinction, because a code that measured only at the hard end would conclude that either repair is adequate and that the cheaper one should win.

Crossing is the interesting event and it is not the end of the story: past it the two curves keep separating, and they separate in the direction that punishes the repair.

The two repairs cross at ε ≈ 0.1 and semi-coarsening ends up worse than doing nothing. Point relaxation falls from 0.966 to 0.204 as the anisotropy is removed and semi-coarsening rises from 0.119 to 0.915 — so on an isotropic problem the coarse grid built for an anisotropic one converges at 0.915 a cycle, which is barely converging, against 0.204 for the smoother it was supposed to rescue. A repair applied where it is not needed is not neutral here; it is the worst option on the page.

Line relaxation is never worse than point relaxation, and at ε = 1 it is exactly equal — 0.204 against 0.204, the same number, because with no anisotropy a y-line solve and a point solve are solving the same thing. That is what makes it the safe repair: it degrades to the thing it replaces rather than to something worse.

And there are two optimal weights, not one. Point relaxation wants ω → 1 in the anisotropic limit and the line and semi-coarsening pair want 2/3 — 0.667, 0.667, 0.668, 0.669 at the four strongest anisotropies. Both converge on 0.800 at ε = 1. So the weight a code should use is not a property of the smoother alone; it depends on what the coarse grid is doing, and the two answers differ by half at the hard end.

The operator, and what is exact about it #

The five-point Laplacian on a k×k interior grid, ordered row by row, is the matrix poisson2d built for the sparsity field two phases ago, and the operator this essay measures is written from a 3×3 array of stencil coefficients instead. The two are checked against each other entry for entry and agree exactly, which is the cheapest available guarantee that a new piece of code has not quietly redefined a problem the site has been measuring since its expansion.

What makes this the model problem is that its eigenvectors are known and its eigenvalues are a formula. The grid function sin(mπx)·sin(lπy), sampled at the interior points, is an eigenvector of the operator, and its eigenvalue is

λ(m, l) = 4 sin²(θx/2) + 4 sin²(θy/2)

with θx = mπh and θy = lπh. That is the one-dimensional answer in each direction, added, because the operator is a sum of two one-dimensional operators that commute. Every measured quantity below has that formula sitting beside it, and the agreement is asserted rather than assumed: applying the matrix to a mode and comparing against the scalar its frequencies predict agrees to 10⁻¹⁵.

The condition number follows without an SVD. The extreme eigenvalues are at (1, 1) and (k, k), so κ = cot²(πh/2) — 103.1 at k = 31, growing like h⁻² exactly as it does in one dimension. That is the number every Krylov bound on this site is stated in, and it is worth having in closed form before anything is measured against it.

The rate survives, and it survives in the way that matters #

Three grid sizes, a factor of nearly eighteen in the number of unknowns, and the V-cycle’s asymptotic convergence factor comes out 0.2032, 0.2016 and 0.2037. The spread across the whole range is 0.0028.

That is the same statement the one-dimensional essay made and it is worth more here, because here there is something to beat. At 3,969 unknowns the operator stores 19,593 entries; a Cholesky factor of it under a good ordering stores several times that and the elimination costs O(N^1.5). The V-cycle costs a fixed number of sweeps and transfers per unknown per cycle, and it needs the same number of cycles at every size. The work to a fixed tolerance is proportional to N and to nothing else.

The hierarchy that does it is five levels deep at k = 31: 31×31, then 15×15, 7×7, 3×3, and a single unknown at the bottom, where the “exact solve” is a division. The unknowns over the whole hierarchy come to 1.30 times the finest grid’s, and the stored entries to 1.48 times — the two complexities that the algebraic essays will have cause to distinguish, and which here are both comfortably small.

What the alternative costs, counted #

The claim that a direct solve is the expensive option is easy to assert and this site already has the machinery to count it, so here it is counted on the same operator.

On the 15×15 grid the operator’s lower triangle holds 645 entries. Its Cholesky factor holds 3,389 in the order the grid arrived in and 1,833 under minimum degree — a fill of 1,188 entries that were zero in the matrix and are not zero in the factor, and the better of the two orderings still nearly triples the storage. Neither ordering changes the answer by a digit; that is the whole content of the ordering essay, and it is why the two numbers can be compared at all.

The V-cycle’s hierarchy on the same grid stores 1.48 times the fine operator’s 1,065 entries across all four of its levels — about 1,576 — which is less than the single Cholesky factor of the finest level alone, and it does not need the factorisation that produces it.

That comparison is at one size and it is the ratio’s behaviour rather than its value that decides the argument: the fill grows superlinearly and the hierarchy’s storage does not. What the direct solve buys for its money is an exact answer in a bounded number of operations, which the V-cycle never offers — five cycles for five digits, and a sixth for another. At the sizes this site can draw the two are comparable. At the sizes the method exists for they are not, and the reason is in the exponent rather than in either constant.

The first number that does not survive: the smoothing factor #

In one dimension, weighted Jacobi’s smoothing factor is minimised at ω = 2/3, where it is exactly ⅓ — not approximately, the two branches of the maximum cross there. That number is quoted in the one-dimensional essays as though it were a property of the method.

It is a property of the method in one dimension. Run the same analysis over the square of frequency pairs and the answer at ω = 2/3 is 2/3, and the best available is 3/5 at ω = 4/5.

The reason is geometric and is visible in the figure below. The coarse grid represents the modes with both frequencies below π/2 — a quarter of the square, not a half — so the smoother is responsible for three quarters of it, and the newly added region contains modes that oscillate in one direction and are smooth in the other. Such a mode is damped by only the part of the operator that acts across it, which is half of the operator, so it is damped half as much.

Both numbers have a closed form, and it is the same expression that will do the whole of the anisotropy argument in the next three essays:

μ = max( |1 − ωc/(1+ε)| , |1 − 2ω| )

with c = ε for point relaxation under full coarsening. At ε = 1 that is max(|1 − ω/2|, |1 − 2ω|), the two branches cross at ω = 4/5, and the value there is 3/5. The scan over 96×96 frequency pairs agrees to nine digits.

So the second dimension costs a factor of about two in the smoothing factor and moves the optimal weight. The V-cycle’s convergence factor is 0.20 rather than the one-dimensional 0.10, which is that cost arriving. The claim that the rate does not grow with the size is untouched; the constant in front of it is worse, and it was worth knowing which of those two the second dimension moves.

The second number that does not survive: the Galerkin identity #

The maturity phase asserted this, in a check that still runs on every build:

The Galerkin coarse operator IS the coarse discretisation — R·A·P is the coarse tridiagonal, entry for entry, to rounding, at every level of a six-level hierarchy.

That is true, it is checked, and the natural reading of it is that the identity is what makes the recursion well defined. If the coarse operator were merely approximately the coarse problem, one would expect the correction it returns to be approximately right, and the whole argument to be an argument about how approximately.

The natural reading is wrong. In two dimensions, bilinear interpolation and full weighting applied to the five-point Laplacian give a coarse operator with nine nonzeros in an interior row:

1/16 × ⎡ −1  −2  −1 ⎤
       ⎢ −2  12  −2 ⎥
       ⎣ −1  −2  −1 ⎦

The fine operator has five entries in that row and the coarse one has nine, and the four diagonal couplings are not rounding — they carry a sixth of the row’s total weight. Every level below the first solves an operator the grid it sits on would never have produced.

And the method converges at 0.20 a cycle regardless.

So the identity was never what the recursion needed. What it needs is the property that does survive, and which is checked here in its place: the coarse operator annihilates a constant. The nine stencil coefficients sum to zero to 10⁻¹¹, so a constant error, interpolated up, corrected and restored, comes back unchanged. That is the property that makes a coarse-grid correction a correction rather than a perturbation, and it holds for the nine-point operator exactly as it did for the tridiagonal.

The one-dimensional identity was a coincidence of the one-dimensional operator — a happy one, and one this site was in danger of reporting as a principle.

What this costs to compute, and why the figures are shaped as they are #

A 31×31 grid is 961 unknowns and a dense Galerkin triple product on it is 961³ ≈ 8.9·10⁸ multiply-adds, per level. Done that way, one figure takes minutes and a draggable one takes an hour, and the phase would have consisted of one static picture.

Done from the 4,681 entries the operator actually has, the same product is a few million operations and the figure builds in under a second. So the depth phase carries a sparse matrix type — { n, m, rows }, each row a sorted list of column-value pairs — and the routines that go with it. Nothing about it is a lesson; sparse.js’s header explains why the sparsity field deliberately does not have one, and that reasoning held for six essays and stops holding at exactly the size where the second dimension starts.

The type earns a second keep, which is the algebraic essays’ rather than this one’s: it is what lets the algebraic code be written against rows and columns with no access to anything that knows what a grid is. That claim is checkable by reading an import list, which is worth more than a paragraph promising it.

What is asserted, and what refuses #

Six claims run on every build, and two of them are refusals of statements this site has itself made.

The operator is the old one. aniso2d(k, k, { eps: 1 }) equals poisson2d(k) entry for entry, and every grid mode is an eigenvector with the eigenvalue its frequencies give, to 10⁻¹⁵.

Every smoother’s measured damping is its symbol. One sweep applied to a pure mode multiplies it by a scalar; the scalar is a formula in the two frequencies; they agree to 10⁻¹⁵ for point relaxation and for both line smoothers, over five modes each.

The rate is flat and the anisotropy destroys it. 0.2032, 0.2016, 0.2037 across a factor of eighteen in size — and 0.951 at ε = 10⁻³ with identical components, which is the anisotropic case.

The one-dimensional smoothing factor is refused in two dimensions. The assertion assertClose(μ(ε = 1, ω = 2/3), 1/3) is fed to rejects() and must fail. If a later change ever made it pass, either the analysis or the coarsening would have quietly become one-dimensional.

The one-dimensional Galerkin identity is refused too. assertThat(coarseEntries === 5) is fed to rejects(). The gate insists the coarse operator has nine.

A grid the coarsening cannot halve is refused rather than clamped. hierarchy2d(16) throws. A hierarchy that silently rounded 16 down to 15 would draw a figure of a problem the caption does not describe, which is the failure the standard pass found in hessenberg-reduction and fixed the same way.

What the second dimension does not change #

Worth saying plainly, because the two failures above are the interesting part and it would be easy to read them as the method being shakier than advertised.

The smoother still smooths. Its factor is 3/5 rather than ⅓, which is a worse constant and is still bounded away from one at every grid size — the property the whole method rests on and the one the first multigrid essay established.

The coarse grid still represents what the smoother leaves. That is what the outlined quarter in the frequency square means, and the fact that it is a quarter rather than a half is exactly the factor of two in the smoothing factor.

The recursion still bottoms out in a division, and the work per cycle is still bounded by a constant times the fine grid’s, because the sizes fall by four at each level in two dimensions rather than by two — a faster geometric series, so a smaller bound on the transfer cost than in one dimension.

What changed is that the coarse problem stopped being the same problem, and it turned out that nothing depended on its being so.

The third dimension, which the closed form already predicts #

The coarse grid holds an eighth of the modes rather than a quarter, so the smoothing factor should move once more in a direction the closed form predicts. It does, and the prediction has a clean shape.

Damped Jacobi on the (2d+1)-point Laplacian has symbol 1 − (2ω/d)·Σₖ sin²(θₖ/2). Over the oscillatory region — where some θₖ ≥ π/2 — that sum ranges from ½, one frequency at π/2 and the rest at zero, to d, all of them at π. So the smoothing factor is max(|1 − ω/d|, |1 − 2ω|), the two branches meet at ω = 2d/(2d+1), and the value there is (2d−1)/(2d+1).

Scanning the symbol numerically rather than trusting the algebra:

Each dimension costs the smoothing factor: a third, three fifths, five sevenths, tending to one. A smoother that removes two-thirds of the oscillatory error on a line removes two sevenths of it in a cube.

In sweeps that is the number worth carrying. A tenfold reduction of the oscillatory error takes ln(0.1)/ln(μ) sweeps — 2.10 on a line, 4.51 on a square, 6.84 in a cube — a factor of 3.3 between the first and the last, before any of the third dimension’s other costs are counted. The V-cycle’s flat convergence rate survives the extra dimension; the work per cycle does not, and this is where that goes.

And the one-dimensional weight is increasingly wrong. ω = 2/3 gives 0.658 in two dimensions against the attainable 0.600, and 0.776 in three against 0.714. Carrying the familiar constant across dimensions costs about a tenth of the smoothing factor each time, and a tenth of a smoothing factor compounds over every sweep of every cycle.

assertTheSmoothingFactorIsTwoDMinusOneOverTwoDPlusOne checks the closed form against a scan of the symbol at all three dimensions, checks it is a minimum in ω rather than a crossing, and requires the cube to need several times the line’s sweeps.

What is left #

The anisotropic operator, which a direction the smoother cannot see takes up and where the real difficulty of the field lives. Everything measured here is on the isotropic Laplacian, and the isotropic Laplacian is the case every component was designed for.

The proof. Two-grid convergence for this problem is provable by Fourier analysis and the V-cycle’s by a longer argument, and neither is run here. What is run is the measurement, checked against a closed form wherever one exists — the same standing this site gives the conjugate-gradient bound.

The proof and the anisotropic operator remain — see above.

And the coarse problem as a matrix equation #

On a rectangle with a separable stencil the whole hierarchy is operating on a Kronecker sum — the same problem written as two n×n operators rather than one n²×n² one, with a direct solution in n³.