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%.

Whether the consequence is binary too #

Normality is binary as a property. Its consequence for the weights is a number, and three operators at 0, 0.086 and 0.518 cannot say what shape that number has. Blending the symmetric blur into the row-normalised one, holding the right-hand side fixed so that only the operator moves:

The misfit is proportional to the asymmetry over the first two decades, with a constant near a half and no threshold anywhere. The same holds at other step counts — 0.6 at four steps, 0.55 at eight, 0.87 at twelve — so the constant drifts slowly upward with the work while the proportionality does not break. Past about two per cent the growth turns superlinear, which is where all three of the operators above sit, and is why three points could not have shown the linear part.

So the asymmetry is the quantity that matters, and it supplies the decision procedure this section was written to say does not exist: measure ‖A − Aᵀ‖/‖A‖, multiply by about one, and that is how far the spectral vocabulary is from describing the method. A matrix a tenth of a per cent asymmetric has weights that miss a function of σ by a twentieth of a per cent, which is nothing.

What that leaves standing is the finding itself, and it is worth separating from the mechanism. The row-normalised blur’s 8.6% is a real effect not because a small departure has a large consequence, but because 8.6% is not small — and it arrived from a modelling decision taken by somebody not thinking about eigenvectors, which is the part of the story the proportionality does not touch. A continuous consequence is easier to reason about than a cliff, and it makes the original point sharper rather than softer: the decision that cost 8.6% could have cost 0.1%, and nothing in the construction was watching.

The proportionality has one more consequence worth spelling out, because it decides who has to care. An operator assembled from a discretisation is asymmetric by whatever the discretisation makes it — a boundary treatment, an upwinding, a non-uniform mesh — and those are per-cent effects, so the spectral reading of a Krylov method on them is wrong by per cents and the vocabulary survives. An operator assembled from a model — a shifted kernel, a directional blur, an advective term — is asymmetric by tens of per cents, and there the vocabulary does not survive at all.

That is a division by provenance rather than by field, and it is checkable in one line before any iteration runs. The essay’s three operators were chosen to span the range and, measured, they span only its upper half: every one of them is past the point where the linear regime ends. A fourth operator at a tenth of a per cent would have been the one showing that the question has a small answer as well as a large one.

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.

“Above one” is a claim about a sign, and a sign measured at one step count is a claim about that step count. It is also the kind of claim a reader will assume grows with the iteration, so both the sign and its trend are worth reading off the slider.

The overshoot does not grow with the step count, which is the first thing the slider corrects — and the two frames so far bracket a minimum rather than showing the start of a climb, so the rest of the range is where the climb actually begins.

Every factor exceeds one, at all six step counts — 1.104, 1.015, 1.038, 1.070, 1.120, 1.203. That is the essay’s claim, and the interesting part is that it holds at two steps as well as sixteen: the overshoot is not something the iteration accumulates, it is there from the start.

And it is not monotone. The minimum is at four steps, 1.015, with two steps above it at 1.104. So a reader would not find the smallest overshoot by taking the fewest steps, and the shape of the dependence is a U rather than a ramp. Nothing in the mechanism suggested that, which is why it is worth having measured.

The two curves approach one from opposite sides. Tikhonov’s ceiling rises 0.8234, 0.8997, 0.9506, 0.9670, 0.9866, 0.9952 — always under one, as its formula σ²/(σ² + λ²) requires — while CG’s factor rises past it. At sixteen steps Tikhonov is 0.0048 below one and CG is 0.203 above, so the overshoot is forty-two times the undershoot. They are not two settings of one kind of filter; one of them is a weighting and the other is not.

One row is the instrument rather than the method. The measured-against-predicted agreement runs 8.13·10⁻¹⁴, 6.15·10⁻¹⁴, 5.20·10⁻¹⁴, 6.88·10⁻¹⁴, 1.01·10⁻¹³ and then 1.16·10⁻¹¹ — a hundredfold jump at the last step. Both routes are still agreeing to eleven digits, so nothing is wrong; but the recurrence-derived prediction is beginning to lose the basis’s orthogonality, and a reader extending this figure to thirty steps should expect that column to be the first thing to go.

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.

A filter that is a matrix function #

Every polynomial in A is a filter on its spectrum, and every Krylov method chooses one. Changing the function being filtered towards turns a solver into a matrix-function evaluator with no change to the recurrence.