Keeping the vectors, and losing the bound

Restarting is a filter measured what a restart does: multiply the starting vector by Π(A − θⱼI), suppressing the directions the run has decided it does not want. The Ritz pairs are discarded at the end of every cycle and rebuilt from scratch in the next one, and the whole content of the method is in the polynomial.

There is another way to do it, and it does the opposite. Keep the k Ritz vectors that are wanted, add the last residual vector, and continue the Lanczos recurrence from there. Nothing is filtered and nothing is rebuilt.

The two are known to be equivalent, and the equivalence is the reason thick restarting is usually described as an implementation detail rather than as a method: same answer, simpler code, no shifts to choose. Everything below takes that seriously enough to measure it, and the first two measurements confirm it exactly. The third does not.

The identity, and the shape it produces #

After m steps, each Ritz vector satisfies

A yᵢ = θᵢ yᵢ + sᵢ v, with sᵢ = βₘ times the last entry of yᵢ’s coefficient vector

so each retained vector is an eigenvector of A up to one shared direction — the residual vector v, which is the same vector for all of them. That is what makes the retention work: the space spanned by {y₁ … y_k, v} is almost invariant, and A maps each yᵢ back into it exactly.

It also fixes the shape of the projected matrix. The retained values sit on the diagonal, their couplings to v sit in one row and one column, and the ordinary Lanczos tridiagonal continues from there.

That is an arrowhead, and this site has an essay about arrowhead matrices in a completely different field: two ends of the same arrow is about what happens when one row’s position decides whether elimination costs n or n³. The same shape arrives here from a restart rather than from an ordering, and for the same underlying reason — one object coupled to everything, and everything else coupled to nothing.

What it buys #

Two things, both measured against restart.js’s filtered restart on the same matrix with the same seed.

The same eigenvalues. At the first two cycles the two methods return Ritz values agreeing to 3.6·10⁻¹⁵ and 1.2·10⁻¹⁴; over six cycles the worst disagreement is 4.8·10⁻¹⁴. That is the equivalence theorem, measured rather than quoted, and the assertion is deliberately placed at the early cycles where the claim is about the algebra rather than about the rounding.

For a third of the products with A. The filtered restart rebuilds the whole m-step recurrence each cycle and then applies p shifts, so it spends m + p products; thick restarting keeps k of the vectors and spends p. Over eight cycles with k = p = 4 that is 36 products against 96.

The identity holds for one restart and not for two #

Now the part that is not in the description.

The residual of a retained vector is supposed to lie entirely along v. Measure the part of it that does not:

The third cycle is the first one whose retained vectors were themselves built out of retained vectors. Before it the identity is exact to rounding; from it the identity is wrong by nine orders of magnitude more than rounding, and stays wrong.

The consequence is that the arrowhead is not the Rayleigh quotient of the basis it names. ‖T − QᵀAQ‖/‖QᵀAQ‖ is 1.6·10⁻¹⁶ and 1.1·10⁻¹⁵ at the first two cycles and 4.2·10⁻⁶ at the third.

That last point is worth dwelling on. Every other failure of a Lanczos-type method on this site has been a lost basis — ghosts, duplicated eigenvalues, orthogonality collapsing geometrically over ten steps. This is not that. The basis here is reorthogonalised twice against everything at every step and is orthonormal to 7·10⁻¹⁶ throughout. What is wrong is the matrix written down about it.

So the stopping rule stops measuring #

The residual bound is βₘ|s_{m,i}|, it comes out of the eigenvectors of the projected matrix, and it costs nothing. It is what a real code stops on, and it is computed inside a model that has stopped describing the operator.

Forty-four orders of magnitude, in the direction that makes the run look finished.

And nothing else reports it. The eigenvalues are correct to 4.3·10⁻¹⁴ at every one of those cycles, which is the accuracy the method is usually asked about; the basis is orthonormal; the residual bound is falling monotonically, which is what a converging run looks like. The only quantity that says anything is the one nobody computes, because computing it costs a product with A per vector and the bound exists to avoid exactly that.

The vectors really are worse #

A residual is a proxy, and a proxy that has just been shown to be untrustworthy is a poor thing to rest a finding on. These matrices are built from a spectrum and an orthogonal basis, so the true eigenvectors are available — the exact ground truth move the site makes wherever it can — and the angle to them can be measured directly.

A factor of 410 in the quantity the method exists to produce, and nothing at all in the quantity anybody checks.

That last column is not a coincidence, and the depth phase already measured why. An eigenvalue is quadratically insensitive to an error in its eigenvector — the Rayleigh quotient of a vector 8.6·10⁻⁶ from an eigenvector is wrong by something like the square of that, divided by a gap — so an eigenvalue error of 10⁻¹⁴ is exactly what a vector error of 10⁻⁶ produces.

What is stored, which does not change #

It is worth being explicit that this is not a memory trade.

The filtered restart keeps m = k + p basis vectors and rebuilds them each cycle. Thick restarting keeps m too: the k retained Ritz vectors, the residual vector, and the p − 1 it generates. Neither method’s storage grows with the number of cycles, which is the property both exist for and which the previous essay asserts of the filtered one.

What differs is entirely in the arithmetic: one method spends m + p products with A a cycle to rebuild what it discarded, and the other spends p because it did not discard it. The retained vectors are not free to keep — they have to stay orthonormal, and they are reorthogonalised against at every step — but that cost is in inner products rather than in products with A, which on a large sparse matrix is the cheap side of the ledger.

One product with A a cycle #

The border entries are the couplings yᵢᵀAv. The cheap method takes them from the previous cycle’s bounds, which is what the identity says they are. Recomputing them takes one product Av, shared by all k retained vectors, and then k inner products.

Seven orders of magnitude of eigenvector, for seven products with A. The arrowhead becomes the Rayleigh quotient again at every cycle, the bound tracks the residual to a part in a thousand while both are falling, and the residual reaches the level everything else on this site reaches.

The eigenvalues are 4.26·10⁻¹⁴ and 4.09·10⁻¹⁴ — the same to rounding either way, which is the last piece of why the defect is invisible.

How big it is depends on a tuning parameter #

The slider is the number of retained vectors, and the gap between the reported bound and the residual at the tenth cycle runs across it like this:

At every setting the bound leaves the residual behind. What the setting decides is how much damage that does: at three retained vectors the cheap run happens to land exactly where the repaired one lands, and at four it is a billion times worse.

A defect whose size varies this way across a tuning parameter is one nobody finds by testing at one setting, which is the practical reason this essay exists. The assertion is therefore written against the part that is invariant — the bound is the residual at the first cycle, and is at least a hundred times below it by the tenth, at every k — rather than against the number the default happens to produce.

How it was found, which is the argument for measuring twice #

The path is worth writing down because no single measurement would have got there.

The first thing drawn was the bound against the cycle count, because that is the stopping rule and the obvious thing to plot. It falls smoothly through forty decades, which looks like a method converging beautifully and is the reason nobody looks further.

The second was the residual recomputed from the returned pair — one product with A per vector, which is exactly the cost the bound exists to avoid, and therefore a thing only a figure would ever pay for. That produced the flat line at 5.7·10⁻⁵ and the question.

The third measurement was the diagnosis, and it had two candidates. This site’s standing explanation for a Krylov method behaving badly is a lost basis, so ‖QᵀQ − I‖ was measured first and came back at 7·10⁻¹⁶ — ruling out the answer that would have been reached for by habit. What was left was the projected matrix itself, and comparing it against the Rayleigh quotient of its own basis is a three-line computation that no working code would ever do, because forming QᵀAQ costs m products with A and defeats the entire purpose of a projection method.

Every step in that chain costs something the method is designed not to spend. That is the general shape of this site’s findings and it is worth naming: a defect that only an expensive measurement can see will survive in a method whose whole point is to avoid that expense.

When to use it anyway #

None of the above says not to use thick restarting. It says what to measure.

The method’s advantages are real and are not affected by any of this: the same eigenvalues as the filtered restart, a third of the products with A, no shift selection, no bulge chase, and a projected matrix whose shape a reader can check by eye. For a run of two or three cycles — which is what a well-separated spectrum needs — the arrowhead identity has not yet broken and the free bound is a bound.

What the measurements say is: recompute the border, at one product with A a cycle, and the free bound becomes trustworthy again along with everything else. That is between two and three per cent of the products a cycle costs at any realistic p, and it removes the only quantity in the method that lies.

If that is too expensive, the fallback is the one this essay’s third measurement used: recompute the residual of the returned pairs when the run stops, once, and compare it against the bound that stopped it. One product with A per wanted eigenvector, paid once, and the whole failure mode is visible.

What this is an instance of #

The site has a name for this shape by now.

mixedcg.js measured a preconditioned conjugate gradient run whose reported residual read 1.1·10⁻¹³ while the true error was 1.6·10⁻⁷, and drew the rule out of it: a quantity that steers may be rounded, and a quantity that measures may not. The projected residual bound here is the same kind of object — a quantity computed inside the method’s own model of itself, used as though it were a measurement of the world.

The difference, and it is the reason this one is worse, is that the mixed-precision case has a knob with the trade written on it. Nobody chooses the stale border; it is what the derivation says the number is, and the derivation is right in exact arithmetic. There is no setting at which a reader is warned.

What the drag does #

The slider is the number of retained vectors, and every position of it shows the same qualitative picture with a different magnitude — which is the point made above. Two things do not change at any position: the reported bound falls without limit, and the residual it names does not.

What is worth watching is the third curve. The repaired run’s residual is never the larger of the two, at any setting, and the assertion says exactly that rather than claiming an improvement — because at three retained vectors there is no improvement to claim, and a figure that asserted one would be asserting something false at one of its own stops.