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.

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.