The zero that means it is finished

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.

What the zero says #

Write 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ₖ₊₁,ₖ, the norm of the part of A·vₖ that is not already in Kₖ.

If that number is zero then A·vₖ is in Kₖ. Every earlier basis vector has the same property by construction, so A·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ₖ, subtracts off ⟨vᵢ, w⟩·vᵢ 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ᵢ,ₖ of the Hessenberg matrix, and they can be anything at all — large, small, zero — without anything being wrong. Only hₖ₊₁,ₖ, 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.

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

Four sizes, thirty-two spectra, and not one exception. At n = 20, 30, 40, 50 and 60 every one of the eight spectra breaks at exactly m, and the worst residual at the breaking step is 6.38·10⁻¹⁶, 8.51·10⁻¹⁶ and 9.1·10⁻¹⁶ at the three intermediate sizes — the rounding level, drifting upwards only as fast as the arithmetic does.

What does move with n is the size of the jump, and it moves the wrong way. The smallest gain the breaking step shows over the step before it is 1.46·10¹⁰ at n = 20, 8.75·10⁹ at n = 40 and 3.4·10⁹ at n = 50 — a factor of four smaller across a factor of 2.5 in the size. The signal is enormous at every one of them, so nothing here is at risk; but it is shrinking, and it shrinks because the residual before the breaking step is itself falling as the larger space does more work. So the cleanest possible termination test — watch for the residual to drop by ten orders in one step — is a test whose margin depends on n, and the direction of the dependence is the one that makes it worse on the large problems where an early exit is worth the most.

This is why every account of Krylov convergence talks about the spectrum rather than about the dimension, and why the rate the condition number predicts has to be measured against a method that does not terminate when the theory says it will. 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 approach, which is a straight line #

A repeated eigenvalue arriving twice from a Krylov recurrence is its own subject; an exactly repeated one in the matrix 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 subdiagonal h₆,₅ 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.

How long “for a while” is #

Close eigenvalues behave, for a while, as though they were identical — the dip says how deep, and running past it says how long. The same family, carried to convergence:

The while is one inequality. Five steps suffice exactly while the tolerance is above the dip, and the dip is 0.176·η — so a cluster of width η is a point for any tolerance above about 0.18η and is not one below it. Every row obeys that: at η = 10⁻⁶ a tolerance of 10⁻⁶ sits above 1.8·10⁻⁷ and costs five steps, while 10⁻⁹ sits below it and costs ten.

Which is a decision procedure rather than an intuition. A caller who knows the cluster width and the tolerance knows before running anything whether the spectrum has five eigenvalues or thirty, and that is not something watching the residual would tell them — the residual after five steps looks like convergence in every row of the table — and a stopping test reading it would stop, which is a stopping test is a race’s hazard arriving from the spectrum’s side rather than the machine’s.

What happens past the threshold is the other half of the answer, and it is gradual. The count does not jump to n; it climbs, and it reaches n only when the cluster is about as wide as the gaps between clusters. At η = 10⁻¹ — a tenth of the spacing — reaching 10⁻¹² takes 28 of the 30 steps, so the clustering is worth nothing there and was worth everything at 10⁻⁸. The continuum runs the whole way from m eigenvalues to n, and the tolerance is what says where on it a particular solve sits.

That also puts a number on why preconditioning is worth what it is worth. Changing the condition number on purpose is usually argued from κ and the √κ rate; the table above argues it from the other end. A preconditioner that gathers a spectrum into five clusters of width 10⁻⁸ has bought a five-step solve at every tolerance down to 10⁻⁹, and one that gathers it into five clusters of width 10⁻¹ has bought almost nothing. The width of the clusters is the quantity, and it is not κ. The same distinction appears one field over as rank is a decision: a spectrum with no gap has no rank, and a spectrum whose clusters are wider than the tolerance has no reduced step count, for the same reason in both cases — the structure the argument needs is below the resolution the arithmetic has.

And where a cluster is genuinely too wide, the repair is not more steps but a different space: restarting is a filter is what a method does when it cannot afford the steps the spectrum demands, and the filter it applies is chosen against exactly the cluster widths measured here.

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.

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.

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₆,₅ 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.

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ₖ₊₁,ₖ 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.

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.

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.

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ₖ₊₁,ₖ 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 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ₖ₊₁,ₖ 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.

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.

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.