A preconditioner that changes sign

The previous essay leaves a solver half-specified. A Toeplitz matrix–vector product costs O(n log n) through a circulant embedding, the matrix is symmetric positive definite, and conjugate gradients therefore applies — so the whole question is the step count, and the step count is governed by the spectrum.

The standard answer is to precondition, and the standard preconditioner is the nearest circulant. The argument for it is elegant, the method is genuinely excellent, and this essay is about the measurement that is missing from every summary of it.

The construction, which is a projection #

Strang’s preconditioner is the simplest thing that could work. Take each diagonal of T and wrap it round: diagonal k of the circulant gets tₖ for k in the near half and tₖ₋ₙ in the far half. At even n the entry at exactly n/2 has two diagonals competing for one slot, and Strang’s own definition sets it to zero rather than choosing, which is followed here because choosing would make the operator depend on a convention nobody stated.

Read as a map from Toeplitz matrices to circulants, it is a projection: it keeps as much of each diagonal as fits and discards the corners. There is a cousin — T. Chan’s preconditioner — which is the circulant closest in Frobenius norm, and averages the two ends of each diagonal instead of truncating one. This essay uses Strang’s because its failure mode is sharper and because that failure mode is the finding.

What preconditioning is for here, and it is not the condition number #

Changing the condition number on purpose is the site’s existing essay on preconditioning, and it measures the obvious thing: incomplete Cholesky takes κ from 48.37 to 5.12 on the model problem, and the conjugate gradient count from 35 to 16.

Circulant preconditioning does something different, and the difference is the whole reason it works as well as it does. It does not principally reduce κ. It clusters the spectrum: C⁻¹T has all but a bounded number of its eigenvalues within a shrinking neighbourhood of 1, and the outliers stay outliers.

That distinction is one this site has already established from the other side. The rate the condition number predicts found that at step 20 a κ = 10⁶ problem sits inside the bound belonging to κ = 100 — 1.7·10⁻² against 3.6·10⁻² — because early convergence is governed by clustering and only the asymptotic rate by κ. That was reported there as a surprise about the bound. Here it is a design principle: build the preconditioner to cluster and the bound’s looseness becomes the method.

Eigenvalues, not singular values, and it is not a quibble #

C⁻¹T is not symmetric, so it has both, and they are different sets of numbers. The figure draws eigenvalues, obtained through the similarity to T^½C⁻¹T^½, which is symmetric because T is positive definite.

Singular values would have been easier to compute and would have hidden the entire finding. They are all positive by construction. The eigenvalues are not: at n = 16, ρ = 0.8 three of them are negative, and a plot of singular values at that size would have shown a tidy clustered spectrum next to a step count of 19 against an unpreconditioned 19.

The first version of the generator did exactly that, and the assertion that caught it was not about the spectrum at all — it was the step count refusing to behave. Which is the useful order for this to happen in, and is why the figure carries the sign of the preconditioner’s smallest eigenvalue in its badge.

The measurement #

Strang’s construction takes each diagonal and wraps it. Nothing in that operation preserves positive definiteness, and at moderate n on a strongly correlated matrix it does not preserve it.

The smallest eigenvalue of C, against the size, at ρ = 0.9:

The sign changes between n = 32 and n = 64. Above it the method is spectacular — 5 steps against 133 at n = 256, and the count falls as the problem grows. Below it the method costs steps: 44 against 37 at n = 32, which is a preconditioner making a well-behaved iteration worse.

At ρ = 0.95 the crossing moves right and the damage is larger:

109 steps against an unpreconditioned 66, at a size and a correlation that are both entirely ordinary. And then, one doubling later, 10 against 179.

Why the theory does not mention it #

Because the theory is asymptotic and correct. Strang proved, and Chan and others sharpened, that for a symbol in the Wiener class the preconditioner is eventually positive definite and the spectrum of C⁻¹T eventually clusters at 1 with O(1) outliers, so the step count is eventually O(1).

Every one of those eventuallys is doing work. The published statement is about the limit and there is a sign change on the way to it — at n = 64 for ρ = 0.9 and past n = 128 for ρ = 0.95 — and the crossing size grows with the correlation, which is to say it grows with the difficulty of the problem the method is for.

Swept over ρ, the crossing size and the limit it is approaching are both readable, and the second is a number the previous essay already computed.

Across ρ = 0.5, 0.6, 0.7, 0.8, 0.9 and 0.95 the preconditioner is definite at every size for the first three, crosses between 16 and 32 at 0.8, between 32 and 64 at 0.9 and between 128 and 256 at 0.95. So the crossing size roughly doubles with each step of ρ towards one, and the last step — 0.9 to 0.95 — moves it by a factor of four.

And the number λ_min© is climbing towards is the minimum of the symbol. At n = 256 it reads 0.3333, 0.25, 0.1765, 0.1111, 0.05263 and 0.0242 across the six, against (1 − ρ)/(1 + ρ) = 0.33333, 0.25, 0.17647, 0.11111, 0.05263 and 0.02564. The first five agree to four figures; the sixth has reached 94% of its limit and is the one whose crossing is still in the range drawn.

That closes the loop with the previous essay, which is about a finite section approaching a symbol’s extremes from inside. The same thing is happening here to a different matrix, and this is what it costs: the symbol’s minimum is positive at every ρ below one, the finite sections approach it from below, and the part of that approach that lies below zero is the region where a method justified by the limit does damage. A sign change is what approaching a small positive number from below looks like — so the failure is not a defect of Strang’s construction, it is the previous essay’s asymptotic gap with a sign attached.

This is the same shape as the previous essay’s finding and the same shape as nested dissection losing to minimum degree at every size this site draws. A statement about a limit is being read as a statement about a computation, and the gap is not a rounding detail — here it is the difference between 10 steps and 134.

What the assertions say #

Three claims, each in the direction it holds, because a single unconditional one would assert something false over half its range. That is the discipline the mixed precision field established at its own threshold and the multigrid field at its own.

The sign change happens once. λ_min© rises monotonically with n and crosses zero exactly once in the range drawn, which is asserted step by step and then as a whole. An oscillating sign would be a different phenomenon and would need a different explanation.

Where C is definite, more size does not cost more steps. Not the count is constant, which is what the summary says and is not what happens at these sizes: at ρ = 0.9 it goes 28, 9, 5 as n doubles twice, still falling. The checkable claim over this range is the direction, against a plain count asserted to be climbing on the same systems.

And where C is indefinite, the repair is not an improvement. This is the assertion the phase was not expecting to write, and it is stated as an assertion rather than as a remark for the reason the multigrid field gives about a coarse-grid correction that made an error worse: an outcome that contradicts the published account is worth a check of its own, because a check is the only thing that will still be true after a refactor.

What clustering measures, and what it does not #

The figure’s badge carries the fraction of the preconditioned spectrum within a tenth of 1, and it is worth being clear about what that number is and is not.

It rises with n — 12.5% at n = 16, 31% at n = 32, 97% at n = 64 for ρ = 0.8 — which is the clustering theorem’s content. And it is not the same claim as definiteness. At n = 32, ρ = 0.8 the preconditioner is already positive definite, λ_min = +0.0798, and only 31% of the spectrum is clustered. The first version of this generator asserted the two were equivalent and was wrong.

Definiteness is when the method stops being harmful. Clustering is when it starts being spectacular. There are sizes in between where it is neither, and those are the sizes at which somebody would first try it.

Where the negative eigenvalue comes from #

It is not mysterious, and seeing where it comes from is what makes the size dependence obvious.

The circulant’s eigenvalues are the transform of its first column, and Strang’s first column is the Toeplitz one truncated at n/2 and wrapped. For the ρ^|i−j| family that column is

1, ρ, ρ², …, ρ^(n/2−1), 0, ρ^(n/2−1), …, ρ², ρ

— the geometric sequence, cut off, and mirrored. Its transform is the geometric sum with the tail removed, and removing the tail of a geometric series is subtracting a term that is small in magnitude and oscillating in sign across frequencies. At the frequency where the truncation error is most negative, the eigenvalue it perturbs is the smallest one — which is the one nearest zero and has the least room.

So two things race. The truncated tail is of size ρ^(n/2), which falls fast with n and slowly with ρ. The smallest eigenvalue it has to perturb is heading for (1 − ρ)/(1 + ρ), which shrinks with ρ. The sign change happens when the first becomes smaller than the second, and

ρ^(n/2)  ≈  (1 − ρ)/(1 + ρ)

solves to n ≈ 56 at ρ = 0.9 and n ≈ 143 at ρ = 0.95, against measured crossings between 32 and 64, and between 128 and 256. Not a fit and not a proof — a two-line estimate with no free parameter in it that lands inside the measured interval at both correlations, which is enough to say the mechanism is understood rather than merely observed.

The estimate also predicts the shape of the dependence, which is the part a table of five numbers cannot. Solving it for n gives

n  ≈  2 · log((1 − ρ)/(1 + ρ)) / log ρ

so the crossing grows without bound as ρ approaches one, and grows like 1/(1 − ρ) there. A family with a correlation of 0.99 does not cross until about n = 1,053. That is not a size nobody runs — it is a perfectly ordinary Toeplitz system — and it is a size at which the unpreconditioned iteration is genuinely painful, which is to say the case where the method is most wanted is the case where it is furthest from working.

It also says immediately why T. Chan’s variant does not have the problem. Averaging the two ends of each diagonal rather than truncating one produces a first column that is a non-negative average of the Toeplitz entries, and a non-negative combination of positive-definite things is positive definite. The two constructions differ by which of two obvious things to do with the corner, and one of them keeps a property the other throws away.

What this costs a reader who tries it #

The practical shape of this is worth stating, because it is not “avoid circulant preconditioning”.

Somebody meeting a Toeplitz system for the first time writes the unpreconditioned iteration, finds it slow, reads that circulant preconditioning fixes it, implements Strang’s — which is the one presented first, because it is the one with the one-line description — and measures. On a small test problem they get 44 steps where they had 37. The natural conclusion is that the implementation is wrong, because the method is famous and the measurement is a factor of 1.2 in the wrong direction, which is exactly the size of a plausible bug.

There is no bug. The measurement is right, the method is right, and the test problem is on the wrong side of a threshold nothing told them about. Scaling the test problem up by a factor of eight — which is the last thing anybody does while debugging — turns 44 into 5.

That is the sequence this essay exists to interrupt, and it is the same sequence the maturity phase’s sparse LU essay describes from the other direction: an unpivoted factorisation that reproduces its matrix to 3.8·10⁻¹⁷, better than every pivoted run, and returns an answer wrong in the fifth digit. In both cases the measurement a reader would naturally take is the one that misleads, and the one that settles it is a measurement nobody would think to take.

The badge, on a figure with no factorisation in it #

residualcheck does not require this essay’s figures to print a residual — neither of them draws a decomposition — and both print one anyway, because the interesting quantities here are of the same kind.

strang-sign-change prints λ_min© at every size on the sweep. That is not a residual and it is the number the whole figure turns on, and putting it in the badge rather than in the caption means it is generated rather than typed. clustered-spectrum prints the clustered fraction, the count of negative eigenvalues and λ_min© together, which is the pair of facts the section above spends four paragraphs separating.

The rule is no decomposition without its residual, and the habit underneath it is that the number a figure rests on belongs on the figure. Where the two come apart — a figure with no factorisation and a load-bearing number — the habit is the one worth keeping.

Two numbers this field now has that the iterative field did not #

Worth recording, because they are the reason the structure was worth a field of its own rather than three essays appended to the iterative one.

A step count with a ceiling. Every convergence claim in the iterative field is about a model problem whose condition number grows with the grid, so every step count there grows too and the whole field is about slowing that growth down. Here κ is bounded above by Szegő’s limit at every size, so the unpreconditioned count plateaus on its own — 34, 34 at ρ = 0.5 — and the preconditioner is competing against something that already stops rather than against something that runs away.

And a preconditioner with a closed form. Incomplete Cholesky is defined by an algorithm: run the factorisation, drop what was not there before. Its own residual ‖A − LLᵀ‖/‖A‖ is 0.0825 on the model problem and is asserted large, which is the one place on this site where a residual is asserted in that direction. Strang’s circulant is defined by a formula on the entries, so its spectrum is a transform away and every claim about it — including its sign — is available before it is applied to anything. That is why the sign change is measurable at all, and it is why an equivalent statement about incomplete Cholesky would be much harder to make.

The other preconditioner, on the whole grid #

T. Chan’s preconditioner is the Frobenius-nearest circulant rather than the truncation, and it is positive definite whenever T is — by construction, because it is an average of positive things rather than a truncation of them. So the finding above is specific to Strang’s variant, and the honest form of the practical advice is not “circulant preconditioning has a sign problem” but “one of the two standard circulant preconditioners does, and it is the one usually presented first”.

That much can be said from the construction. What the construction does not say is what the choice costs, and the number usually offered for it is measured at one point: at ρ = 0.95 and n = 256 both preconditioners are definite and the counts are nine steps against ten, from which “the safer one costs a step” is a fair reading of that row and a poor summary. Sixteen combinations, same tolerance, iteration cap 400:

Chan’s column is between 7 and 10 at every one of the sixteen. Four steps of total range, across a family whose unpreconditioned count runs from 29 to 204 — it does not move with the size and it does not move with the correlation. That is what the clustering theorem promises, delivered without a caveat, and it is worth pausing on because the theorem is usually stated for both variants and only one of them keeps it here.

Strang’s column runs from 3 to 400, and the 400 is the iteration cap rather than a convergence. At ρ = 0.99 and n = 256 the wrapped construction has not converged in twice the unpreconditioned count. Every row on which it loses has λ_min < 0, without exception in the grid.

So the trade is not symmetric, and “one step” understates both halves of it. Where Strang wins it wins by a factor of two — three steps against seven at ρ = 0.7 — not by a step. Where it loses it loses by a factor of forty. Trading a possible factor of two for a possible factor of forty is not a trade, and the practical advice that follows is not the balanced one the single row suggests: use the averaged circulant.

And the sign is not the criterion, which is the part the grid shows and no single row could. The row at ρ = 0.9, n = 64 has λ_min = +0.017 — positive, definite, inside the theorem’s hypotheses — and takes 26 steps against Chan’s 10. Two rows below it λ_min is +0.051 and the count is 7. The failure does not switch off at zero; it fades out over a neighbourhood of zero, because a preconditioner with an eigenvalue barely above zero divides by a small number just as surely as a negative one divides by a wrong-signed one. A code that guarded Strang’s construction with min(eig(C)) > 0 would have accepted that row and spent two and a half times the steps.

The last row is also the fitted threshold behaving exactly as the size the rank does not notice would lead one to expect a threshold not to: n★ ≈ 1056 at ρ = 0.99, so no size on this grid rescues the wrapped construction there, and the crossing this essay’s title is about is simply off the right-hand edge.

assertTheGridSaysUseTheAveragedOne measures all sixteen rows, requires Chan’s count to stay single-figure and Strang’s to span two orders, requires every indefinite row to lose, and requires at least one definite row to lose badly — because that last one is the claim a sign test would get wrong.

What is left #

The non-symmetric case. Everything here rests on T being positive definite, which is what makes C⁻¹T’s eigenvalues real and makes conjugate gradients applicable at all. A non-symmetric Toeplitz system needs GMRES, and the spectrum predicts nothing there — this site has a cyclic shift with every eigenvalue on the unit circle on which GMRES makes no progress for n − 1 steps. A clustering argument in that setting is not obviously an argument at all.

And the two-dimensional case, where the operator is block-Toeplitz with Toeplitz blocks and the natural preconditioner is block-circulant with circulant blocks. The clustering theorem is weaker there in a way that is known to be essential rather than technical, and it is the case every image deblurring problem actually is — including the one the regularisation field builds, which uses a one-dimensional blur precisely to avoid needing any of this.

A preconditioner nobody can build #

This page’s preconditioner is built from the matrix’s structure. Where the operator is a subroutine there is no structure to read, and the choice narrows to what can be described without looking.