One sequence and two recurrences

This site’s oldest argument is that two algorithms can be one derivation and two computations. Two Gram–Schmidts reorders the same projections and moves ‖QᵀQ − I‖ by fourteen orders of magnitude. The road that squares the problem forms AᵀA, which is algebraically harmless and numerically the end of the matter.

Here is the same shape in a pair nobody teaches as a pair.

Both CGLS and LSQR minimise ‖b − Ax‖ over the Krylov space generated by AᵀA from Aᵀb, at every step. The minimiser over a space is unique. So their iterates are equal, not similar, and nothing in the algebra can prefer one.

What each recurrence keeps #

CGLS is conjugate gradients applied to the normal equations without forming them. It carries a residual r, a normal residual s = Aᵀr, a search direction p, and the scalar γ = ‖s‖². Each step is one product with A, one with Aᵀ, and three vector updates.

LSQR carries the Golub–Kahan bidiagonalisation instead. The recurrence produces a lower bidiagonal matrix Bₖ one column at a time, and a plane rotation is applied to it as each column appears, so that the growing least-squares problem stays solved without ever being re-solved. What is stored is three vectors and two scalars, and the cost per step is the same: one product with A, one with Aᵀ.

Two implementations of one bidiagonalisation stand behind this collection — this one and the hybrid method’s, which keeps its bases because its subject is the projected matrix. The α and β they produce are asserted to agree, and do, to zero difference in every entry.

The measurement #

A 60×12 matrix built to a prescribed condition number, a consistent right-hand side, and the question: how many steps to reach a relative error of 10⁻⁶?

Four seeds, no overlap between the two groups at any condition number, and a separation that grows monotonically with κ. That is not a matrix being unlucky.

Two things about this table are worth being explicit about, because both are limitations rather than findings.

The Krylov space is exhausted at twelve. A twelve-column matrix has a Krylov space of dimension at most twelve, so in exact arithmetic both methods terminate at step twelve with the exact answer. Every step in the table past the twelfth is a step neither method should be taking, and both take them, which is the CG is not a direct method finding this site has asserted since its foundation. This whole comparison lives in the region where the algebra has already finished.

And the target is chosen against a floor. The forward error cannot go below about κu, so at κ = 10¹⁰ the floor is 1.1·10⁻⁶ and the target of 10⁻⁶ is close to it. Asking for 10⁻⁸ would measure the floor rather than the methods — the assertion refuses to draw it — and that boundary belongs to the error field’s identity rather than to either recurrence.

The penalty is a function of two things and the essay’s headline names one of them. Across targets of 10⁻³, 10⁻⁴, 10⁻⁵ and 10⁻⁶ the worst ratio reads 1.57, 1.63, 1.84 and 1.90, and at κ = 10⁶ the two methods take 39/44, 40/45, 41/45 and 42/47 — five steps apart at every target, or about twelve per cent.

So the second recurrence costs nothing worth measuring on a well-conditioned problem however accurately it is solved, and its cost on a badly conditioned one grows with the accuracy demanded: 58, 66, 89 and 99 extra steps at κ = 10¹⁰ as the target tightens by three decades. The conditioning decides whether there is a penalty and the target decides how large it is.

Iterations, not digits #

A table of step counts says one recurrence gets there sooner. It does not say whether it gets further, and those are different things to buy — so here is the best relative error either method ever reaches, over four seeds, median:

at κ = 10⁴, 7.5·10⁻¹⁴ for the bidiagonalisation against 2.9·10⁻¹⁴ for the normal equations. At 10⁶, 3.3·10⁻¹² against 5.6·10⁻¹². At 10⁸, 1.1·10⁻⁹ against 4.3·10⁻¹⁰. At 10¹⁰, 1.5·10⁻⁷ against 8.8·10⁻⁸. At 10¹², 8.7·10⁻⁶ against 2.1·10⁻⁵.

The ratio runs 0.38, 1.70, 0.40, 0.60, 2.41 — no trend, and no consistent winner: the normal-equations route is ahead at three of the five condition numbers. Both floors sit at about a tenth of κ·u, which is the problem’s floor rather than either recurrence’s, and neither is nearer it than the other.

So the whole of the advantage is a rate. It buys iterations and not digits, and a run that has converged has nothing left to gain from the choice. That is worth stating beside the step counts, because a 1.9× step advantage at κ = 10¹⁰ invites the inference that the same factor buys accuracy somewhere, and it does not: at that conditioning both methods stop at 10⁻⁷, one of them a hundred steps sooner.

Which is also why the step comparison stops where it does #

One reading to rule out before the section closes, because the floors invite it. The two floors are not equal because the two methods are the same method — they are not the same method in floating point, which is the whole page. They are equal because the floor belongs to the problem: a forward error cannot fall below about κ·u whatever computes it, and both recurrences are backward stable enough to reach that. Two different computations arriving at one floor is a statement about the floor.

What would distinguish them at the floor is a method that is not backward stable, and this collection has one to hand: forming the normal equations squares the condition number, so its floor is κ²·u rather than κ·u and it stops six orders earlier at κ = 10⁶. That is the comparison in which a floor does separate two routes, and it is the reason CGLS is written the way it is — it is conjugate gradients on the normal equations without forming them, and the difference between that and forming them is the difference between sharing this floor and having a much worse one.

The section above notes that the target of 10⁻⁶ sits close to the floor at κ = 10¹⁰ and that asking for 10⁻⁸ would measure the floor rather than the methods. The floors measure how close that boundary was.

At κ = 10¹² the floor is 8.7·10⁻⁶ — above the 10⁻⁶ target the step counts are measured against. Neither method reaches it, at any number of steps, so the comparison has nothing to report and the table cannot be extended by one more row.

That is a hard stop rather than an affordability limit, and it is worth distinguishing from the other kind. Most sweeps on this site end where the arithmetic gets expensive — a dense reference decomposition, a Gram matrix, a headless browser. This one ends where the question stops existing: past κ = 10¹⁰ there is no accuracy at which both methods can be compared, because the accuracy they would be compared at is one neither can deliver.

And it locates what the comparison is really about. The two recurrences differ only in the region between where the algebra has finished — twelve steps, since the Krylov space is exhausted — and where the arithmetic has finished, at the floor. Both of those boundaries belong to the problem. The methods differ in how fast they cross the gap between them, and the gap closes as κ rises: at 10⁴ it is a hundred steps wide, at 10¹² it has closed entirely.

What a step of each one is #

The two recurrences are the same length and store the same amount, which is why the comparison is about arithmetic rather than about cost per step.

CGLS. Update the iterate along the search direction, update the residual by the same multiple of Ap, form the new normal residual Aᵀr, and take the ratio of two consecutive values of ‖Aᵀr‖² as the coefficient that makes the new direction conjugate to the old one. Two products with the matrix, three vector updates, and one scalar ratio that carries the whole conjugacy.

LSQR. Extend the bidiagonalisation by one column — a product with A, orthogonalise against the previous left vector, normalise; a product with Aᵀ, orthogonalise against the previous right vector, normalise. That produces one new α and one new β. Then apply a plane rotation chosen so that the new β is annihilated, which keeps the growing bidiagonal least-squares problem in solved form, and step the iterate along the rotated direction.

Neither keeps a basis. Both are short recurrences of fixed cost. The difference is which quantities the arithmetic is done on: LSQR’s α and β are norms of vectors the size of the problem’s own scale, and CGLS’s γ is a squared norm and its coefficients are ratios of squared norms.

One suspect can be eliminated straight away, and it is the usual one. The residual CGLS propagates — updated by a recurrence rather than recomputed from x — could drift away from the true residual, and on this problem it does not: at κ = 10¹⁰ the two agree to a part in forty at every step drawn, and the only place they part company is at the stagnation floor, where the propagated value reads 7.1·10⁻¹⁶ against a recomputed 1.1·10⁻¹⁵. That is a factor of 0.65 at the bottom of a run that has already finished, and it is not what costs a hundred steps.

Where the difference comes from #

Both methods generate a basis for the same space and neither reorthogonalises. So the obvious place to look is the quantity this site measures more often than any other.

Measured over the first four vectors of each, and that restriction is the honest part: past twelve vectors both bases are made of rounding error and both orthogonality errors are of order one, so there is nothing to compare. Over the range where the comparison means anything the bidiagonal basis is the better one at every condition number, by a factor of about five.

A factor of five in a basis and a factor of two in a step count is not a derivation, and it is not offered as one. What it is, is the same quantity separating the two methods that separates every other pair on this site: orthogonal is a number, and the method whose number is smaller is the one that finishes first.

Why the normal-equation recurrence is the one that suffers #

The name is the clue, and it is worth stating carefully because the obvious version of the argument is wrong.

CGLS does not form AᵀA. That is the whole reason it is written the way it is, and this site has an essay about what forming it would cost: κ(AᵀA) = κ(A)², and below ε = √u the product is exactly singular. Nothing in CGLS is computed at κ².

What CGLS does do is monitor a squared quantity. Its recurrence coefficients are built from γ = ‖Aᵀr‖², and its conjugacy relations are relations in the AᵀA inner product. The vectors it generates are therefore conjugate with respect to a matrix whose condition number is κ², even though that matrix never exists — and conjugacy with respect to an operator is exactly the property that rounding error destroys in proportion to that operator’s conditioning.

The bidiagonalisation avoids that by keeping the two sides apart: its α come from Aᵀu and its β from Av, each of them a κ-level quantity, and the least-squares problem it solves at each step is bidiagonal rather than normal.

The precision does what the conditioning does #

The separation is a property of the matrix and the arithmetic, and this site has a slider for the second of those. One matrix, one target, three precisions:

Twenty-one significand bits fewer, on the same κ = 10⁶ problem, and the separation widens from 11% to 46% — reproducing what raising κ by four decades did at double precision. At 24 bits the target of 10⁻⁴ is below the forward-error floor and neither method reaches it, which is the boundary again.

That immediately suggests a single governing parameter, and this site has one to hand: κu, the product that decides the mixed-precision field’s threshold and the interval field’s boundary. It is the natural guess and it is wrong.

The second case has the larger κu and much the smaller separation. A quantity that orders two measurements the wrong way round is not the quantity that governs them, and that is asserted here rather than being left as a sentence — the assertion states the counterexample and would fail if some later change made κu fit after all.

What it does not change #

Three things, and they matter because a factor of two in an iteration count is easy to over-read.

The answers agree. At the step each reaches its own best, the two iterates differ by less than either differs from the truth. Neither method is returning a wrong answer; one is taking longer to arrive.

The floor is the same. Both bottom out at about κu, which is where the problem puts them. No recurrence can do better and neither does.

And at moderate conditioning there is nothing in it. At κ = 10² the two agree to 3.3·10⁻¹³ at step six and cost the same number of steps. A reader whose problems are well conditioned has been choosing between these two methods on documentation and taste for thirty years, and has been right to.

Where the difference does not appear at all #

The two methods that this essay separates are the two methods the regularisation field uses, and on that field’s own problem the separation is invisible.

On the 64-point deconvolution at 1% noise both reach their best error at step 20, and both reach 0.142587. Their iterates agree to 5.3·10⁻¹⁶ at step five, 9.2·10⁻¹⁶ at step ten, and only 7.9·10⁻⁶ at step twenty — by which point the run is over, because the error has turned and everything after it is worse.

The reason is worth stating because it is the practical answer to which should be used. Semi- convergence stops an ill-posed run after a couple of dozen steps, and a couple of dozen steps is before either recurrence has degraded. The comparison in this essay needs a hundred steps to become visible, and no ill-posed problem gets a hundred useful steps.

So the pair separates on consistent, ill-conditioned problems run to convergence, and not on ill-posed ones stopped early. Those are different jobs, and the difference between the two recurrences matters for exactly one of them.

And what it says about a rule this site keeps finding #

The pattern is now familiar enough to be worth naming. Classical against modified Gram–Schmidt. The normal equations against a QR factorisation. Elimination with and without a row interchange. A reduction tree against a Householder sweep. And now two recurrences for one sequence of iterates.

In every one of those pairs the algebra is silent and the arithmetic is not, and in every one of them the difference is invisible until the conditioning is bad enough. That is what makes the class dangerous rather than merely interesting: a method chosen on well-conditioned test problems, where every member of the pair agrees to twelve digits, carries a factor that only appears on the problems that were hard enough to need a numerical method in the first place.

Which one to write, then #

The answer this measurement supports is narrower than a recommendation, and it is worth stating in the narrow form.

If the problem is well conditioned, the two methods are the same method and the choice is documentation. If it is ill-posed and stopped early, the two are the same method for the length of the run that matters. If it is consistent, ill-conditioned, and run to convergence — a least-squares fit with a nearly rank-deficient design, solved iteratively because the matrix is too large to factorise — the bidiagonal recurrence costs half as many products with A, at four seeds and at every condition number above 10⁶.

That third case is the one the method was published for, and it is also the case where the alternative was already available and already looks fine on any test problem short of it. Which is the shape of every entry in this class: the pair agrees exactly where a comparison is easy to run, and separates exactly where running it is expensive.

What the drag does #

The slider is the accuracy asked for. A tighter target costs both methods more steps and costs the normal-equation recurrence proportionally more of them: at 10⁻³ the ratio at κ = 10¹⁰ is smaller than at 10⁻⁶, because the early steps of the two runs are nearly identical and the divergence accumulates.

The slider stops at 10⁻⁶ and the reason is the floor. At κ = 10¹⁰ the forward error cannot reliably go below 1.1·10⁻⁶ whatever either method does, so a target of 10⁻⁷ is a target neither reaches and the figure would be drawing a boundary of the problem rather than a difference between two methods. The generator refuses it, which is the assertion doing what the site’s assertions are for.

A third thing the same recurrence computes #

One Krylov space, several algorithms. Change what is asked of the small projected matrix and the same recurrence returns a matrix function rather than a solve — and the convergence is superlinear rather than governed by a rate.