The basis decides what a filter is

The regularisation field has one vocabulary and it is spectral. A method decides how much of each singular component of the data to believe; the weights it applies are filter factors; the parameter decides how sharply they fall; and every method in the field is one line of arithmetic away from every other, because they all hand the same sum a different f.

A parameter that counts steps extended that vocabulary to an iterative method, and did it properly: the weights of a CGLS iterate are 1 − Π(1 − σ²/θⱼ), a polynomial in σ² whose roots are the Ritz values, and the two routes to them — measured off the iterate, predicted from the recurrence — agree to 10⁻¹³.

This essay is about what that argument rests on, which turns out to be a property of the operator that nobody checks.

The method it is checked against is GMRES, and the reason for choosing that one is not that it is unusual. It is the standard iterative method for a square system that is not symmetric, it is routinely applied to deblurring and to every other discretised inverse problem whose operator happens not to be symmetric, and the sentence “the iteration count acts as a regularisation parameter” is applied to it in exactly the words it is applied to conjugate gradients.

Two spaces, and what each one makes of σ #

CGLS and LSQR generate the Krylov space of AᵀA, starting from Aᵀb. So an iterate is p(AᵀA)Aᵀb; AᵀA’s eigenvectors are exactly A’s right singular vectors; the coordinate along vₖ is σₖ p(σₖ²)(uₖᵀb); and the weight — the coordinate divided by what an unregularised solve would put there — is σₖ² p(σₖ²). A polynomial in σ, and an even one.

GMRES generates the Krylov space of A, starting from b. An iterate is q(A)b, the polynomial acts on A’s eigenvalues in A’s eigenvectors, and if those coincide with the singular vectors the weight is σₖ q(σₖ). Also a polynomial in σ — with odd powers in it.

Two consequences, and the second is the essay.

The two are not the same kind of object even when both are functions of σ. σq(σ) is not bounded by one, is not positive, and has no cut-off; it is not a filter factor in the sense the regularisation field means, and no reordering makes it one.

And GMRES’s is a function of σ only if the eigenvectors are the singular vectors — that is, only if A is normal. For anything else the polynomial acts in a basis that has nothing to do with the singular basis, and there is no reason for the weights to be a function of σ at all.

Turning it into a number, twice #

Every vector has coordinates in the singular basis, so measuring them proves nothing on its own — krylovreg.js says exactly that about its own measurement and builds a second route for the reason. Here the second route is a fit: take the measured weights, fit the best polynomial of the right form by least squares, and report how far it misses relative to the size of the weights.

The right form is where the first version of this measurement went wrong. Fitting even powers only, GMRES’s weights on an exactly symmetric operator miss by 0.247 — which reads as evidence about non-normality and is evidence about parity. Each method has to be fitted against the space its own derivation puts it in.

Four steps, sixty-four unknowns, 1% noise. The left column does not move. The right column moves by twelve orders of magnitude between the first two rows.

Eight and a half per cent, and where it comes from #

The middle row is the finding, and it took a moment to believe.

regular.js builds its blur by normalising every row to sum to one. The reason is stated where it is built and is a good one: it makes the operator a blur rather than a blur-and-scale, so that a constant signal survives it and the ill-posedness is entirely in the smoothing. It is a decision about what the test problem means, taken by somebody not thinking about eigenvectors.

It also destroys the symmetry. The rows near the two ends of the domain have less kernel to sum, so they are divided by smaller numbers, and the matrix comes out 8.6% asymmetric.

That is an amount nobody would mention. It moves GMRES’s weights from lying on a function of σ at the level of rounding to missing one by 4.4% — and at eight steps by 6.3%. The property being lost is binary, and the quantity measuring how much of it is left is not the quantity that matters.

The third operator, and the one it replaced #

The shifted blur moves the kernel’s centre one place off the diagonal. That keeps everything that makes the problem ill-posed — Toeplitz, rows summing to one, singular values decaying exponentially with the smallest consecutive ratio 0.318 against the symmetric problem’s 0.346 — and takes the asymmetry to 0.518 and κ to 5.2·10¹³ against 5.7·10¹².

The first version of that operator narrowed one side of the kernel instead of moving it. It produced a matrix 68% asymmetric and — measured — a condition number of 254. Ten orders of magnitude better conditioned than the problem it was supposed to be a variant of, because a narrower kernel blurs less and the blurring is the ill-posedness.

A test problem that differs in more than the property under test settles nothing, and a comparison run on that operator would have shown GMRES doing well for a reason that had nothing to do with the argument. It is recorded in the library beside the function that replaced it.

And GMRES is not the worse method #

The half of this that stops it being a verdict.

GMRES reaches an answer within 12% of LSQR’s for a tenth of the products — one product a step against two, and a best step of three or four against twenty. Where the matrix is only available through its action on a vector, that is the whole of what an iterative method is for.

What it costs is that the good iterate arrives and departs. By the end of a thirty-two-step run GMRES is at 1.8·10² on the shifted problem — a thousand times its own best — where LSQR is at 0.31. A method that is best at step four and useless at step thirty needs a stopping rule far more urgently than one that is best at step twenty, and the rules available are the ones this field already scored against an oracle, none of which reached it.

What the weights look like when they are not a function of σ #

The misfit is a summary. The weights themselves say the same thing more plainly.

On the shifted operator eight of the sixty-three usable components come back with a negative weight — the method has subtracted part of a component of the data from the answer — and one comes back at 1.90, nearly doubled. Neither is a defect: the iterate minimises a residual over a space and is under no obligation to be a reweighting of anything.

But “keeps the components the data supports and discards the rest” is not a description of a computation that multiplies one of them by −0.13. LSQR’s weights, on the same operator and the same data, stay in [0, 1.054] with none negative.

The overshoot above one, for the record, is not the interesting part. The CG filter overshoots too — measured at 1.20 on this site — because a polynomial that is one at the converged Ritz values and zero far below them cannot be monotone in between. That is a property the two methods share. The sign changes are not.

The same non-normality, in a place this site already looked #

The iterative field has an essay about GMRES whose finding is that the spectrum does not predict the convergence: for a non-normal matrix, any convergence curve is compatible with any spectrum, so the eigenvalues — the obvious thing to look at — say nothing about how fast the residual falls.

That is this essay’s finding from the other side, and the quantity underneath both is the same one. The departure from normality, measured as ‖AᵀA − AAᵀ‖/‖A‖², is exactly zero on the symmetric blur, 1.9·10⁻² on the row-normalised one and 2.7·10⁻² on the shifted one — and the property it destroys is that the eigenvectors of A are an orthonormal basis in which everything can be read off.

When they are, GMRES’s polynomial acts diagonally in a basis the singular value decomposition also uses, its convergence is governed by the eigenvalues, and its weights are a function of σ. When they are not, both statements fail together, and they fail for one reason rather than two.

That is worth knowing in the direction it is usually needed. Nobody sets out to build a non-normal operator; the middle row of the table above got its asymmetry from a normalisation intended to make the test problem cleaner. So the question a reader should carry away is not is my operator pathological but is it exactly normal, and the answer for a discretised integral operator with boundary conditions on it is almost always no.

Why this matters more than a taxonomy #

The practical consequence is not that GMRES should be avoided. It is that a sentence true of one Krylov method is applied to all of them, and the sentence is load-bearing.

“The iteration count is a regularisation parameter” means something specific: each extra step admits another band of singular components, so stopping is choosing a cut-off, so the parameter is comparable with a truncation index and with a Tikhonov λ. That is exactly what four knobs and one floor measured, and it is why four methods sharing no arithmetic land within 3% of each other.

For GMRES on a non-normal operator none of it follows. Each step admits another polynomial in A, and what that does to the singular components is not describable by a cut-off. The step count is still a knob — the error still turns, and turning it still trades noise against signal — but it is not the same knob, and the identification that makes the combination field’s reading work is unavailable.

The repair, and what it costs #

There is an obvious fix, and this site has an essay about why nobody takes it.

Apply the Arnoldi method to AᵀA rather than to A. That matrix is symmetric whatever A is, so its eigenvectors are the right singular vectors, the polynomial acts in the singular basis, and the weights become a function of σ again — measured at 6.5·10⁻¹³ against an even polynomial after eight steps on the shifted operator, which is the level of rounding.

The price is on the label. κ(AᵀA) on that problem is 1.0·10¹⁷ against κ(A) = 5.2·10¹³ — the square, and past 1/u, so the product is numerically singular. The road that squares the problem is about exactly this trade, from the least-squares field, and its conclusion applies unchanged: the repair that restores the description is the one that destroys the conditioning.

Which is what LSQR already is. Its Krylov space is the space of AᵀA, generated without ever forming it, and that is the whole reason it is written as a bidiagonalisation rather than as three lines of matrix products. The choice on offer is not between a describable method and a fast one; it is between building the space of AᵀA carefully, building it carelessly, or building a different space and giving up the description.

What the fit cannot see #

The test has a limitation, and refusing to state it would make the essay weaker than it is.

A misfit at the level of rounding is evidence for a function of σ. A large misfit is evidence against a low-degree polynomial, which is not the same as evidence against any spectral description — a truncation’s weights are a step function, an entirely legitimate filter, and no degree-12 polynomial fits a step either. That case is refused explicitly in the library: feeding the fit a truncation’s own weights and claiming they lie on a polynomial has to fail, and it does.

The parity distinction is refused in the same way, in both directions. Claiming GMRES’s weights lie on an even polynomial fails on every operator including the symmetric one; claiming LSQR’s lie on a general polynomial of the same degree fails too, at 2.1·10⁻². Two spaces, neither containing the other, and a gate that would notice if some later change collapsed them into one.

What the drag does #

The slider is the asymmetry of the operator, and it has three positions because there are three operators — a blur scaled by one number, the same blur row-normalised, and one with its kernel moved. Nothing else about the problem changes between them: the same signal, the same seed, the same noise level, the same exponential spectrum.

The bidiagonal weights do not notice any of it, at any position, because their parity is a property of AᵀA and A’s symmetry does not enter. The Arnoldi weights go from the level of rounding to 4.4% at the first step of the slider and to 53% at the second.

There is one thing the slider does that is worth watching for and is not the argument: at higher step counts the bidiagonal misfit itself drifts — 10⁻¹⁴ at two steps, 10⁻⁹ at twelve, 10⁻⁵ at sixteen. That is the same expiry date the filter identity has, the description degrading with the orthogonality of the recurrence underneath it, and the assertion is written in two halves either side of twelve steps rather than quietly restricted to the range where the simple version holds.