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

The zero that means it is finished

Every Krylov method ends by dividing by a number the previous step produced, and when that number is zero the recurrence stops. In Arnoldi the stop is the answer — the subspace has closed, the solution is inside it, and the residual is at the unit roundoff. The literature calls it a lucky breakdown, and the adjective is doing real work.
31 min read 28 figures Where the algebra stops being the arithmeticMeasured, not assertedTwo routes to a numberAssertions that reject

Worth reading first: The rate the condition number predicts · The spectrum that predicts nothing · An eigenvalue that arrives twice.

Conjugate gradients, GMRES, MINRES, LSQR, Lanczos, Arnoldi. Every one of them builds a basis by taking a vector, multiplying it by A, subtracting off what is already in the space, and then dividing by the norm of what is left. The last step is the one this essay is about, because the norm of what is left can be zero.

When it is, there is no next basis vector, the recurrence has nowhere to go, and the method stops one step short of the step it was about to take. A first reading says a division by zero has been narrowly avoided and something has gone wrong. The opposite is true, and the size of the opposite is worth stating precisely: at that moment the answer is exact.

The step at which the Arnoldi recurrence breaks down, against the number of distinct eigenvalues, n = 30Every matrix here is 30×30 and every one of them breaks down before step 18. The step is the number of distinct eigenvalues in the spectrum, exactly, at every one of the 8 spectra drawn — the points lie on the diagonal and the horizontal line at n = 30, which is the bound every course states, is off the top of the picture. The relative residual at the breakdown step is at most 5.64·10⁻¹⁶; one step earlier it is at least 3.5·10¹⁰ times larger.024681012141618024681012141618distinct eigenvalues in the spectrumstep the recurrence stops atthe step is m, not nn = 30 throughoutspectra drawn8every one breaking at m8worst residual at the breakdown5.6·10⁻¹⁶smallest gain over the step before3.5·10¹⁰an invariant subspace contains the answerand its dimension is what the method costs
Fig. 1 Thirty-by-thirty matrices, and the step at which the recurrence stops. The bound every course states is the horizontal line at n = 30, and it is off the top of this picture at every one of the eight spectra drawn.

What the zero says

Write K_k for the space spanned by b, Ab, …, A^(k−1)b. The Arnoldi recurrence produces an orthonormal basis for it, and the number it divides by at step k is h_{k+1,k}, the norm of the part of A·v_k that is not already in K_k.

If that number is zero then A·v_k is in K_k. Every earlier basis vector has the same property by construction, so A·K_k ⊆ K_k: the space is invariant. And an invariant subspace containing b contains A⁻¹b, because A restricted to it is an invertible map of the space to itself. The solution is inside the space the method has already built.

So GMRES, which minimises the residual over exactly that space, finds a residual of zero. Not approximately zero — zero in the algebra, and at the unit roundoff in the arithmetic. On the eight spectra in the figure above the largest relative residual at the stopping step is 5.6·10⁻¹⁶, and the smallest gain over the step before it is a factor of 3.5·10¹⁰.

That last number is the one worth holding onto. The step before the breakdown is not nearly finished. It is 10¹⁰ away.

Why the division is there at all

It is worth being precise about which division this is, because a Krylov method contains several and only one of them has this property.

The recurrence at step k computes w = A·v_k, subtracts off ⟨v_i, w⟩·v_i for every i ≤ k, and then normalises what is left. The subtractions are the projections; the normalisation is the division. The numbers subtracted are the entries h_{i,k} of the Hessenberg matrix, and they can be anything at all — large, small, zero — without anything being wrong. Only h_{k+1,k}, the norm of the remainder, is a norm, and a norm is zero exactly when the vector is.

That is what makes the Arnoldi breakdown unambiguous. There is nothing to interpret: the vector being normalised is either the zero vector or it is not, and if it is, the space is invariant. No other reading is available. The next essay is about a recurrence whose divisor is an inner product of two different vectors rather than a norm of one, and every ambiguity in it comes from that single change.

The other divisions in a Krylov solver — the α of conjugate gradients, the step length of a line search, the entry a triangular solve divides by — are divisions by quantities that can be small without meaning anything, and this collection has essays about two of them already.

Two Krylov methods against products with A, at a kernel shift of 0Two error curves against the number of products with A, on a logarithmic vertical axis. The Arnoldi method reaches 0.1477 after 3 products and is 98.9 by the end of the run. The bidiagonal method reaches 0.1426 after 40 and degrades far more slowly.16111621263136414651566110⁻¹110¹10²products with Arelative errorArnoldi's best: 3Arnoldibidiagonalwhat a step buysArnoldi's best0.15products to reach it3bidiagonal's best0.14products to reach it40a tenth of the work to the same answerand no time at all spent there
Fig. 2 One of those other divisions, seen where it matters: the point at which a Krylov least-squares method turns from fitting the answer to fitting the noise.

Where it happens, which is not where the textbook says

The bound everybody learns is that GMRES converges in at most n steps, because K_n is the whole space. It is true, it is almost never useful, and there is a sharper statement that costs nothing extra.

The Krylov space closes as soon as it contains a full set of eigenvectors for the eigenvalues b actually has components along. If A has m distinct eigenvalues, that happens at step m — whatever n is. The minimal polynomial of A restricted to the space b generates has degree m, and the degree of that polynomial is the dimension at which the space stops growing.

The hero above is that claim tested rather than restated. Each point is a 30×30 matrix built to have exactly m distinct eigenvalues, m running from two to sixteen, and every point sits on the diagonal. The recurrence stops at m. At m = 2 it stops after two steps on a thirty-dimensional problem.

The step at which the Arnoldi recurrence breaks down, against the number of distinct eigenvalues, n = 60Every matrix here is 60×60 and every one of them breaks down before step 18. The step is the number of distinct eigenvalues in the spectrum, exactly, at every one of the 8 spectra drawn — the points lie on the diagonal and the horizontal line at n = 60, which is the bound every course states, is off the top of the picture. The relative residual at the breakdown step is at most 8.34·10⁻¹⁶; one step earlier it is at least 1.6·10¹⁰ times larger.024681012141618024681012141618distinct eigenvalues in the spectrumstep the recurrence stops atthe step is m, not nn = 60 throughoutspectra drawn8every one breaking at m8worst residual at the breakdown8.3·10⁻¹⁶smallest gain over the step before1.6·10¹⁰an invariant subspace contains the answerand its dimension is what the method costs
Fig. 3 And at n = 60 the picture is the same picture. The size of the matrix is the bound that moves and does not matter; the number of distinct eigenvalues is the one that does.

This is why every account of Krylov convergence talks about the spectrum rather than about the dimension. A matrix with a thousand unknowns and six distinct eigenvalues is a six-step problem. It is also why clustering is worth so much: eigenvalues that are not identical but are close behave, for a while, as though they were.

The spectrum of C⁻¹T against T's own, n = 16, ρ = 0.8Two sorted spectra on a linear vertical axis. The unpreconditioned Toeplitz matrix spreads its eigenvalues across its whole range; the preconditioned operator puts 13% of them within a tenth of 1, with 3 below zero.0481216-3-2-101234567891011121314151617181920212223242526index, sortedeigenvalueT aloneC⁻¹Teigenvalues, not singular valuesfraction within 0.1 of 10.13eigenvalues below zero3λ_min of the preconditioner-0.075the dashed line is 1, where the cluster formsthe solid line is zero
Fig. 4 A spectrum that clusters rather than repeats. The essays on structure reach this from the other side — a preconditioner is worth having when it produces a picture like this one.
GMRES on the Laplacian and on the cyclic shift, both 12×12A 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.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
Fig. 5 The counterweight, and it is one this collection has already made: a spectrum alone does not predict what GMRES will do on a nonsymmetric matrix. What the breakdown argument gives is a bound on the number of steps, not a rate for the steps before it.

The approach, which is a straight line

An exactly repeated eigenvalue is a measure-zero event. Every matrix anybody hands a solver has n distinct eigenvalues, so the breakdown above never happens and the whole argument is about a case that does not occur.

That reading is wrong in a way that is measurable, and the measurement is the second half of this essay. Take the same construction and replace each repeated eigenvalue with a cluster of width η. At η = 0 the spectrum has m distinct values; at any η > 0 it has n. Nothing closes, nothing is invariant, and the breakdown does not happen. What happens instead is a dip, and the dip is exactly proportional to η.

The Arnoldi subdiagonal and the GMRES residual against the width of an eigenvalue cluster, n = 30, 5 clustersA matrix of size 30 whose spectrum is 5 clusters of width η. At η = 0 the Krylov space closes at dimension 5 and the recurrence breaks down there. As η opens up, h_{6,5} rises with a fitted slope of 1.00 and the relative residual at step 5 with a slope of 1.00 — both exactly proportional to η, at a fixed ratio of 102 between them. The residual at step 4, one short of the breakdown, is 0.0238 and stays inside a per cent of that across the whole sweep: the step before the event carries no information about it.10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10²cluster width ηmagnitudeh at step 5residual at 5residual at 4fitted over ten decadesslope of h against η1slope of the residual1h ÷ residual, at every stop102residual at step 40.024the zero arrives in proportionand the step before it learns nothing
Fig. 6 Ten decades of cluster width, and three curves. Two of them have slope one. The third does not move.

The subdiagonal h_{6,5} is 18·η, to three digits, at every stop from 10⁻⁴ down to 10⁻¹². The GMRES relative residual at step 5 is 0.176·η over the same range. Their ratio is 102 and it is the same 102 at every stop.

So the breakdown is not a special case sitting at the end of a continuum — it is the continuum’s limit, and the continuum is a straight line of slope one. A cluster of width 10⁻⁸ gives a residual of 1.8·10⁻⁹ after five steps of a thirty-dimensional problem. That is not a breakdown and it is not different from one in any way a caller would notice.

And the third curve is the finding

The dashed line in that figure is the residual at step 4 — one short of the dip. It is 2.38·10⁻² and it stays inside a per cent of 2.38·10⁻² across the whole sweep.

Ten decades of change in the quantity the method is about to meet, and the observable a caller is watching does not move.

That is a shape this collection keeps meeting from different directions. A residual history is what every stopping rule reads, every progress bar plots and every convergence claim is made from, and it is silent about the structure that is about to end the iteration. The information arrives entirely in the last step.

The Arnoldi subdiagonal and the GMRES residual against the width of an eigenvalue cluster, n = 30, 2 clustersA matrix of size 30 whose spectrum is 2 clusters of width η. At η = 0 the Krylov space closes at dimension 2 and the recurrence breaks down there. As η opens up, h_{3,2} rises with a fitted slope of 1.00 and the relative residual at step 2 with a slope of 1.00 — both exactly proportional to η, at a fixed ratio of 4.08 between them. The residual at step 1, one short of the breakdown, is 0.3299 and stays inside a per cent of that across the whole sweep: the step before the event carries no information about it.10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10²cluster width ηmagnitudeh at step 2residual at 2residual at 1fitted over ten decadesslope of h against η1slope of the residual1h ÷ residual, at every stop4.1residual at step 10.33the zero arrives in proportionand the step before it learns nothing
Fig. 7 At two clusters the same two slopes and the same flat line, with the flat line higher: a two-dimensional Krylov space has removed less of b before the last step than a five-dimensional one has.
The Arnoldi subdiagonal and the GMRES residual against the width of an eigenvalue cluster, n = 30, 12 clustersA matrix of size 30 whose spectrum is 12 clusters of width η. At η = 0 the Krylov space closes at dimension 12 and the recurrence breaks down there. As η opens up, h_{13,12} rises with a fitted slope of 1.00 and the relative residual at step 12 with a slope of 1.00 — both exactly proportional to η, at a fixed ratio of 82000 between them. The residual at step 11, one short of the breakdown, is 7.576·10⁻⁵ and stays inside a per cent of that across the whole sweep: the step before the event carries no information about it.10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10²cluster width ηmagnitudeh at step 12residual at 12residual at 11fitted over ten decadesslope of h against η1slope of the residual1h ÷ residual, at every stop8.2·10⁴residual at step 117.6·10⁻⁵the zero arrives in proportionand the step before it learns nothing
Fig. 8 And at twelve, lower. The number of clusters moves where the dip is and how much has been removed before it, and moves neither slope.

What a cluster is worth, stated as a number

The proportionality gives a usable rule, and it is worth writing down in the form a practitioner would use it.

Suppose a preconditioner has been applied and the spectrum of the preconditioned matrix falls into m groups, each of relative width η. Then after m steps the relative residual is about 0.18·η — a fifth of the cluster width — and getting further requires either a better preconditioner or many more steps, because the next m steps buy the next factor of the same size rather than the next order of magnitude.

Two things follow that are not obvious from the more usual statement about eigenvalue clustering.

The number of clusters costs steps and the tightness of them costs digits, and they are not interchangeable. Halving the number of clusters halves the work; halving their width buys a factor of two in the answer for no work at all. A preconditioner that turns twelve loose clusters into five loose ones and a preconditioner that turns twelve loose clusters into twelve tight ones are doing different things, and a summary that reports only “the spectrum is better clustered” cannot tell them apart.

And the constant is not universal. The 0.18 above is for this family; the ratio between the subdiagonal and the residual — 102 here — is a property of how the clusters are arranged and of where b points. What is universal is the slope: both quantities are proportional to η, at every stop, on every construction tried. A slope of one is a statement about the mechanism; the constant in front of it is a statement about the matrix.

Two filters on one sum, λ = 0.001The weight each term of the solution is given, against its index. Truncation is a step: one for the first 32 terms and zero after. Tikhonov is σ²/(σ² + λ²), which falls smoothly through the same place. The unregularised solution is the constant one, which is why it divides noise by a σ of 1.7·10⁻¹³.081624324048566400.250.50.751index kfilter factor fₖno regularisation: fₖ = 1truncationTikhonovthe same sum, three weightsTikhonov, relative error0.16truncation, relative error0.13no filter at all5.5·10⁸both filters are one expression with a different weightfₖ = 1 is the catastrophe
Fig. 9 The same relationship read as a filter: which components of the answer a k-step method has recovered and which it has not, which is what a Krylov space is doing when it closes.

What the arithmetic does to all of this

Everything above is exact-arithmetic reasoning, and this site’s premise is that exact-arithmetic reasoning is where the interesting part starts rather than ends.

Two things happen in floating point and they pull in opposite directions.

The zero is not zero. At η = 0 the measured h_{6,5} is about 10⁻¹⁴ rather than 0 — it is u·‖A‖·√n, which is what an orthogonalisation leaves behind. So a program cannot test for zero; it has to test for small, and small compared to what is a decision somebody makes. The last essay of this collection’s current run is about that decision and about the fact that both of its directions are failures.

And orthogonality is lost long before the space closes. The Arnoldi basis is orthogonal in the algebra and drifts in the arithmetic, and the drift arrives on a schedule this site has already measured: it begins exactly when the first Ritz value converges.

‖QᵀQ − I‖ at every step of a Lanczos run, n = 24The 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.0481216202410⁻¹⁷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. 10 The loss of orthogonality in a Lanczos basis. It is not a slow decay; it is a collapse, and it starts where the first eigenvalue arrives.
Copies of each eigenvalue after 80 steps on a 40×40 matrix with a simple spectrumA bar per eigenvalue that came back more than once, showing how many times. The matrix has 40 distinct eigenvalues by construction; the run returned 25 extra copies of 13 of them, the most-copied arriving 5 times. Every copy is accurate to 1.9·10⁻⁸ relative, which is why nothing but the true spectrum could detect them.λ = 105 timesλ = 9.55 timesλ = 95 timesλ = 8.55 timesλ = 2.952 timesλ = 2.92 timeseigenvalues that arrived more than once — the matrix has 40 distinct onesa spectrum with the wrong multiplicitiesextra copies, no reorthogonalisation25extra copies, full reorthogonalisation0worst relative error among the copies1.9·10⁻⁸steps taken of 80 asked for, full40no arithmetic error was madeevery one of these is right to eight digits
Fig. 11 And the consequence a caller sees: eigenvalues that arrive twice. A method whose basis has stopped being a basis returns copies, and every copy is accurate.
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. 12 What repairing it costs, which is less than it looks, because the repair only has to touch the directions the theory names.

The two effects together are why the measured picture needs a reorthogonalisation pass to be about what it claims to be about. Without one, the h_{k+1,k} at a near-breakdown is the loss of orthogonality rather than the quantity being measured, and the twelve-decade line bends at about 10⁻⁸ — which is √u, and a fact about Gram–Schmidt rather than about breakdown.

What it is worth, in the methods that use it

The lucky breakdown is not a curiosity. Three things in daily use depend on it.

Deflation. When a subdiagonal entry becomes negligible the problem splits, and the algorithm continues on the smaller block. That is the same event as the one here, read as a licence rather than as a stop, and it is what makes the QR algorithm cost what it costs.

The width of a block of 4, at every stepThree step functions against the step count. A block of 4 independent columns keeps its width at every step. A block with one column copied from another loses that column before the first product with A and keeps 3 for the rest of the run. A block whose columns all lie in the span of two eigenvectors starts at 2 and ends at 0 after one step, having found the whole subspace.01234567891001234stepcolumns in the blockindependent columnsone column repeatedinside an invariant planewhen a block narrowsrandom: width at the end4repeated: width at the end3invariant: vectors found2the mechanism worksand nothing a random start does sets it off
Fig. 13 A block method’s deflation: a direction that has become numerically dependent on what is already there is dropped, and the block shrinks.

Restarting. A restarted Arnoldi keeps the part of the space that has converged and throws away the rest, which is only sound because the converged part is very nearly invariant.

What 3 shifts do to each direction of the starting vectorAmplification of each eigen-component against its eigenvalue, on a logarithmic vertical axis. The points are measured — the component of the vector along each eigenvector, before and after — and the curve is |Π(λ − θ)| with the roots at the 3 discarded Ritz values. They agree to 3.9·10⁻¹¹. The 3 wanted directions are amplified by at least 1 and everything else by at most 0.802.1357910⁻⁵10⁻⁴10⁻³10⁻²10⁻¹1eigenvalue λamplificationkeptdiscardeda filter on the starting vectormeasured against the polynomial3.9·10⁻¹¹worst kept direction1best discarded direction0.8the roots are the discarded Ritz valuesand the vertical lines are where they sit
Fig. 14 A restart is a polynomial filter on the starting vector, and what it is filtering towards is an invariant subspace.
12 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 1.08·10⁻¹³ across 12 cycles and 140 products, and the true error reaches 7.11·10⁻¹⁵. The basis is 8 vectors at every cycle and never grows.0183654729010812610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹products with Asizeresidual boundtrue errorbounded memorybasis vectors kept8products with A140worst error in the k wanted7.1·10⁻¹⁵the bound is free and the error is notand the basis never grows
Fig. 15 And what the filtering buys, cycle by cycle.
The reported bound and the residual it bounds, over 10 cyclesThree quantities against the cycle count on a logarithmic vertical axis. The residual bound the method reports falls without limit, reaching 1.73·10⁻³⁰. The residual it claims to bound stops at 2.44·10⁻¹⁴ and does not move. Recomputing the arrowhead's border entries, at one extra product with A a cycle, takes the residual to 2.48·10⁻¹⁴.1234567891010⁻⁵⁵10⁻⁴⁹10⁻⁴³10⁻³⁷10⁻³¹10⁻²⁵10⁻¹⁹10⁻¹³10⁻⁷10⁻¹cyclesizethe residualrepairedthe reported boundwhat the stopping rule readsreported at the last cycle1.7·10⁻³⁰the residual there2.4·10⁻¹⁴with the border recomputed2.5·10⁻¹⁴products, cheap and repaired9a bound with nothing under itand one product a cycle to fix it
Fig. 16 Keeping more of the previous space than the theory strictly licenses, and the bound that goes with it.

Preconditioning. The whole point of a preconditioner is to make the spectrum look like the picture at the top of this essay — a few clusters instead of a spread — which is a statement about how close to a lucky breakdown the recurrence can be brought.

The condition number a sketched preconditioner leaves, against the condition number it was givenTwo curves against κ(A), both axes logarithmic. The matrix's own condition number climbs the diagonal from 100 to 10¹⁰; κ(AR⁻¹), where R comes from a QR of a 4n-row sketch, is 2.2284 at every one of them — the same number to ten digits, not a similar one. The reason is two lines of algebra: with G = SU the preconditioned singular values are those of (GᵀG)⁻¹, which has no spectrum of A in it at all.10²10⁴10⁶10⁸10¹⁰110²10⁴10⁶10⁸10¹⁰condition number of the matrixcondition number seen by the iterationκ(A), unchangedκ(AR⁻¹)a bound with no κ(A) in itκ(AR⁻¹), every κ(A)2.2κ(SU), the other route2.2κ(A), across the sweep10⁸κ(AR⁻¹), across the sweep1the sketch never sees the spectrumand the spectrum cancels out of the answer
Fig. 17 A preconditioner changing the spectrum it is applied to, which is the same axis this essay’s cluster width sits on.
LSQR steps to a normal residual of 10⁻¹², with and without a sketched preconditionerTwo curves against κ(A). The unpreconditioned solve takes 15, 22, 32, 59, 150 steps as the conditioning worsens, and at the worst one it does not reach the tolerance inside 150. The preconditioned solve takes 10, 10, 10, 10, 10 — flat across eight decades. Both return the same vector, to the accuracy the problem allows: the preconditioned answer and the direct one differ by about κ(A)·u at every κ, 5.7·10⁻⁶ at the worst. The SAME sketch used as an answer instead of as a preconditioner is 4.5 away.10²10⁴10⁶10⁸10¹⁰110¹10²10³condition number of the matrixLSQR steps to 10⁻¹²no preconditionersketchedthe cost, not the answersteps, preconditioned10steps, without150answer gap, preconditioned5.7·10⁻⁶κ(A)·u at the worst κ1.1·10⁻⁶answer gap, sketch-and-solve4.5one sketch, two things to do with itand only one of them moves the answer
Fig. 18 And the steps that follow from it.
Conjugate gradients at κ = 103, 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κ bound85 steps40×40, spectrum spread evenly in logbound permits 448
Fig. 19 The oldest member of the family, on a well-conditioned problem. Conjugate gradients has the same breakdown available to it and reaches it for the same reason.
The residual basis of conjugate gradients, at κ = 106A 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. 20 And the same method on an ill-conditioned one, where its residual vectors stop being orthogonal — the failure that the reorthogonalisation above exists to prevent.

Two routes to the same step

The stopping step is claimed above in two ways that share no arithmetic, which is this collection’s standing requirement for a number.

The first is the recurrence itself: run Arnoldi and see where h_{k+1,k} falls to the roundoff. The second is the construction: the matrix was built with a prescribed spectrum, so the number of distinct eigenvalues is known before anything runs. The two agree at every one of the eight spectra and at both sizes.

They could disagree. If the right-hand side happened to have no component along one eigenvector, the recurrence would stop earlier than the spectrum predicts, and that would be a fact about b rather than about A. The construction avoids it by giving b a component along everything, and the agreement is therefore evidence rather than a tautology.

The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 21 The habit this rests on, in its strongest form elsewhere on this site: a problem whose answer is known in closed form, so an error is measured rather than estimated.
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: κ · u30×30, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 22 And the reason a residual of 10⁻¹⁶ is not the end of the argument: a small residual is not a small error, and the gap between them belongs to the problem.

The refusal

An assertion that cannot fail proves nothing, so the machinery behind this essay is fed the case it must refuse: an Arnoldi recurrence started from the zero vector.

It is not a contrived input. A restarted method whose filter has removed everything, a matrix-free operator returning zero, a right-hand side that is exactly in the null space of a preconditioner — all of them produce it. What makes it worth refusing rather than tolerating is what the code does otherwise: the first basis vector becomes 0/0, every subsequent step orthogonalises a vector of NaNs against a basis of NaNs, and h_{k+1,k} comes back as NaN, which is not greater than any tolerance and therefore reads as a breakdown at step one.

A recurrence that reports finishing immediately, on a problem it has not touched, is the failure mode this whole essay is about arriving in its least useful form.

Past the cliff: what each method returnsLoss of orthogonality at four condition numbers past κ²u = 1. Cholesky QR returns a factor whose implied Q is 0.37 away from orthogonal at κ = 10⁸ and 1.49 at 10¹² — and at 3·10¹⁰ it refuses outright, on a pivot of -6.5·10⁻¹⁷. The sweep and the tree return a usable factorisation of every one of these matrices.‖QᵀQ − I‖ of the implied Qκ = 10⁸, Cholesky0.37κ = 10⁸, sweep8.5·10⁻⁹κ = 10¹⁰, Cholesky1.3κ = 10¹⁰, sweep3·10⁻⁷κ = 3·10¹⁰, Choleskyrefusedκ = 3·10¹⁰, sweep6.7·10⁻⁶κ = 10¹², Cholesky1.5κ = 10¹², sweep1.4·10⁻⁴the safe run is the one that failsrefusals in the range1wholly non-orthogonal returns2the pivot it refused on-6.5·10⁻¹⁷one of these outcomes is safeand it is the refusal
Fig. 23 The general shape, from the essay on communication: a routine that cannot do what was asked has two options, and one of them is worse than an error.

And what it means for a stopping rule

A stopping rule reads the residual and decides whether to continue. Everything above says that on a problem with clustered eigenvalues the residual history is nearly flat for m − 1 steps and then falls by ten orders in one.

A rule that stops when progress has stalled therefore stops immediately before the answer, on exactly the problems a preconditioner has been applied to. That is not a hypothetical: “stagnation” detection is a standard feature of iterative solver libraries, it is usually implemented as a comparison of consecutive residuals, and the picture above is what it is looking at.

The repair is not subtle and it is the reason GMRES implementations do what they do: the bound the method has on its own residual is computable at every step from the Hessenberg matrix, without forming the iterate at all, and it says nothing about stagnation because it is a bound rather than a history. A rule built on it stops when the bound is small rather than when the history is flat, and the difference between the two is the difference between finishing and stopping one step short.

Adding 0.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.02505007501000250275300325350375additionsrunning totalexactstochasticnearesta thousand additionshalf an ulp at 2561.5·10⁻⁵moves, round to nearest1000moves, stochastic1000relative error, nearest1.7·10⁻⁵relative error, stochastic8.6·10⁻⁷24 significand bits, unbounded exponenta flat line is not a small error
Fig. 24 Stagnation of a different kind, from the arithmetic field: a sum that stops moving because the addend has fallen below the ulp of the running total. The shape a history makes when it stalls is not evidence about why.
Conjugate gradients on an ill-posed problem at 0.10% noiseTwo curves against the step count on a logarithmic vertical axis. The relative residual falls at every one of the 120 steps without exception. The error against the true signal falls to 0.1050 at step 44 and then climbs, reaching 0.56 by the end — 5.3 times its best value.015304560759010512010⁻³10⁻²10⁻¹1steprelative sizeleast error: 44discrepancy stop: 27errorresidualthe knob is an integerleast error, at step44error there0.11error at step 1200.56the residual falls at every stepthe error turns and keeps rising
Fig. 25 And a third: an iteration whose residual keeps falling while its error turns round. A stopping rule reads the first curve and is asked about the second.
The spectrum of C⁻¹T against T's own, n = 64, ρ = 0.8Two sorted spectra on a linear vertical axis. The unpreconditioned Toeplitz matrix spreads its eigenvalues across its whole range; the preconditioned operator puts 97% of them within a tenth of 1, with 0 below zero.08162432404856640123456789index, sortedeigenvalueT aloneC⁻¹Teigenvalues, not singular valuesfraction within 0.1 of 10.97eigenvalues below zero0λ_min of the preconditioner0.11the dashed line is 1, where the cluster formsthe solid line is zero
Fig. 26 And the same clustering at four times the size, where the number of clusters has not changed and the number of eigenvalues has.
The step at which the Arnoldi recurrence breaks down, against the number of distinct eigenvalues, n = 45Every matrix here is 45×45 and every one of them breaks down before step 18. The step is the number of distinct eigenvalues in the spectrum, exactly, at every one of the 8 spectra drawn — the points lie on the diagonal and the horizontal line at n = 45, which is the bound every course states, is off the top of the picture. The relative residual at the breakdown step is at most 9.65·10⁻¹⁶; one step earlier it is at least 1.7·10¹⁰ times larger.024681012141618024681012141618distinct eigenvalues in the spectrumstep the recurrence stops atthe step is m, not nn = 45 throughoutspectra drawn8every one breaking at m8worst residual at the breakdown9.6·10⁻¹⁶smallest gain over the step before1.7·10¹⁰an invariant subspace contains the answerand its dimension is what the method costs
Fig. 27 At forty-five the bound the textbook gives has moved half as far again, and the step the recurrence stops at has not moved at all.

What is next

The zero here means an object was found. The next essay is about the same arithmetic event in a method that is only two lines different, where it means nothing was found at all — no invariant subspace, no converged anything, and no next term. The event is identical, the numbers look identical, and the two meanings are as far apart as two meanings get.

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. 28 The pair of recurrences that essay is about, from the side this collection has already looked at them: one sequence, two ways of generating it, and the difference is not in the algebra.

Where this goes next

0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417 The rate the condition number predicts Conjugate gradients converge at a rate governed by the square root of the condition number. That is a bound rather than an estimate, it is provable, and it is loose enough that provisioning iterations from it wastes nine out of ten. 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. λ = 105 timesλ = 9.55 timesλ = 95 timesλ = 8.55 timesλ = 2.952 timesλ = 2.92 timeseigenvalues that arrived more than once — the matrix has 40 distinct onesa spectrum with the wrong multiplicitiesextra copies, no reorthogonalisation25extra copies, full reorthogonalisation0worst relative error among the copies1.9·10⁻⁸steps taken of 80 asked for, full40no arithmetic error was madeevery one of these is right to eight digits An eigenvalue that arrives twice A matrix with forty distinct eigenvalues, handed to Lanczos for eighty steps, returns twenty-five extra copies of thirteen of them — the largest arriving five times. Every copy is accurate to 1.9·10⁻⁸ relative. No arithmetic error was made, nothing overflowed, and a caller counting eigenvalues gets the wrong multiplicity from a computation in which no individual number is wrong.

Reads more easily once this is understood

Essays that name this one as worth reading first.

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.

ArnoldiClustered spectrumConvergence rateFinite terminationGMRESInvariant subspaceKrylov subspaceLoss of orthogonalityResidual