Two dimensions, and the cluster that thins

This is the third time this site has taken a one-dimensional result into two dimensions, and the third time the load-bearing part has been the part that did not survive.

The multigrid field lost the Galerkin identity: in one dimension the coarse operator RAP is the coarse discretisation entry for entry, and in two a five-point operator produces a nine-point coarse one. The smoothing analysis lost ω = ⅔ and μ = ⅓, which turned out to be 4/5 and 3/5 in the plane. Now the structure field loses the cluster.

What makes each of these worth writing rather than noting is that the method survives every time. The two-dimensional multigrid works; the coarse operator is simply a different one. The circulant preconditioner below works too. What stops being true is the sentence people repeat about it.

What a BTTB matrix is, and why this one #

Discretise a convolution on an mx × my grid and order the unknowns row by row. The matrix is block Toeplitz with Toeplitz blocks: constant along block diagonals, and each block constant along its own diagonals. It has n = mx·my rows and is described by (2mx − 1)(2my − 1) numbers, so the compression is the same kind the field’s first essay was built on, and the same trade a hierarchy makes with its blocks: a large object described by a small number of numbers, with the description rather than the entries deciding what is affordable.

The kernel drawn here is ρ^(|k₁| + |k₂|) — the Kac–Murdock–Szegő family’s outer product with itself. That choice is what makes every comparison below honest: the two-dimensional problem is the one-dimensional problem twice, so anything that changes between them is a property of the dimension rather than of a different matrix. A reader suspicious that the second dimension has been handed a harder problem can check the kernel and see that it has not.

The symbol multiplies, so the difficulty squares #

Sum the double series and it factorises: f₂(θ₁, θ₂) = f(θ₁)·f(θ₂), with f the one-dimensional symbol this site already has. Checked at three arbitrary angle pairs to 10⁻¹⁴, because a claim of the form “it factorises” is exactly the kind that is true of the formula somebody wrote down and false of the matrix somebody assembled.

The range of a product is the product of the ranges, so the limiting condition number is

κ → ((1+ρ)/(1−ρ))⁴

the fourth power rather than the second — 6,561 at ρ = 0.8 where the one-dimensional family’s limit is 81.

Measured against that limit, κ of the finite grids is 581, 1196, 1794 and 2337 at grid sides 4 to 10, reaching 8.9%, 18%, 27% and 36% of 6,561. The breadth phase’s finding survives into two dimensions and gets worse: the one-dimensional 8×8 section reached 52% of its limit, and the 8×8 grid — sixty-four unknowns, the same count — reaches 27% of its.

So the structure is intact and the problem is harder, and both halves are what the symbol says.

What the decay rate changes, and what it does not #

Everything above is drawn at one ρ. The slider carries five, and running it separates two quantities that the single-ρ reading leaves joined: how hard the problem is, and how much the preconditioner is worth on it.

The unpreconditioned column has moved by half and the preconditioned one has barely moved at all. That continues to the end of the slider.

Reading the largest grid down the slider:

The preconditioner is worth 2.4 times at ρ = 0.5 and 7.5 times at ρ = 0.95, and the reason is that only one of the two columns is a function of ρ. The unpreconditioned count triples across the slider; the preconditioned count reads 18, 21, 21, 21, 19 — inside the noise of an integer step count.

That is worth having beside this essay’s negative result, because the two are easily confused. What fails in two dimensions is the clustering theorem — the share of the spectrum near one falls from 56% to 17%, and the step count stops being independent of the size. What does not fail is the preconditioner’s usefulness, which improves as the problem gets harder, exactly as it does in one dimension. A method can lose the property it was sold on and remain the right thing to use.

The preconditioner generalises exactly #

Averaging along each axis in turn gives a matrix that is block circulant with circulant blocks, and everything the one-dimensional construction had comes with it.

It is diagonalised by the two-dimensional transform, exactly. Two routes to one set of numbers — the transform of the first column, and a symmetric eigensolve of the assembled matrix — agree to 4.4·10⁻¹⁴ at three grid sizes, and solving with it undoes multiplying by it to 1.3·10⁻¹⁴. That is what licenses the preconditioner solve to be three lines: transform, divide, transform back.

It is still the minimiser, and that is checked against separability rather than assuming it. The entries of the averaged BCCB are compared against the plain entrywise average of A over each residue class of the index pair — a different computation, agreeing to 3.3·10⁻¹⁶ — because separability is precisely the property a reader would be right to suspect of doing all the work.

And it is still positive definite, for the same one-line reason: every entry is a convex combination of entries of A’s own diagonals — the circulant that cannot be indefinite, whose argument never mentions how many axes the index runs along.

The two-dimensional transform is dft rather than fft here, deliberately. The grids are 4 to 10 a side, most of which are not powers of two, and fft refuses a length that is not one rather than padding silently — because padding computes the transform of a larger grid and calls it this one’s, which is the exact answer to a nearby problem with nobody told how near. The direct transform is O(m²) a line and invisible at these sizes.

And then the cluster thins #

The one-dimensional clustering theorem says all but O(1) of the preconditioned eigenvalues fall into a shrinking neighbourhood of one. Its two-dimensional analogue says all but O(mx + my) do — and at n = mx·my that is a vanishing fraction and a growing count.

Both halves are visible, and which one appears depends on how wide a neighbourhood is asked about.

Asked about the tighter band the one-dimensional essays use — within a tenth of one — the count is 5, 1, 0, 0: the cluster is empty at the two larger grids, which is why both windows are reported rather than the flattering one. A figure drawn only at the tight band would show nothing and suggest the preconditioner had failed; a figure drawn only at the wide one would show a cluster growing and suggest it was fine.

What conjugate gradients pays for is the count outside, not the share inside. Each eigenvalue away from the cluster costs at most one extra step — the rate the condition number predicts is a bound on a polynomial, and an outlier is what a polynomial has to spend a root on — which is why the preconditioned count climbs — 10, 18, 20, 21 — while the one-dimensional count at the same unknown counts does not move.

The comparison is the finding #

Without the one-dimensional column, the climb above is a statement about sizes rather than about dimensions. With it:

Same kernel, same averaging, same transform, same conjugate gradient code, same right-hand side construction, same unknown counts. The only difference is how many axes the indices run along.

And the preconditioner is still very much worth having: 21 steps against 115 at a hundred unknowns. The claim that fails is not “the preconditioner works”. It is “the step count stops depending on the size”, which is what the method is sold on and what makes an O(n log n) solver an O(n log n) solver rather than an O(n^1.5 log n) one.

At ρ = 0.5 and the smallest grid, the preconditioned and unpreconditioned counts are both 10 — sixteen unknowns at that correlation are solved before either method has anything to distinguish it. The assertion behind the figure was written as helps at every grid, and that slider stop failed it; it now says does not hurt at any grid, and helps on at least one, which is what the measurement supports. That was found by the gate that renders every stop of every slider, which is the second time this phase that a range has been wrong in a way no essay could reach.

Why it thins, in one paragraph #

The one-dimensional argument goes: T − C is small in a norm and low rank apart from a small perturbation, so C⁻¹T is the identity plus a small-plus-low-rank matrix, and low rank means a bounded number of eigenvalues can escape the cluster.

The seam that the difference lives on is the boundary between the wrapped and unwrapped parts — a distinction of the same kind as the one the kernel with nothing to compress turns on — and in one dimension that boundary is two points — the ends of the vector. In two dimensions it is the edge of the grid: 4m points on an m×m grid, growing with the side. The rank of the difference grows with the boundary, and the boundary of a square is one-dimensional, so the count of escapees is O(m) and the count of unknowns is O(m²).

That is the same shape as the multigrid field’s loss, and worth naming: the identity in one dimension was a boundary effect nobody had to count, and in two dimensions the boundary is large enough to count.

How fast it actually grows #

The paragraph above says the escapee count is O(m), and the summary sentence at the end of this essay turns that into a claim about the step count: it grows like the side of the grid. The first is an argument about rank and it is right. The second is a step from it, and it is the kind of step this site is supposed to measure rather than take.

The measurement is cheap, because a step count needs no eigenvalues — the conjugate gradient loop is a dense matrix–vector product and nothing else, so it runs far past the 10×10 the spectrum figures stop at. Four times the side, sixteen times the unknowns:

m       n      unpreconditioned    preconditioned    if linear in m
 8      64           73                  20                20
12     144          167                  26                30
16     256          250                  28                40
20     400          333                  30                50
24     576          437                  30                60
28     784          467                  32                70
32   1,024          600                  31                80

Four times the side buys 1.55 times the step count. Linear growth would have predicted 80 and the measurement is 31, so the sentence overpredicts by a factor of 2.6 at the far end of this range and by more beyond it. Fitted as a power of m, the exponent is 0.31; at ρ = 0.8 it is 0.17, and the counts there are 22, 26, 28, 29, 28, 28, 29 — almost flat over the same range. Fitted as a + b·log m, every point on both curves lands within 1.7 steps.

So the growth is real, and it is logarithmic in the side rather than proportional to it. In the unknowns that is log√n, which is to say the preconditioned solve is O(n log n · log n) rather than the O(n^1.5 log n) the linear reading implies — a difference between a method that is essentially optimal and one that is not.

Where the step from the rank argument goes wrong is the phrase each eigenvalue away from the cluster costs roughly one extra step. That is an upper bound, and conjugate gradients beats it whenever the outliers are themselves grouped: a polynomial that is small on a tight bunch of outliers costs one root for the bunch rather than one per eigenvalue. The escapees here come from one mechanism at one scale — the wrap at the boundary of the grid — so they are not m independent perturbations, and counting them as though they were is what produces the linear estimate.

That does not rescue the one-dimensional sentence, and it is worth saying which claim is left standing. The count is not independent of the size in two dimensions: it moves from 20 to 31 over this range and it is still moving. What it is not is proportional to the side, and the difference between “grows” and “grows like m” is the difference between a caveat and a change of complexity class.

The unpreconditioned column is the control that says the measurement is capable of seeing growth. It runs 73 to 600 over the same grids, a fitted exponent of 1.46 in the side — so the instrument reports a steep climb where there is one, and reports a logarithm here.

What the transform had to be checked for #

Three things in this file are exact claims about a computation, and each is checked by a route that does not share arithmetic with the thing it checks — which is the site’s two-routes-to-a-number habit applied where the temptation to skip it is strongest, because everything here is “obviously” the one-dimensional case twice.

The transform diagonalises the BCCB matrix. The eigenvalues from the two-dimensional transform of the first column, against a symmetric eigensolve of the assembled matrix: agreeing to 4.4·10⁻¹⁴ at three grid sizes. Without that, an index error in the row-major flattening would produce a set of plausible positive numbers and a preconditioner that solves with the wrong matrix.

Solving undoes multiplying. A vector multiplied by the assembled BCCB and then solved back through the transform returns to itself to 1.3·10⁻¹⁴, which is the check that the divide-by-eigenvalue step has its conjugations the right way round.

And the averaged BCCB is the entrywise average of A. The construction applies the one-dimensional weighting along each axis in turn, which is the minimiser if the separability argument holds. Rather than assume it, the entries are compared against the plain average of A over each residue class of the index pair — a different computation with no weights in it — agreeing to 3.3·10⁻¹⁶.

The transform itself is dft rather than fft, and that is a decision rather than laziness. These grids are 4 to 10 a side and most of those are not powers of two; fft refuses a length that is not one rather than padding, because padding computes the transform of a longer signal and returns it as this one’s. The direct transform costs O(m²) a line, which at these sizes is invisible, and the alternative would have been an embedding whose correctness is exactly the kind of thing this section exists to check.

What this essay is not entitled to say #

The grids stop at 10×10 because the preconditioned spectrum is computed by a dense symmetric eigensolve of a 100×100 matrix, and the figure needs every eigenvalue rather than a few.

The step counts alone would run much further — the conjugate gradient loop is a dense matrix–vector product per step and nothing else — and a claim about how the cluster behaves at 40×40 is not one these measurements support. What they do support is a trend across four grids in both quantities at once, with the count inside growing like the side and the count outside like the area, which is the shape the theorem predicts and is checkable at this size.

There is also a boundary of kind rather than of size. Everything here uses a separable kernel, chosen so the two-dimensional problem is provably the one-dimensional problem twice. A genuinely two-dimensional kernel — a rotated Gaussian blur, say — has a symbol that does not factorise, and every closed-form statement in this essay would become a measurement. Whether the cluster thins the same way there is an open question this essay does not touch, and the honest reading of it is as the best case rather than the general one.

What survives, and how much of it #

Three things do, and it is worth being precise, because “the theorem is weaker” is not a useful summary on its own.

The preconditioner is definite in every dimension, for a reason that does not know how many axes there are: an average of positive numbers is positive. That is the previous essay’s one-line argument, and nothing in it mentions the index set.

The solve is still O(n log n), because the two-dimensional transform is, and a preconditioned step still costs an unpreconditioned one plus a transform. The method’s cost per step is unchanged; what changed is how many steps there are.

And the clustering still happens. The cluster is real, it grows, and the spectrum stays bounded away from zero — the smallest eigenvalue of C⁻¹A is 0.27 at every grid drawn, against an upper end that grows from 3.6 to 16.4. What changed is the growth rate of what is outside the cluster, and that is the quantity nobody quotes.

So the correct summary of this essay is not that circulant preconditioning fails in two dimensions. It is that a step count sold as independent of the size is independent of the size in one dimension and is not in two — measured, growing like the logarithm of the side rather than like the side, which is still a very good method and is not the method the sentence describes.

What is left #

Non-separable kernels. Everything here uses ρ^(|k₁|+|k₂|), whose symbol factorises, and the factorisation was used to state the limit in closed form. A genuinely two-dimensional kernel — a Gaussian blur with a rotation in it, say — has a symbol that does not factorise, and every closed-form statement above would have to be replaced by a measurement.

The better two-dimensional preconditioners, which exist precisely because of what this essay measures: the block-circulant-with-circulant-block family is the naive lift, and the practical answers are block preconditioners built level by level rather than by averaging in both directions at once.

And the sizes. The grids here go to 10×10 because the preconditioned spectrum is computed by a dense symmetric eigensolve, which is a hundred-by-hundred problem at the largest grid drawn. The step counts alone would run much further; the spectrum is what limits the figure, and a claim about how the cluster behaves at 40×40 is not one this essay is entitled to make.