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_{k−n} 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.

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.

What is left #

T. Chan’s preconditioner, which is the Frobenius-nearest circulant rather than the truncation, and which is positive definite whenever T is — by construction, because it is an average of positive things rather than a truncation of them. So the whole 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”. Measuring the two side by side is a phase’s work rather than a section’s, and this field does not have it.

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.