The answer that arrives when the space runs out

One picture in this collection carries two different events, and it has been read as carrying one. It plots ‖QᵀQ − I‖ for two Krylov bases against the number of vectors held, with the distance from the dominant eigenvalue of a damped chain drawn on the same axes, and the reading taken from it is that the cheap recurrence loses its basis while its answer goes on improving.

That reading is right, and it is the finding a Krylov space for a problem that is not linear was written to make. What it does not say is that the improving curve has a second discontinuity in it, that the second one is not a property of the recurrence at all, and that on the chain length the picture is drawn for it lands just outside the frame.

The two events are worth separating before either is described. The first is a loss of orthogonality: the second-order recurrence orthogonalises a residual whose auxiliary half has been carried along uncorrected, so the stored coordinates decay, and they decay after about the same number of vectors whatever the size of the problem. The second is a loss of room: the second-order basis holds vectors of length n, so once it holds n of them there is nothing left in ℝⁿ for the next one to be. The quantity that decides the first is the number of steps. The quantity that decides the second is the length of the chain, and no figure in this collection had ever varied it against the number of steps.

The number that settles which is which is the one that picture never moves. Both of its curves are drawn against the number of vectors held, and the chain’s length enters only as the size of the problem underneath, fixed at forty in every drawing of it this collection had. Measured at twelve, sixteen, twenty and twenty-four masses instead, one of the two events moves with the chain and the other stays where it is — and which of them moves decides what a machine-precision reading at the end of a long run is a measurement of.

The second cliff is at the chain’s own length, and it is off this page #

Twenty vectors is where the drawing stops, and forty masses is the chain it is drawn on. Carried past the edge of the frame on that same chain, the distance to the dominant eigenvalue reads 5.24·10⁻⁵ at thirty-eight vectors, 4.71·10⁻⁵ at thirty-nine, and 1.90·10⁻¹⁵ at forty. One step, ten orders. Nothing about the curve up to thirty-nine anticipates it: over the three vectors before the drop the error falls by a factor of 2.3, from 1.08·10⁻⁴ at thirty-six, which is the same unhurried order as the drawn part of the figure, where sixteen more vectors buy a factor of nineteen.

Forty is the number of masses. It is also the length of every vector the second-order recurrence stores, and therefore the number of them that can be linearly independent. At m = 40 the recurrence holds forty vectors in a forty-dimensional space, and if they are independent their span is the whole of it — at which point the subspace the method is projecting onto is not an approximation of anywhere, it is everywhere, and the small problem it forms is the large problem in different coordinates.

That is a claim about arithmetic as well as about algebra, so it is measured rather than argued, and it is measured at four chain lengths rather than at the one where it was noticed. The generator draws twenty vectors, so a chain short enough for its own exhaustion point to fall inside that range puts both events on one page.

What the projection is once the space has run out #

The mechanism is a determinant, and it is short enough to write down.

The second-order route projects the problem rather than the matrix: from a basis Q it forms Mₘ = QᵀMQ, Cₘ = QᵀCQ and Kₘ = QᵀKQ, and asks for the eigenvalues of that m × m quadratic, which is the step a matrix that depends on its own eigenvalue supplies the machinery for. When m = n the matrix Q is square, and for any λ

det(λ²Mₘ + λCₘ + Kₘ) = det(Qᵀ(λ²M + λC + K)Q) = det(Q)² · det(λ²M + λC + K).

The two polynomials differ by the constant det(Q)², so they have the same roots — all 2n of them, not merely the dominant one. The projected quadratic at m = n is not a good approximation to the original; it is the original with a factor in front of its characteristic polynomial.

The condition for that is det(Q) ≠ 0. It is not orthonormality, not near-orthonormality, and not any bound on ‖QᵀQ − I‖. A basis whose columns are badly skewed still spans, and spanning is the whole of what the identity needs, which is the distinction a basis that is the same subspace and not the same thing is about, arriving here on the side where it costs nothing.

Measured, the whole spectrum does arrive at once. At m = n the projected quadratic’s 2n Ritz values match the chain’s closed form to a worst relative distance, taken over every exact eigenvalue, of 8.91·10⁻¹⁴ at twelve masses, 2.83·10⁻¹⁴ at sixteen, 5.42·10⁻¹³ at twenty and 4.12·10⁻¹² at twenty-four. One step earlier, at m = n − 1, the same quantity reads 0.551, 0.573, 1.04 and 0.227. The transition is one step wide and it takes the entire spectrum with it, which is what distinguishes a space running out from a Ritz value converging: a converging Ritz value is one number getting better, and this is 2n numbers becoming right together.

Termination that survives a basis of 1.39, where another method’s does not #

That last pairing is the sharpest number on the page and it is worth stating on its own. At sixteen masses, the step at which the second-order route returns the dominant eigenvalue to 5.96·10⁻¹⁶ is the step at which its own basis reads ‖QᵀQ − I‖ = 1.390. For scale, orthogonal is a number fixes what that means: zero is exactly orthonormal, and about one is what an arbitrary set of unit vectors gives. The coordinates are useless as a frame at the moment the answer becomes exact.

This collection has the opposite result for the method most readers reach for first. An orthogonalisation nobody calls one measures conjugate gradients run to exactly n steps on the model problem at n = 30 and finds the residual basis at ‖RᵀR − I‖ = 1.41 — the same neighbourhood as 1.390 — and there the finite-termination property is not degraded but absent: on the matrix built at κ = 10⁶ where the property is actually under pressure, the loss reads 5.29 and the residual after all thirty steps is 5.1·10⁻⁴.

Two methods, two nearly equal readings of the same measurement, opposite consequences. The reason is that the two terminations are derived from different things. Conjugate gradients terminate in n steps because the residuals are mutually orthogonal, and an argument that needs orthogonality is spent exactly when orthogonality is. The second-order route terminates because Q is square and invertible, and rounding does not change the rank of a square matrix whose smallest singular value is well clear of zero. Measured at m = n, that singular value is 1.000 at twelve masses, 0.1305 at sixteen, 0.03319 at twenty and 0.01318 at twenty-four. Skewed, and nowhere near singular.

The two cliffs cross at sixteen masses #

With four chain lengths measured, the relation between the two events can be stated rather than guessed at.

The orthogonality cliff barely moves. Taking the first subspace size at which ‖QᵀQ − I‖ exceeds one, it is thirteen vectors at twelve masses, sixteen at sixteen, seventeen at twenty and sixteen at twenty-four — a spread of four steps across a factor of two in problem size, which is the finding that essay already reports as a collapse arriving at about the same count whatever the chain’s length.

The exhaustion cliff moves with the chain: twelve, sixteen, twenty, twenty-four, exactly. So the two cross. At twelve masses the space runs out first and the answer is exact while the basis is still orthogonal to five digits. At sixteen they arrive on the same step. From twenty masses upward the basis fails first and the exactness comes later — three steps later at twenty, eight at twenty-four, twenty-three at forty. Above the crossing, every reading of exactness is taken through a basis that is no longer one, which is why the two have looked like a single event: on the chain they were first drawn for, one always follows the other.

Arnoldi falls off the same cliff at twice the step #

The comparison method exhausts its space too, and the size of its space is the whole point of the storage claim.

Arnoldi on the linearisation holds vectors of length 2n, so its Krylov space is all of ℝ²ⁿ at m = 2n, and its own plunge is there. Measured on the four chains, the distance from the dominant eigenvalue at m = 2n is 4.77·10⁻¹⁶, 5.96·10⁻¹⁶, 4.17·10⁻¹⁵ and 1.19·10⁻¹⁵. At m = n, the step where the second-order route has finished, the same runs read 7.45·10⁻³, 3.46·10⁻⁴, 3.55·10⁻³ and 8.53·10⁻⁴ — two to four digits short, at twice the storage.

The cliff is cleanest on the shortest chain, where ordinary convergence has not had time to get close first: at twelve masses the linearised route reads 1.73·10⁻⁶ at twenty-two vectors, 9.53·10⁻⁹ at twenty-three and 4.77·10⁻¹⁶ at twenty-four. On the longer chains the approach is gentler, because 2n steps of a Krylov method on a 2n-dimensional problem is enough for genuine convergence as well, and the two effects arrive together.

So the storage saving and the termination point are one fact rather than two. The factor of two in vector length is the factor of two in steps to exactness, and it has to be, because both are the dimension of the space each method’s vectors live in.

The reading that has to be corrected #

The essay that introduced the recurrence carries an extended table to sixty-four vectors on the chain of forty, and one row of it needs a different explanation from the one it was given.

That row is forty-four vectors, where the second-order route reads 2.38·10⁻¹⁶ against a basis whose ‖QᵀQ − I‖ is 2.83, and the sentence built on it is that past the collapse the cheap recurrence does not merely survive but wins by ten orders. Forty-four is past forty, and forty is the length of the vectors it stores. The plunge on that chain is at forty, measured above at 1.90·10⁻¹⁵ against 4.71·10⁻⁵ one step earlier, so by forty-four the second-order span has been the whole of ℝ⁴⁰ for four steps and its projected quadratic has been the chain’s own quadratic for four steps. The number is real. The cause named for it is not: what the row records is a method that has run out of room, not a method that has converged.

The sizes that table samples are the reason the distinction was invisible. They step from thirty-eight to forty-four, so the exhaustion point falls between two rows and the first reading taken after it is the first reading that could have been anything else. The assertion the site runs behind that table requires the machine-precision reading to arrive at or before fifty vectors, and that requirement is met by exhaustion and by nothing else — every sampled size below forty is at 10⁻⁵ or worse, and the two that table prints above forty, forty-four and sixty-four, read 2.38·10⁻¹⁶ and 8.69·10⁻¹⁵.

The comparison in the same row needs the same correction. At sixty-four vectors the linearised route is 2.41·10⁻⁶ from the eigenvalue, which is genuinely ten orders behind — and sixty-four is sixteen short of the eighty dimensions that route’s vectors have. It is a method with a fifth of its space still unused, measured against one that finished twenty-four steps earlier. Read at any subspace size below forty the ordering is a comparison of two unfinished runs, and that essay’s own caution that the ordering is not monotone applies exactly there.

What is not corrected is most of that essay, and it is worth being specific. The storage factor of two and the Ritz-value factor of four are exact by construction and stand. The comparison at six to twenty-six vectors, which is the essay’s own table and its own headline claim, sits entirely below the exhaustion point on a chain of forty and is unaffected — 1.4·10⁻² against 7.3·10⁻² at six vectors is a measurement of two unfinished runs at equal work, which is what it says it is. The orthogonality collapse stands, at 1.347 by twenty vectors, and so does the observation that no quantity the method returns reports it. What changes is one row’s explanation and the sentence drawn from it.

What it costs, and where it stops being true #

Exactness at m = n is a statement about the algebra, and the arithmetic charges for it.

The price is the conditioning of Q. At m = n its largest singular value is 1.00 at twelve masses and 1.41 at the other three, and its smallest is the number given above, so the condition number is 1.00 at twelve masses, 10.8 at sixteen, 42.6 at twenty and 107 at twenty-four — and the accuracy of the exact answer follows it, at 8.91·10⁻¹⁴, 2.83·10⁻¹⁴, 5.42·10⁻¹³ and 4.12·10⁻¹² over the whole spectrum. At twenty-four masses the square of the condition number against the unit roundoff is about 1.3·10⁻¹², the same order as the 4.12·10⁻¹² measured, which is what a congruence by an ill-conditioned Q would be expected to cost. So the termination is exact and the answer is not, and the gap grows with the chain.

The larger limit is the one that decides whether any of this is useful. A damped structure of the kind this field’s methods exist for has n in the hundreds of thousands, so exhaustion is unreachable by many orders and the second cliff is a fact about small problems only. That is precisely why the correction matters rather than being a detail: a machine-precision reading attributed to convergence is a claim that scales, and one attributed to exhaustion is a claim that does not. On a real problem the useful subspace sizes are the ones a shift puts within a handful of steps, which is well below the collapse and further below the chain’s length, and the rate that governs them is the unhurried one the drawn part of every figure here shows — two digits in six vectors, the same rate a rate that does not notice the size insists is the quantity to report rather than an iteration count.

None of which is an argument against the recurrence, and it is worth saying where the honest comparison sits. The rate the condition number predicts makes the general case that a guarantee stated for a limit is a poor guide to the sizes anybody runs, and finite termination is the extreme instance of it: a property that holds at exactly one subspace size, reachable on a chain of twenty-four and unreachable on a structure of two hundred thousand. What a caller gets on the large problem is the measured rate and nothing else, and that rate is flat in the size of the problem — at six vectors the distance to the dominant eigenvalue is 0.01043 at twelve masses, 0.01305 at twenty-four and 0.01379 at forty.

And the recurrence does not announce the moment it happens. It carries a breakdown test: if the new residual’s norm falls below 10⁻¹³ the run stops, because a zero residual is what a Krylov space closing looks like. In exact arithmetic the residual at step n + 1 is exactly zero, since a vector orthogonalised against an orthonormal basis of the whole space is the zero vector. Measured, that norm is 0.317 at sixteen masses and 0.468 at twenty-four, against a threshold of 10⁻¹³. The test never fires, the recurrence goes on producing vectors, and it produces them past the point where a seventeenth vector of length sixteen can be independent of the previous ones — the smallest singular value of Q at m = n + 1 is 2.94·10⁻¹⁵⁷. The loss of orthogonality is what disarms the detector: the sweep subtracts components along vectors that are no longer perpendicular, so what it returns is not the projection onto anything and its norm is not a measurement of what is left. The same defect hides both the exhaustion and the signal that would have reported it.

What follows for a code that uses this recurrence #

A stopping test cannot be the answer’s own improvement. The distance to the dominant eigenvalue falls smoothly for thirty-nine steps and then falls ten orders in one, so no extrapolation from the history predicts the step it arrives at, and no plateau precedes it. The residual the method reports is the general form of this: a quantity computed from the projected problem is consistent with the projected problem, and the projected problem is the thing that changed.

A restart resets the counter that matters. Restarting is a filter keeps the subspace small deliberately, and small means below n by a wide margin, so a restarted run never reaches exhaustion and never sees the plunge. The discipline is right and the reason to know what it forgoes is that a short chain used as a test problem will show the plunge and a restarted production run will not, which is a difference in the problem rather than in the code.

A bound computed from the projection inherits the projection’s standing. Keeping the vectors, and losing the bound measures a residual bound that reaches 9.4·10⁻⁴¹ while the residual it bounds sits at 5.7·10⁻⁵. At m = n the projected quadratic is the problem, so its residuals are the problem’s residuals and the bound is finally honest — for exactly the one subspace size nobody can afford.

Ritz vectors are still not orthogonal. Nothing here repairs the coordinates, so everything that essay says fails silently downstream of a collapsed basis still fails: repeated Ritz values where the spectrum is not repeated, of the kind an eigenvalue that arrives twice measures for Lanczos, and a projection that is not a projection. Exactness of the eigenvalues at m = n is a statement about roots of a polynomial and carries nothing about the vectors.

And the subspace was never the thing that was damaged. The span of the second-order basis is the top block of the linearised Krylov space by construction, and that stays true whatever the stored coordinates do — the same separation the plane survives what its vectors do not makes between a subspace and any basis for it. Here it is what makes the whole result possible: at m = n the span is all of ℝⁿ because it cannot be anything else, and reading that as a method converging is reading a property of the room as a property of the search.