What a cheap preconditioner has to leave alone

A preconditioner that arrives past the answer measured the one preconditioner for a stopped conjugate gradient run whose effect can be written down exactly. P = AᵀA + αI, applied through the singular value decomposition of the blur, makes the first preconditioned step a Tikhonov solution at $\lambda^2 = \alpha$. On the 64-point deconvolution at 1% noise that gives it an optimum: at α = 10⁻³ the best iterate is the first, at 0.1412, and below the oracle’s $\lambda^2$ every shift is worse than the last, without limit. The shift turned out to be a second regularisation parameter.

Nobody runs that preconditioner. Applying it costs the decomposition that would have solved the problem outright, and it was chosen because its arithmetic could be followed to the last digit rather than because anybody could afford it. What a deconvolution of any real size uses instead is an approximation of the blur that a fast transform diagonalises — the Fourier transform when the signal is treated as periodic, the cosine transform when it is reflected at its ends — so that applying the preconditioner costs a transform, a division and a transform back. The matrix that is one row is the account of why a circulant hands over its whole spectrum for the price of one transform, and why that makes it the natural preconditioner for anything that is nearly one.

Such an approximation does not commute with AᵀA. Its eigenvectors are cosines rather than the blur’s singular vectors, so its first step is not a Tikhonov solution of anything, and no filter factor describes what it does to a run. The exact essay ended on exactly this case and left two questions: whether a preconditioner of that kind still has a best parameter on a stopped iteration, and whether that parameter still sits where the oracle’s regularisation parameter does.

It does, and only if the approximation is made invertible in a way the exact preconditioner never needed. Shifted by αI, as before, neither of two standard approximations beats the unpreconditioned run at any shift. Truncated instead — every eigenvalue of the approximation below a level τ passed through rather than divided — both reach the plain run’s floor in a quarter of its steps, and the smallest τ that keeps the floor sits at a third to a half of the oracle’s λ, at three noise levels two decades apart.

Two approximations that disagree with the blur only at its ends #

The blur is a Gaussian of width 2.5 grid points with each row normalised to sum to one, so that the operator keeps a constant signal constant. In the middle of the interval a row is the kernel itself. In the first and last few rows the kernel runs off the end of the signal, and the normalisation rescales what is left onto the points that remain — a choice about the boundary that is natural to write down and that no fast transform diagonalises.

The two standard approximations handle the ends differently, and neither handles them the way A does. The Fourier one treats the signal as periodic: the part of the kernel that runs off the right end reappears at the left, so its matrix is a circulant and the Fourier transform diagonalises it exactly. The cosine one reflects the signal about each end: the part of the kernel that runs off folds back onto the points it came from, and the cosine transform diagonalises the result. Both keep a constant constant, as A does, and both agree with A to within 10⁻¹⁰ in every row at least sixteen points from an end, where the only difference left is the far tail of a reflected or wrapped kernel.

The picture shows a disagreement that is local, and in the norm that matters to a solver it is not small. The Fourier approximation is 0.497 of ‖A‖ away from the blur — the periodic wrap couples the first point to the last with a weight A does not have at all — and the cosine approximation 0.117 away. Both differences have rank 22 once singular values below 10⁻¹⁰ of ‖A‖ are discarded, which is to say that a perturbation confined to fourteen rows spreads across a third of the dimensions of the problem. A limit the matrix never reaches makes the structural version of the same point about Toeplitz families: a statement about the infinite family is not a statement about the finite section, and the ends of the section are where the two part company.

For a well-posed system this is a standard trade and a good one. The circulant that cannot be indefinite measured it, and a circulant that differs from its matrix in a corner of low rank still brings conjugate gradients down to seven or ten steps, because a low-rank disagreement adds a handful of outlying eigenvalues and the iteration spends a step on each. A stopped iteration on an ill-posed problem cannot read the same disagreement the same way. The directions a stopped run is supposed to leave untouched are the directions in which the approximation and the blur disagree most, and what the preconditioner does to those directions is the whole of the question.

The shift that served the exact preconditioner #

The first construction to try is the one that worked. The exact preconditioner was (AᵀA + αI)⁻¹. Its cheap copy is (MᵀM + αI)⁻¹, with M the approximation, applied as a transform, a division of the k-th coefficient by $\mu_k^2 + \alpha$, and a transform back — where $\mu_k$ is the approximation’s k-th eigenvalue.

At α = 10⁻³ the exact preconditioner delivered its best answer in a single step. The copy built from the cosine approximation is least at step 16, at 0.749, more than five times the floor the plain run reaches with no preconditioner at all. And that is not an unlucky shift. Across the whole sweep, neither approximation does better than leaving the iteration alone:

The best entry in either approximation’s column is at the largest shift, where the preconditioner hardly preconditions: at α = 10⁻¹ it saves eight of the plain run’s twenty steps and loses a per cent and a half of accuracy doing it. Every smaller shift costs more, and the cost is not bounded by what the exact preconditioner paid. At α = 10⁻⁴, where the exact one’s first step was still within 13% of the best Tikhonov answer, the cosine copy’s best iterate is 3.36 — more than three times the size of the signal it was meant to recover.

The Fourier column has a second shape worth reading. From α = 10⁻⁴ down its best iterate is the first, at an error just above one, which is the error of an answer that has barely left zero. The first residual in those runs is a hundred times the noise, as large as the data itself: the preconditioned search direction points almost entirely into directions where the approximation is wrong, conjugate gradients take a tiny step along it, and every later step moves the iterate further from the truth. The run has not fitted any of the signal and is already past its best.

Dividing by what it gets wrong #

The mechanism is in the division. The shifted copy scales the k-th cosine coefficient by $1/(\mu_k^2 + \alpha)$. For the leading cosines $\mu_k$ is close to the blur’s own singular values and the scaling does what the exact preconditioner did. For the high cosines $\mu_k$ falls as far as 1.2·10⁻¹³, and the scaling approaches 1/α: ten thousand at α = 10⁻⁴, a million at 10⁻⁶. On the approximation itself that is harmless. Its high cosines are directions it maps to nearly nothing, and dividing nearly nothing by a small number is a preconditioner doing exactly its job.

On A it is not harmless, because A does not map those cosines to nearly nothing. Near the ends A and M differ, and the difference is not smooth: it is a renormalised kernel against a folded or wrapped one, a jump in the row weights at the boundary, and a jump is made of high frequencies. So AᵀA applied to a high cosine carries a component, localised at the ends, whose size is set by the disagreement rather than by the blur’s spectrum, and the preconditioner multiplies that component by 1/α. The noise in the data has components along the same directions. A stopped run preconditioned this way builds its first steps out of amplified boundary disagreement and amplified noise, and the table is what that looks like from outside.

This is the reverse of what kept the exact shifted preconditioner safe down to $\lambda^2$. There, the directions divided by α were the blur’s own small singular vectors, the division was a Tikhonov filter, and the iteration refined within it. Here the divided directions belong to the approximation, and the division acts on an operator that does not share them. The basis decides what a filter is found the same fragility in the vocabulary itself: an 8.6% asymmetry in the blur was enough to take GMRES’s weights off any function of σ. A preconditioner’s eigenvectors are a basis too, and a division performed in the wrong basis is not a filter at all.

Two spectra with the same cluster #

The ordinary way to judge a preconditioner is to look at its spectrum, and here the ordinary look does not see the difference between a construction that works and one that does not.

Shifted at α = 10⁻⁴ or truncated at τ = 10⁻² — the same threshold, stated in the units of a singular value — the cosine approximation gathers about the same number of eigenvalues near one: twenty-three for the shift and twenty-one for truncation. By the measure the rate the condition number predicts and the essays after it read a preconditioner on, these are the same preconditioner, and the shift is marginally the better one.

The spectra differ in their two tails. At the top the shift pushes its two outliers to 10.9 against truncation’s 5.58; at τ = 10⁻³ the gap is 320 against 112, and for the Fourier approximation at τ = 10⁻² it is 645 against 76. At the bottom the shift leaves 30 eigenvalues below 10⁻³ and truncation leaves 39. Those nine directions are the ones the shift lifted off the floor by dividing them, and a lifted direction is one conjugate gradients will reach within a few steps. On a well-posed system reaching every direction quickly is the point of preconditioning. On this one those nine directions are made of noise and boundary error, and the directions left near zero are the ones a stopped run never reaches — which is where a stopped run’s regularisation comes from.

So the cluster near one, the quantity a preconditioner is built and advertised to produce, is the one part of the spectrum that says nothing here. Changing the condition number on purpose describes preconditioning as choosing a spectrum rather than inheriting one; on a stopped iteration the part worth choosing is the bottom of it.

Leaving the small eigenvalues alone #

The repair keeps the approximation and changes how it is inverted. Where $|\mu_k| \ge \tau$ the coefficient is divided by $\mu_k^2$, as before. Where $|\mu_k| < \tau$ it is multiplied by one — not divided, not shifted, simply passed through as if there were no preconditioner at all. The preconditioner then accelerates the directions the approximation can be trusted on and leaves the rest to the plain iteration’s slow pace, which on an ill-posed problem is the pace that regularises. It is the construction the literature on preconditioned regularisation adopted for these approximations, and the measurement shows why nothing gentler will do.

Truncated, both approximations reach the floor. The cosine run’s best is 0.1412, at step 5, with τ = 10⁻²; the Fourier run’s is 0.1390, also at step 5, at the same τ. Both are a shade below the plain run’s 0.1426 at step 20, and the Fourier one is below the best Tikhonov solution’s 0.1405 — a margin of the size four knobs and one floor found separating regularisers that share no arithmetic, and not evidence of anything beyond that. From τ = 0.3 down to τ = 10⁻² the cosine run’s best stays within 2% of the floor while its best step falls from fifteen to five.

Below 10⁻² the floor is lost. Half a decade further down the cosine run’s best is 0.234, and at τ = 10⁻³ it is 0.484. The failure is the one the shift showed at every level, arriving late: at a small enough τ the directions being divided include ones the approximation gets wrong, and dividing them lets the noise and the boundary disagreement in. What truncation changes is not whether that happens but where. With a shift it begins at the first α that divides anything the approximation is wrong about, which is every α small enough to precondition at all. With truncation it begins at the τ below which the approximation’s trusted directions stop being the blur’s.

Four times further from the blur, and no worse for it #

The Fourier approximation is the sharper test of that reading, because it is four times further from A than the cosine one — 0.497 of ‖A‖ against 0.117 — and the whole of the disagreement is in a place the shift divides.

Shifted, the Fourier approximation is the worse of the two at every small shift, and its outliers are enormous by comparison: 645 against the cosine copy’s 10.9 at the same threshold, and 46,235 against 320 at τ = 10⁻³. Truncated, it is the better of the two at 1% noise and slightly the worse at 0.1%, where it reaches 0.1050 at step 17 against the cosine run’s 0.1039 at step 8. A factor of four in the size of the disagreement with A barely shows in the truncated runs, and dominates the shifted ones. What decides how the preconditioner behaves is not how wrong the approximation is but whether the preconditioner divides by the part that is wrong.

That is worth knowing about a family of constructions whose whole appeal is that their disagreement with the operator is not zero and cannot be made zero at the price. It is also a sharper form of the robustness two dimensions and the cluster that thins and its neighbours measured for circulant preconditioners on well-posed systems, where a boundary disagreement costs a few outlying steps. On a stopped run it costs nothing if it is passed through, and every digit if it is divided.

Where the floor stops being kept #

A shift for the exact preconditioner was a regularisation parameter whose optimum sat at the oracle’s $\lambda^2$. The truncation level is a parameter of the same kind, and the number worth having is its edge: the smallest τ that still keeps the floor, which is also where the run needs the fewest steps to reach it.

Read on a grid of twenty values of τ a decade, the edge at 1% noise is τ = 1.0·10⁻², against an oracle λ of 2.2·10⁻² — a ratio of 0.45. At 0.1% noise it is 1.3·10⁻³ against 3.5·10⁻³, a ratio of 0.36. At 0.01% it is 5.0·10⁻⁴ against 1.1·10⁻³, a ratio of 0.45 again. Across two decades of noise the edge moves with the oracle and stays at a third to a half of it, and the step count at the edge is 5, 10 and 12 where the plain run needs 20, 43 and 74. Half a decade below each edge the best error is 0.234, 0.452 and 0.106 — a loss of 64%, of a factor of four, and of 9%, so how sharply the floor is lost past the edge depends on the noise, and where it is lost does not.

That the edge tracks λ rather than $\lambda^2$ is a matter of units rather than a new law. τ is compared with $\mu_k$, which stands in for a singular value; the shift was compared with $\mu_k^2$, and a shift at $\lambda^2$ is a threshold at λ. Why the edge sits below λ rather than at it has a reading that fits the three ratios. Tikhonov’s filter $\sigma^2/(\sigma^2 + \lambda^2)$ keeps a fifth of a component whose singular value is λ/2 and a tenth of one at λ/3, so at the edge the preconditioner has begun accelerating components the best penalty already damps to between a tenth and a fifth of themselves. That is consistent with the measurement and is not established by it. The approximation’s eigenvalues are not the blur’s singular values, and the correspondence between a cosine above τ and a singular vector above τ is loosest exactly where the two operators disagree.

The low-noise sweep also shows what the edge is not. It is not the τ with the best error: at 0.01% noise every τ from 0.3 to 10⁻³ reaches 0.0975 to 0.0979, and the best of them is chosen by the fourth significant figure. The accuracy is flat and the work is not — the steps fall from fifty-four at the top of that stretch to twelve at its edge — so the edge is the τ worth running at, and it is worth running at because it is the cheapest point that costs nothing.

Three amendments for an approximate preconditioner #

Set beside the exact preconditioner, the two cheap ones change the account of preconditioning a stopped run in three places.

How an approximation is made invertible is part of the preconditioner. A shift and a truncation at the same threshold cluster the same number of eigenvalues near one and differ by a factor of five in the best answer the run can reach. On a well-posed system the two would be interchangeable; on this one only truncation is a preconditioner at all.

A disagreement with the operator is harmless until it is divided. The Fourier approximation is four times further from the blur than the cosine one, and truncated it is as good. The disagreement matters exactly where the construction multiplies it by the reciprocal of something small.

The parameter is still a regularisation parameter. Its edge moves with the noise, at a third to a half of the oracle’s λ over two decades, and above the edge it trades steps for nothing while below it every half decade costs accuracy. That is a parameter that counts steps read one level up: the step count regularised the plain run, the shift regularised the exact preconditioned one, and the truncation level regularises the preconditioner itself.

What the measurement does not cover #

The disagreement between these approximations and the blur is confined to its ends, because the blur is the same kernel in every interior row. A spatially varying blur, or a preconditioner built from an incomplete factorisation, disagrees with the operator in the interior as well, and there the directions the approximation gets wrong need not be its high cosines. Whether truncation still separates the trustworthy directions from the rest when the disagreement is not a boundary effect is not measured.

The problem is one-dimensional, with one signal and white noise. In two dimensions the boundary is a smaller share of the unknowns and a larger number of them, and the rank of the disagreement grows with the side of the image rather than staying at twenty-two. The transforms here are applied as dense orthogonal matrices, which produce the same vectors a fast transform would; the n log n cost is the standard one and is not timed.

And every best step and every edge is chosen against the truth, which a real deconvolution does not have. An edge at a third to a half of the oracle’s λ is a statement about where the right τ lies, not a way of finding it.

Still open: can τ be read off the data, and can the run be stopped without the truth? #

A rule for τ. The edge sits at a fixed fraction of a quantity nobody knows. A rule that places τ from the data alone would have to find, in the approximation’s own coefficients, the index where the data stops carrying signal — the job the Picard plot does for the singular values in where the answer stops being in the data, and does with a systematic overshoot. Whether the same crossing read off cosine coefficients lands above or below the edge measured here is the first measurement a practical code would need.

A stopping rule on the truncated run. A run that reaches its best in five steps instead of twenty has fewer steps near its best to stop at. A stopping rule that follows the run it is given measures what that does to the discrepancy principle, over forty draws of the noise and on both sides of the edge.

A disagreement in the interior. The measurement that would say whether truncation is a general repair or a boundary one is the same sweep on a blur whose width varies across the interval, approximated by the blur of average width.