Numerical Linear Algebra
where the algebra stops being the arithmetic
Iterating, instead of factorising

The number that is re-derived

GMRES prints a residual it never computes from its answer either. On the matrix that sends a conjugate gradient recurrence 7.3·10¹⁰ wrong, and on two others chosen to be worse, its number is never more than a factor of 2.86 out — while the basis it is computed from has lost orthogonality entirely. The disease is not iterative methods, and it is not floating point.
28 min read 25 figures Identical algebra, different arithmeticMeasured, not assertedWhose fault is itTwo routes to a number

Worth reading first: The residual the method reports · The spectrum that predicts nothing · Orthogonal is a number.

The essay on what a method reports ends with a nine-and-a-half-order gap between a printed residual and a real one, and the natural conclusion to draw from it is the wrong one.

The wrong conclusion is that iterative methods print residuals that cannot be trusted. The right one is narrower and far more useful, and getting to it requires a method in the same family that does not have the problem.

GMRES is that method. It also prints a residual it never computes from its answer. Its number is honest.

What GMRES reports and what its answer's residual is, on a matrix with one eigenvalue at 1e-14GMRES also prints a residual it never computes from its answer — but it does not carry one. The number falls out of the Givens rotations at each step as the residual of the projected least-squares problem, re-derived from the whole basis rather than updated from the last value. Over 40 steps the two never separate by more than a factor of 1.80, and they do that while ‖VᵀV − I‖ for the Arnoldi basis is 1.41 — which is to say the basis has stopped being orthogonal and the residual computed from it is still honest.051015202530354010⁻³10⁻¹10¹GMRES steprelative residualthe two curves are the same curvea number re-derived, not carriedworst reported/actual factor1.8at step40‖VᵀV − I‖ of the basis1.4reported at the last step0.041actual at the last step0.074the same family of methodsand only one of them lies
Fig. 1 GMRES on the matrix that sends a conjugate gradient recurrence 7.3·10¹⁰ wrong. The reported residual and the answer’s are the same curve.

Where GMRES’s number comes from

GMRES builds an Arnoldi basis and solves a small least-squares problem in it. The residual it reports is not updated by a recurrence at all: at each step, the Givens rotations that triangularise the Hessenberg matrix carry the right-hand side along with them, and the last entry of the rotated right-hand side is the residual norm of the projected problem. Exactly, as a consequence of the rotations being orthogonal.

So the quantity is re-derived at every step from the whole of what the method currently holds, rather than updated from its own previous value. That is the difference from a conjugate gradient residual, and it is the only difference that matters here.

GMRES on the Laplacian and on the cyclic shift, both 40×40A semi-logarithmic plot of relative residual against iteration. One curve falls steadily; the other is flat at one for every step until the last, where it drops to zero.071421283510⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1iteration‖r‖ / ‖b‖Laplaciancyclic shiftno progress at allevery eigenvalue of the shift is on the unit circleand it predicts nothing
Fig. 2 GMRES’s residual is monotone by construction, from the essay that established it — and the spectrum does not predict it, which is a separate fact about the same curve.
The residual basis of conjugate gradients, at κ = 105A semi-logarithmic plot against iteration count showing the loss of orthogonality among the residual vectors rising while the relative residual falls.05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here
Fig. 3 The basis both methods build, and what happens to its orthogonality, from the field that measured it.

Three matrices chosen to break it

The comparison is only worth having if the matrices are hostile, so all three are.

The diagonal one with an eigenvalue at 10⁻¹⁴. The same matrix the recurrence essay is built on, where the conjugate gradient recurrence goes 7.3·10¹⁰ wrong. Worst reported-to-actual factor: 1.80.

A bidiagonal matrix with a superdiagonal thirty times its diagonal. Strongly non-normal, the kind of operator whose eigenvalues predict nothing about its behaviour. Worst factor: 2.86.

A single Jordan block. The most degenerate matrix in the subject: one eigenvalue, one eigenvector, and a Krylov space that cannot resolve it. Worst factor: 2.00.

Not one of the three produces a discrepancy worth a sentence. The largest is a factor of under three, at one step out of forty, and the two curves are visually one curve on all three matrices.

What GMRES reports and what its answer's residual is, on a matrix with non-normal, superdiagonal 30GMRES also prints a residual it never computes from its answer — but it does not carry one. The number falls out of the Givens rotations at each step as the residual of the projected least-squares problem, re-derived from the whole basis rather than updated from the last value. Over 40 steps the two never separate by more than a factor of 2.86, and they do that while ‖VᵀV − I‖ for the Arnoldi basis is 4.48·10⁻¹⁰ — which is to say the basis has stopped being orthogonal and the residual computed from it is still honest.051015202530354010⁻⁶10⁻⁴10⁻²1GMRES steprelative residualthe two curves are the same curvea number re-derived, not carriedworst reported/actual factor2.9at step38‖VᵀV − I‖ of the basis4.5·10⁻¹⁰reported at the last step7.9·10⁻⁵actual at the last step1.6·10⁻⁴the same family of methodsand only one of them lies
Fig. 4 The non-normal case, where the worst factor over forty steps is 2.86 and both curves fall together.
What GMRES reports and what its answer's residual is, on a matrix with a single Jordan blockGMRES also prints a residual it never computes from its answer — but it does not carry one. The number falls out of the Givens rotations at each step as the residual of the projected least-squares problem, re-derived from the whole basis rather than updated from the last value. Over 40 steps the two never separate by more than a factor of 2.00, and they do that while ‖VᵀV − I‖ for the Arnoldi basis is 1.41 — which is to say the basis has stopped being orthogonal and the residual computed from it is still honest.051015202530354010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹GMRES steprelative residualthe two curves are the same curvea number re-derived, not carriedworst reported/actual factor2at step40‖VᵀV − I‖ of the basis1.4reported at the last step1.5·10⁻¹⁵actual at the last step3·10⁻¹⁵the same family of methodsand only one of them lies
Fig. 5 And a single Jordan block, where GMRES reaches 1.5·10⁻¹⁵ and the vector it holds is at 3.0·10⁻¹⁵.

While the basis is gone

Here is the part that makes this a measurement rather than a reassurance.

The Arnoldi basis is built by modified Gram–Schmidt with no reorthogonalisation, which is what GMRES does. On two of the three matrices, ‖VᵀV − I‖ after forty steps is 1.41.

That is not a loss of orthogonality, it is the absence of it: a value of order one means the basis vectors are no longer even approximately mutually orthogonal, and the projection GMRES is computing is a projection onto something other than the space it thinks it is projecting onto.

And the residual it derives from that projection is still within a factor of two of the truth.

The reason is a theorem rather than luck — modified-Gram–Schmidt GMRES is backward stable, and the loss of orthogonality is exactly compensated by what it does to the projected problem — but the measurement is what this collection is for, and the measurement says: the basis is destroyed and the number is fine.

‖QᵀQ − I‖ at every step of a Lanczos run, n = 40The departure from orthogonality of the Lanczos basis, against the step, on a logarithmic vertical axis. It sits at the level of rounding for 13 steps and then climbs by a factor of about seventeen a step for ten steps running, twelve orders of magnitude, before saturating. The dashed line marks the step at which the first Ritz value converged, which is step 11.051015202530354010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹step‖QᵀQ − I‖first Ritz value converges√uno reorthogonalisationfullthe two events are one eventstep orthogonality crosses √u14step the first Ritz value converges11consecutive steps at ≥5× growth10a drift would grow like a square rootthis grows like a geometric series and then stops
Fig. 6 Orthogonality collapsing rather than decaying, from the essay that established the shape. The basis in this essay is on the far side of that collapse.
Reversing the columns of a 6×6 matrix: what each orthogonalisation returnsBoth grids print the entries of Q₁P − Q₂, where Q₁ orthogonalises A and Q₂ orthogonalises A with its columns reversed. If the answer permuted with the columns, both grids would be zero. The polar factor's is: ‖U₁P − U₂‖_F = 1.53·10⁻¹⁵. QR's is not: ‖Q₁P − Q₂‖_F = 2.353 on matrices whose own Frobenius norm is 2.449, so the two answers are essentially unrelated. Gram–Schmidt and Householder both build Q one column at a time, and the first column is treated differently from the last.polar: U₁P − U₂QR: Q₁P − Q₂0000000000000000000000000000000000000.56-0.15-0.21-0.330.06-0.64-0.430.310.310.120.20.73-0.130.64-0.031-0.0028-0.420.075-0.18-0.43-0.130.21-0.29-0.420.640.630.670.62-0.52-0.420.26-0.520.022-0.0870.2-0.24the same columns, reordered‖U₁P − U₂‖_F1.5·10⁻¹⁵‖Q₁P − Q₂‖_F2.4‖Q‖_F, for scale2.4κ of the matrix10the column space did not moveand one of the two answers did
Fig. 7 And the two orders of the same algebra, from the field that made orthogonality a number. GMRES uses the better of the two and still ends at ‖VᵀV − I‖ of order one.
What each reorthogonalisation costs and what it leaves, n = 40Two quantities for three variants on one logarithmic axis: the number of projections performed, and the departure from orthogonality left behind. Doing nothing costs nothing and leaves a basis that is not one. Full reorthogonalisation costs 820 projections and leaves 1.9·10⁻¹⁵. Selective costs 137 — 6.0 times fewer — and leaves 10·10⁻⁸, which is √u and is enough.010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹count, and ‖QᵀQ − I‖√u — semi-orthogonalitynoneselectivefullprojections‖QᵀQ − I‖two counters, three methodsghosts, none4ghosts, selective0selective is cheaper by6the cheap repair leaves a basis that is not orthogonaland that is the claim, not a shortfall
Fig. 8 What buying the orthogonality back costs, from the essay that priced it. GMRES does not buy it, and this essay is why it does not have to for the purpose of its stopping test.

What the pair actually says

Two methods. The same family, the same Krylov space, the same matrix, the same arithmetic. One of them reports a residual nine and a half orders out and the other reports one within a factor of two.

The difference is not accuracy per step. It is not the conditioning, which is identical. It is not the orthogonality, which is worse in the method that gets it right. It is:

one method carries a vector forward and the other re-derives a number.

A carried quantity accumulates every rounding that has ever entered it, and nothing downstream will ever compare it against anything. A re-derived quantity has no history: whatever was wrong with the last step’s version of it is not in this step’s, because this step’s was computed from scratch.

That is the sentence this field is organised around, and it is worth having it in the form where the two sides are two implementations of the same idea rather than two different subjects.

The residual conjugate gradients reports and the residual of the vector it holds, along one runA 50×50 diagonal matrix with one eigenvalue at 10⁻¹⁴ and the rest between 1 and 50, so every product with it is exact to a rounding and nothing below can be blamed on the operator. The lower curve is r ← r − αAp, updated by the recurrence, which is what a stopping test reads. The upper curve is ‖b − Ax‖/‖b‖, recomputed from the iterate. They start as the same vector and end 7.3·10¹⁰ apart, with the reported one at 6.92·10⁻²¹ — below the unit roundoff of 1.11·10⁻¹⁶, which is the shortest proof available that it is not the residual of anything.01428425670849811210⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best6.9·10⁻²¹the answer's, at its best5.1·10⁻¹⁰unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10¹³iterations drawn110the recurrence remembers every roundingand the stopping test is written in it
Fig. 9 The carried version, from the essay that measured it. Same matrix, same space, same arithmetic.
What residual replacement costs and what it buys, at three periodsEach row is the same conjugate gradient run with the recomputed residual assigned back into the recurrence every k steps. Without it, the reported residual reaches 6.92·10⁻²¹ and the answer's stalls at 5.07·10⁻¹⁰. Replacing every 5 steps costs 30 extra matrix–vector products on top of 150 — 20 per cent — and brings the answer's residual to 8.31·10⁻¹⁷, with the two residuals then agreeing to a factor of 2.8.never replace5.07·10⁻¹⁰replace every 251.5·10⁻¹⁶replace every 101.34·10⁻¹⁶replace every 58.31·10⁻¹⁷reported 6.92·10⁻²¹ · 0 extra productsreported 5.3·10⁻¹⁷ · 6 extra productsreported 1.08·10⁻¹⁶ · 15 extra productsreported 2.95·10⁻¹⁷ · 30 extra productsthe residual of the answer the run returnsone line, at three pricesnever: the answer's residual5.1·10⁻¹⁰never: what it reported6.9·10⁻²¹every 5: the answer's8.3·10⁻¹⁷every 5: what it reported3·10⁻¹⁷extra products for that30ask the matrix againand the recurrence forgets what it did
Fig. 10 And the repair that turns the first into the second, which is exactly “recompute rather than carry” applied every few steps.

The price GMRES pays for it

An honest comparison has to say what the re-derivation costs, because it is not free and the reason conjugate gradients uses a recurrence is a good one.

GMRES keeps the whole basis. Storage grows with the step number and so does the work per step: at step k it orthogonalises against k vectors. That is why it is restarted in practice, and restarting brings its own essay’s worth of problems.

Conjugate gradients keeps three vectors and does a fixed amount of work per step, for as many steps as anybody wants. On a symmetric matrix that is an enormous prize, and it is bought exactly by the short recurrences — including the one for the residual.

So the two methods are not a good one and a bad one. They are two points on a trade this collection has drawn before: a short recurrence is cheap and forgets nothing, and a long one is expensive and starts fresh. The residual gap is one more thing on the bill for the short recurrence, and it is usually the smallest item on it.

Steps to a relative error of 10⁻⁶, on one 60×12 problem at five condition numbersTwo curves against the condition number on a logarithmic horizontal axis, with four seeds drawn at each. The two recurrences compute the same iterates in exact arithmetic. At κ = 10² they cost 16 and 16 steps; at κ = 10¹⁰ they cost 110 and 209.10²10⁴10⁶10⁸10¹⁰0285684112140168196224condition numberstepsnormal equationsbidiagonalisationone sequence, two costssteps at κ = 10², both16at κ = 10⁶, ratio1.1at κ = 10¹⁰, ratio1.9the same iterates in the algebraand twice the work at κ = 10¹⁰
Fig. 11 The trade in its general form, from the essay that established it: what a short recurrence buys and what it costs.
8 restarts keeping 4 of 8, on a 40×40 matrixThe residual bound of the worst wanted eigenvalue and its true error, against the number of products with A. The bound falls from 1.21 to 2.43·10⁻¹⁴ across 8 cycles and 92 products, and the true error reaches 1.42·10⁻¹⁴. The basis is 8 vectors at every cycle and never grows.01224364860728410⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹products with Asizeresidual boundtrue errorbounded memorybasis vectors kept8products with A92worst error in the k wanted1.4·10⁻¹⁴the bound is free and the error is notand the basis never grows
Fig. 12 And what restarting does, from the field that measured it — the usual way of paying for a long recurrence, and not a free one.
The reported bound and the residual it bounds, over 6 cyclesThree quantities against the cycle count on a logarithmic vertical axis. The residual bound the method reports falls without limit, reaching 1.41·10⁻²². The residual it claims to bound stops at 5.68·10⁻⁵ and does not move. Recomputing the arrowhead's border entries, at one extra product with A a cycle, takes the residual to 3.8·10⁻¹⁴.12345610⁻⁵⁵10⁻⁴⁹10⁻⁴³10⁻³⁷10⁻³¹10⁻²⁵10⁻¹⁹10⁻¹³10⁻⁷10⁻¹cyclesizethe residualrepairedthe reported boundwhat the stopping rule readsreported at the last cycle1.4·10⁻²²the residual there5.7·10⁻⁵with the border recomputed3.8·10⁻¹⁴products, cheap and repaired5a bound with nothing under itand one product a cycle to fix it
Fig. 13 Along with the bound that stops being a bound once vectors are carried between cycles, which is the same distinction one level up.

What “re-derived” means, precisely

The word is doing a lot of work in this essay, so it is worth pinning down what distinguishes the two cases, because “computed from scratch” is not quite it.

GMRES’s residual at step k is not computed from scratch either. The Givens rotations are applied incrementally — the rotations from earlier steps are stored and reapplied — so there is a recurrence in the implementation. What is different is what the recurrence is on.

The conjugate gradient recurrence is on the residual vector itself: the quantity being reported is the object being updated, so an error in it is an error in the report, permanently.

The GMRES recurrence is on the factorisation of the projected problem. The reported number is the last entry of a rotated right-hand side, and it is a consequence of a triangularisation whose own errors are bounded by the orthogonality of the rotations — which is exact to a rounding, whatever has happened to the Arnoldi basis. An error in the basis changes which problem is being projected onto, and the reported residual is then the correct residual of that slightly different problem — which is precisely a backward-error statement, and precisely why the number stays honest.

So the distinction is not “recurrence or not”. It is whether the reported quantity is the state of a recurrence, or a function of a state whose own errors are backward errors. The first accumulates into the report. The second accumulates into the problem, where this collection has a whole vocabulary for it.

Backward error with and without the structure, for Levinson and for elimination, n = 12Both solvers are backward stable on this family by the usual measure: the smallest perturbation of any kind that explains either answer stays below 3.1·10⁻¹⁷ at every ρ, and the two curves are indistinguishable. Insisting the perturbation be a symmetric Toeplitz matrix — the kind of object the problem was posed with — lifts both by orders, and at ρ = 0.999 Levinson's rises by a factor of 2.48·10⁶. Neither number is alarming in absolute terms; what is worth carrying is that the reassuring one is the one that is reported, and the two are not measuring the same thing.10⁻³10⁻²10⁻¹110⁻¹⁹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹1 − ρbackward errorLevinson, Toeplitzelimination, ToeplitzLevinson, any kindelimination, any kindat ρ = 0.999Levinson, any perturbation3.1·10⁻¹⁷Levinson, Toeplitz only7.6·10⁻¹¹the ratio between them2.5·10⁶diagonal defect of the first0.97the number that is reportedand the number that is asked about
Fig. 14 The vocabulary in question, from the essay that established it. An error that lands in the problem is a backward error and is describable; an error that lands in the description is not.
One Householder reflection, and the sign that keeps it stableA vector, the mirror line through the origin, and its reflection landing on the negative first axis at the same distance from the origin.the mirrorx, length 4.000Hx = (-4.000, 0)v = x − αe₁safe sign: α = −‖x‖, so v is formed from a sum and nothing cancelsunsafe sign: α = +‖x‖ gives ‖v‖ only 35.1% of ‖x‖ + ‖x‖ — the digits go‖HᵀH − I‖5·10⁻¹⁶‖Hx‖ − ‖x‖8.9·10⁻¹⁶second component2.2·10⁻¹⁶built from a unit vectororthogonality is structural
Fig. 15 An orthogonal transformation and what it costs in accuracy, from the field that measured it. The rotations that produce the GMRES number are of this kind, which is why their contribution is bounded.

The cost of asking, once

There is one number this pair of essays has not priced and it is the one a reader would act on, so it is worth putting down.

Recomputing ‖b − Ax‖ from the returned vector is one matrix–vector product. On a run of a hundred and fifty iterations that is two thirds of one per cent of the work, paid once, at the end, after the iteration has stopped. It does not change the answer, it does not change the search directions, and it cannot fail.

What it buys is that the number printed beside an answer is a number about that answer. On the conjugate gradient runs measured here the difference between the two is a factor of 7.3·10¹⁰; on almost every ordinary run it is a factor of one, and the product was wasted. That asymmetry is the whole argument: a check that is free when it is unnecessary and decisive when it is not is a check worth making unconditionally, which is the same reasoning the collection’s own rule about drawing a residual beside a factorisation rests on.

And a third quantity, which is neither

There is a number GMRES does not report and could not, and it is worth naming so the essay does not read as an endorsement.

GMRES’s residual is honest about the vector it would return. It says nothing about the error, and the collection’s standing warning applies here exactly as it does everywhere else: a relative residual of 10⁻¹⁵ on the Jordan block is a statement about ‖b − Ax‖ and not about ‖x − x*‖, and on a matrix with a condition number of 10¹⁰ the second can be ten orders larger than the first.

So the ledger for a GMRES run has three quantities in it, and the essays in this field have now separated all three:

  • what the method printed — honest here, and not in the conjugate gradient case;
  • the residual of the answer — which is what the printed number should equal, and does;
  • the error of the answer — which neither quantity measures and which the condition number amplifies from the first.
What the solver reported, and what it achievedFinal relative residual and final relative error against the working precision, on a logarithmic vertical axis. The residual sits at the stopping tolerance at every precision — the method converged, by every test available to it. The error runs from 5.26·10⁻¹⁴ at 53 bits to 0.00804 at 8, tracking the unit roundoff, which is drawn beside it.614223038465410⁻¹³10⁻¹⁰10⁻⁷10⁻⁴working significand bitsrelative sizeerrorresidualuthe gap the solver cannot seeerror ÷ residual at 24 bits4.1·10⁵error ÷ residual at 16 bits3.6·10⁷error ÷ residual at 8 bits1.9·10¹⁰every convergence test passesand the answer is wrong in proportion to u
Fig. 16 The third quantity, from the essay that established it. A method that reports its residual honestly is still not reporting its error.
Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 17 And the amplifier between them, which is the problem’s property and not the method’s.

Which methods are on which side

The pair in this essay is two points, and it is worth saying what the rest of the family looks like, because the classification is usable rather than anecdotal.

Methods whose reported residual is carried. Conjugate gradients, the whole family of short recurrences built on it, and every method whose implementation updates a residual vector because it cannot afford a second product: BiCG, BiCGStab, CGS, and conjugate gradients on the normal equations. All of them have a residual that is the state of a recurrence, all of them are subject to the mechanism the recurrence essay measures, and the size of the effect is set by how far the iterates travel rather than by anything about the method.

Methods whose reported residual is re-derived. GMRES and its relatives, MINRES, and anything else that reports the residual of a projected problem obtained by an orthogonal factorisation. Their number is a backward-error statement about a nearby problem, and the loss of orthogonality moves which nearby problem rather than corrupting the number.

The classification is not about symmetry, it is not about storage, and it is not about which method is better. It is about a single implementation decision that is usually made for cost reasons and is not usually written down as an accuracy decision at all.

And the useful consequence is the one a reader can act on: if a solver’s reported residual comes from a short recurrence, the number it prints is a claim about the run and not about the answer. One extra product at the end converts the claim; nothing else needs to change.

Steps to a relative error of 10⁻⁶, on one 60×12 problem at five condition numbersTwo curves against the condition number on a logarithmic horizontal axis, with four seeds drawn at each. The two recurrences compute the same iterates in exact arithmetic. At κ = 10² they cost 16 and 16 steps; at κ = 10¹⁰ they cost 110 and 209.10²10⁴10⁶10⁸10¹⁰0285684112140168196224condition numberstepsnormal equationsbidiagonalisationone sequence, two costssteps at κ = 10², both16at κ = 10⁶, ratio1.1at κ = 10¹⁰, ratio1.9the same iterates in the algebraand twice the work at κ = 10¹⁰
Fig. 18 The two-sided recurrence, from the essay that measured what a short recurrence costs. It is on the carried side of this classification, and for the same reason.
BiCG iterations on one 12×12 system against the distance of the shadow vector from a breakdownThe matrix, the right-hand side and the answer are the same at every stop. The only thing that moves is r̃₀, the second starting vector, which the method requires and for which every account gives the same non-reason. At η = 0.01 from the surface where the second divisor vanishes, BiCG is the direct method it is advertised as and finishes in 12 steps on 12 unknowns. The steps then run 12, 12, 13, 17, 20, 24, 80, 80, 80, 80 as η falls, and at 10⁻¹¹ the method has not converged after 80. Wherever it does finish it finishes at the same accuracy — the cost is the guarantee, not the answer.10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹0122436486072distance of the shadow vector from the breakdownBiCG steps to 10⁻¹⁵n = 12, where it should endone matrix, one right-hand sideρ₂ ÷ η, at every stop0.28steps at the far stop12steps at the near stop80residual history over n steps0.27the answer does not movethe guarantee does
Fig. 19 And a choice that changes the arithmetic and not the algebra, from the same essay. The carried-versus- re-derived decision is another of those, made in the implementation and not in the method.
Conjugate gradients at κ = 104, against the bound κ permitsA semi-logarithmic plot of the relative A-norm error against iteration count. The measured curve falls below a smooth dashed curve showing the classical condition-number bound.0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417
Fig. 20 The rate the condition number predicts, on the family that is on the carried side.

What a reader should take from the pair

The two essays together are a claim about implementations rather than about methods, and it is worth separating what a reader can act on from what is merely interesting.

What is interesting is the mechanism: a residual carried in a vector accumulates every rounding that enters it, and a residual re-derived from an orthogonal factorisation does not. That is a fact about arithmetic and it does not depend on anybody’s code.

What is actionable is narrower. If a solver reports a residual and the method behind it uses a short recurrence, the number it prints is a claim about the run rather than about the answer, and one extra product converts it into the second. If the method holds a full basis and derives the number from a projection, the printed number is already a statement about the vector that will be returned, and there is nothing to do.

Neither of those requires knowing the theory. Both require knowing which method is underneath, which is the one thing a caller usually does know and the one thing a residual printed on its own does not say. The collection’s rule about drawing a factorisation with its residual attached is the same idea one level down: a number and the thing it describes should not be separable by a reader skimming, and here they are separated by an implementation decision nobody wrote down as one.

The refusal

The claim under test is the over-general one: that a Krylov method’s reported residual cannot be trusted.

It is fed the conjugate gradient run on the same matrix, where the factor between the reported number and the answer’s is 7.3·10¹⁰. Fed that, the assertion that GMRES’s worst factor is under five is false about the wrong method, and the paired assertion — that on the same matrix the recurrence is a million times further out than the projection ever is — passes.

The two assertions together are the essay. Neither of them alone says anything: “GMRES is accurate” is a reassurance, and “conjugate gradients is not” is a scare. The measurement is the ratio between them on one matrix, and it is 7.3·10¹⁰ against 1.80.

What GMRES reports and what its answer's residual is, on a matrix with one eigenvalue at 1e-14GMRES also prints a residual it never computes from its answer — but it does not carry one. The number falls out of the Givens rotations at each step as the residual of the projected least-squares problem, re-derived from the whole basis rather than updated from the last value. Over 40 steps the two never separate by more than a factor of 1.80, and they do that while ‖VᵀV − I‖ for the Arnoldi basis is 1.41 — which is to say the basis has stopped being orthogonal and the residual computed from it is still honest.051015202530354010⁻³10⁻¹10¹GMRES steprelative residualthe two curves are the same curvea number re-derived, not carriedworst reported/actual factor1.8at step40‖VᵀV − I‖ of the basis1.4reported at the last step0.041actual at the last step0.074the same family of methodsand only one of them lies
Fig. 21 The pair on one matrix, once more. Everything in this essay is a reading of this figure beside the one in the essay before it.
How far the two residuals drift apart, against how large the iterates got on the waySix runs, differing only in the smallest eigenvalue of the matrix, which sets how far the iterates travel: from 9309 to 9.31·10¹³, ten decades of it. The upper line is the standard bound, u·‖A‖·max‖x‖, which every run stays inside. The lower series is what the drift actually does, and its fitted slope is 0.507 where the bound's is 1 — so the roundings are a random walk rather than a sum, the bound is correct, and it is loose by a square root that reaches 9·10⁴ across the sweep.10³10⁶10⁹10¹²10¹⁵10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1largest iterate on the way, ‖x‖difference between the two residualsthe bound, linear in ‖x‖what it does, slope 0.51a bound of one, a walk of a halffitted slope0.51the bound's slope1smallest gap measured4.6·10⁻¹⁵largest gap measured5.1·10⁻¹⁰share of the bound, near end6.6·10⁻⁴share of the bound, far end7.3·10⁻⁹the bound is rightand loose by a square root
Fig. 22 And the mechanism the other method has and this one does not: an accumulation whose size is set by how far the iterates travelled.
‖Aᵏ‖ for a 6×6 matrix whose spectral radius is 0.8, with 2 above the diagonalTwo curves against the power. The norm of Aᵏ rises to 19800 at step 24 before turning over and decaying to 2.5·10⁻⁵ by step 160; ρᵏ = 0.8ᵏ, drawn beside it, falls from the start. The two horizontal lines are the Kreiss bracket: the peak is at least K = 6757 and at most e·n·K = 1.1·10⁵, both computed from the resolvent norms outside the unit circle and not from the powers at all.027548110813510⁻⁶10⁻⁴10⁻²110²10⁴10⁶power‖Aᵏ‖Kreiss constant 6760e · n · K‖Aᵏ‖ρᵏtwo routes to one peakspectral radius0.8peak of ‖Aᵏ‖2·10⁴Kreiss constant6757e · n · K1.1·10⁵everything here decays in the endand one of these curves says how much first
Fig. 23 Non-normality in its own right, from the field that measures it — the property that makes the bidiagonal matrix here a hard case rather than an arbitrary one.
Adding 1 to 256, a thousand times, at 24 significand bitsThree lines against the number of additions: the exact total rising steadily, a staircase following it, and a flat line at the starting value.02505007501000250500750100012501500additionsrunning totalexactstochasticnearesta thousand additionshalf an ulp at 2561.5·10⁻⁵moves, round to nearest1000moves, stochastic1000relative error, nearest0relative error, stochastic024 significand bits, unbounded exponenta flat line is not a small error
Fig. 24 And an iteration that stops improving, from the arithmetic field. GMRES on the Jordan block does not stagnate, which is the last thing this essay checks and the least surprising.
Implied orthogonality of four factorisations of a 256×6 matrixFour curves against the condition number, both axes logarithmic. Cholesky QR's implied orthogonality has a fitted slope of 1.96 and reaches 0.169; the same routine run twice has a slope of 0.79 and reaches 2.32·10⁻¹⁰, which is where the Householder sweep and the reduction tree are.10²10³10⁴10⁵10⁶10⁷10⁸10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²condition number‖QᵀQ − I‖one passsweeptwicetreewhat the second pass removesfitted slope, one pass2fitted slope, two passes0.79rounds, two passes6rounds, the sweep18one pass squares the condition numberand two do not
Fig. 25 What a second orthogonalisation pass buys, from the field that priced it. GMRES does not take one and this essay is why its stopping test survives that.

Where this goes next

01428425670849811210⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best6.9·10⁻²¹the answer's, at its best5.1·10⁻¹⁰unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10¹³iterations drawn110the recurrence remembers every roundingand the stopping test is written in it The residual the method reports Conjugate gradients prints a relative residual of 6.9·10⁻²¹. The unit roundoff is 1.1·10⁻¹⁶, so that is not a small residual and not a large one — it is not a residual. The vector the method is holding at that step has ‖b − Ax‖/‖b‖ = 5.1·10⁻¹⁰, and nothing in the run says so. 02468101210⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1iteration‖r‖ / ‖b‖Laplaciancyclic shiftno progress at allevery eigenvalue of the shift is on the unit circleand it predicts nothing The spectrum that predicts nothing For a symmetric matrix the eigenvalues govern how fast an iteration converges. Drop symmetry and they stop governing anything — there is a matrix whose eigenvalues are as evenly spread as eigenvalues can be, on which GMRES makes no progress at all until the last possible step. 05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here An orthogonalisation nobody calls one Conjugate gradients are derived as a minimisation and behave as an orthogonalisation, which is why the finite-termination property in every textbook is not a property the method has in floating point.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Arnoldi iterationBackward errorGivens rotationGMRESKrylov subspaceModified Gram–SchmidtNon-normalityOrthogonalityStopping criterion