A basis that is the same subspace and not the same thing

Exact at the points that were named builds a reduced model by projecting onto V = span{(σᵢI − A)⁻¹B}. The interpolation conditions that follow are conditions on V, not on the vectors that span it, so in exact arithmetic any basis will do and the choice is free.

It is not free, and the way it is not free is worth an essay because the rate at which the obvious basis fails is set by something the user chose without information.

The obvious basis #

If all r points are taken at one σ, the subspace is spanned by the moments

(σI − A)⁻¹B,  (σI − A)⁻²B,  …,  (σI − A)⁻ʳB,

which is what the derivation writes down, and what an implementation writes first. It is one solve per vector with the same factorisation reused, which makes it the cheapest possible way to obtain r vectors, and it corresponds to matching r + 1 derivatives of H at the single point σ rather than values at r different ones. That is a legitimate thing to want — it is what the series that has to be squared back is doing for a matrix function — and it is not what fails here.

Its condition number, measured on a twenty-four-state model at σ = 20:

vectors    2      3      4      5      6      7      8
κ₂        11    1.4e2  2.7e3  6.9e4  2.6e6  1.2e8  7.7e9

multiplying by between thirteen and sixty-six per vector. The same eight solves spent at eight points spread across the model’s poles and orthogonalised as they are built stay at κ₂ = 1.0.

σ = 20 is one shift, and the shift turns out to matter by five orders of magnitude.

The worst shift is fifty thousand times worse than the best, and the best is not the largest. Across σ = 0.5, 3, 120 and 800 the moment basis reaches 2.55·10¹², 7.08·10¹¹, 5.04·10⁷ and 3.51·10⁸ — falling by five orders and then rising again, with an interior optimum somewhere near a hundred.

The orthogonalised basis is at exactly 1 in all four, which is what makes this a comparison rather than a tuning exercise. There is a real optimum in the shift and finding it is worth five orders — and it buys a basis that is still unusable, because 5·10⁷ on eight vectors in double precision has already spent half the significand. So the shift is a parameter with a large effect and no useful setting, which is the strongest form the essay’s argument can take: the failure is not that σ was chosen badly, and choosing it well does not rescue the method.

That also explains why the number in the table above is not worth quoting precisely. κ₂ = 7.7·10⁹ at eight vectors is a fact about σ = 20; the same measurement at σ = 120 gives 5.04·10⁷ and at σ = 0.5 gives 2.55·10¹². What survives across all of them is the rate — between thirteen and sixty-six per vector — and the fact that the orthogonalised alternative does not move.

Why it collapses, and it is not the same reason as last time #

This site has an essay about a basis that loses independence and it is an orthogonalisation nobody calls one. The sequence there is v, Av, A²v, … and the reason is the power method: repeated application of A converges to the dominant eigenvector, so every new vector is nearly the last one, at a rate given by the ratio of the two largest eigenvalues.

The mechanism here is the same and the ratio is a different quantity, which is the point.

(σI − A)⁻¹ has eigenvalues 1/(σ − λₖ), so repeated application converges to the eigenvector whose λ is nearest σ. The rate is the ratio of the two smallest |σ − λ|. So the number that governs how fast the basis collapses is a ratio of distances from the interpolation point to the poles — and the interpolation point is the user’s choice, made before anything about the poles is known.

The user’s own choice decides how fast their basis becomes unusable. Measured across the range: at σ = 800, past every pole, every solve is nearly B/σ and the collapse is fastest; at σ = 3, near the bottom of the spectrum, it is slower; there is no placement at which eight moments at one point stay independent.

That is a sharper statement than the power-method version. There the rate is a property of the matrix and nothing can be done about it inside the method. Here it is a property of a parameter, and what can be done about it is not to reorthogonalise harder but to spend the solves somewhere else.

Spending them somewhere else #

Eight solves at eight different points span a different subspace — the rational Krylov space rather than the polynomial one at a single shift — and it is a better subspace as well as a better conditioned one.

Better conditioned, because the vectors are solves at points that are far apart, so nothing makes them parallel. Measured at κ₂ = 1.0 for all eight, which is what orthogonalisation performed as the basis is built gives when there is genuinely a direction to add each time.

Better as an approximation, because eight distinct interpolation points give sixteen conditions spread over the frequency range, where eight moments at one point give sixteen conditions all at the same place. The error curve in the previous essay has r zeros distributed over four decades; the single-point version has one zero of multiplicity 2r and rises monotonically away from it in both directions.

The cost is identical. Eight solves either way. The single-point version can reuse one factorisation across all eight, which on a dense model is a real saving and on a sparse one — where each solve is an iterative solve anyway — is not.

The orthogonalisation, and where it goes #

“Just orthogonalise the basis” is the natural response to the figure above, and it is worth measuring rather than arguing about, because the argument against it does not survive.

The argument was that by the eighth moment the eight vectors span a space that is numerically seven-dimensional, so a QR returns an orthogonal Q whose eighth column is rounding. The rank says otherwise. At eight vectors the moment basis’s smallest singular value is 1.3·10⁻¹⁰ of its largest — ill conditioned, and six orders clear of rank-deficient. Its numerical rank at 10⁻¹⁴ is eight, and stays full until twelve.

And the outcome says otherwise too. Building the two-sided projection three ways on the same model and measuring the interpolation error at σ:

A QR afterwards recovers the interpolation to the rounding level at every size — including twelve and fourteen, where the basis genuinely is rank deficient and κ has reached 1/u. The raw basis is four to six orders worse and stops improving past ten, which is the collapse the figure draws; but the collapse belongs to the basis, and the projection stops inheriting it the moment the basis is orthogonalised.

The reason the noise directions do not spoil the answer is that the interpolation conditions are conditions on the subspace. A subspace containing the first ten moments interpolates whether or not its eleventh direction means anything: the noise adds a direction the reduced model can ignore, and a projection onto an orthonormal basis has no ill-conditioned r × r inverse left to compute.

So of the two reasons this essay gives for spending the solves at spread points, the conditioning one does not survive and the approximation one does. Eight distinct points give sixteen conditions spread over the frequency range where eight moments give sixteen at one place, and that is a statement about which subspace, which no orthogonalisation can change. The rational basis is the better choice for the second reason, at a cost the essay attributes to the first.

The correction has a practical edge to it, because the two reasons recommend different things when the solves cannot be spent freely. A single factorisation reused across r moments is genuinely cheaper than r factorisations at r shifts on a dense model, and if conditioning were the objection then the cheap route would be unusable and the expensive one forced. It is not: the cheap route plus one QR — O(nr²), less than a single solve — gives a projection as good as any. What it does not give is a better subspace, and whether that matters depends on whether the model is wanted accurate near one frequency or across a range.

So the decision reduces to a modelling question rather than a numerical one, which is the useful place for a decision to be. Match many derivatives at one point when the application lives near that point; spread the points when it does not; and in either case orthogonalise, because the basis’s condition number stops being anybody’s problem the moment somebody does.

Orthogonalising as the basis is built is different only in that it fails loudly: the norm of the new vector after subtraction is a number the routine can look at, and a construction that finds it below a threshold refuses rather than continuing. The implementation behind these figures does that, and it is what stops the sixth interpolation point in the previous essay when two points are both past the last pole.

The subtraction is done twice, which is two Gram–Schmidts applied as a policy rather than as a comparison: one pass is not enough when the vectors are nearly parallel, and nearly parallel is the ordinary case here rather than the pathological one — the same policy doing it twice prices in the communication field, for a different reason and at the same cost.

What a projection needs that a spanning set does not #

The distinction the figure turns on is easy to state and easy to lose, so it is worth one paragraph of its own.

A spanning set is enough to define a subspace. A projection onto that subspace requires solving with the basis: the reduced matrices are WᵀAV, WᵀB and CV, and the normalisation that makes WᵀV the identity involves inverting an r × r matrix built from the two bases. When V has condition number 10⁹, that r × r matrix inherits it, and the reduced model’s coefficients are computed from a system that has lost nine digits before the model exists.

So the interpolation conditions still hold for the subspace, and the model that comes out is not the model the subspace defines. Measured on this basis at eight vectors, the biorthogonality ‖WᵀV − I‖ that ought to be 10⁻¹⁶ reaches 10⁻⁷ — which is exactly the number of digits the condition number predicts, and which no assertion about spans would have caught.

This is the same distinction the plane survives what its vectors do not draws in the spectral field, arriving at the opposite conclusion because the question is different. There the subspace is the answer and the vectors are incidental, so a collapsing basis costs nothing. Here the vectors are the arithmetic, so it costs everything.

The refusal, and what it is worth #

The implementation refuses a new vector whose norm after subtraction is below 10⁻¹³ of the right-hand side’s. That threshold is a decision and it is worth saying what it is a decision about.

It is not a tolerance on accuracy. It is the point at which the direction being added is indistinguishable from the rounding in the subtraction that produced it — so accepting it would produce a basis vector whose content is the arithmetic error of the previous seven. Everything downstream would then be computed from a matrix one of whose columns means nothing, and nothing downstream would notice: the projection would succeed, the reduced model would have the right size, and its behaviour would depend on the last bits of the last solve.

That is why the check is a refusal rather than a warning. This site’s habit is that an assertion which has never rejected anything proves nothing, and this one rejects: it is what stops the six-point construction in exact at the points that were named when two points are both past the last pole, and it is the reason that essay could report the failure as a property of the point placement rather than as a mysteriously poor model.

The measurement that is not available #

There is a temptation to compare the two bases by the accuracy of the reduced models they produce, and it is worth saying why that comparison cannot be made cleanly.

They do not produce the same reduced model. The single-point basis and the spread basis span different subspaces, so they interpolate at different places, so their errors differ for a reason that has nothing to do with conditioning. Any accuracy comparison between them is a comparison of two interpolation strategies with a conditioning difference confounded into it.

What can be measured cleanly is the conditioning, and that is what the figure shows. The conclusion it supports is narrow and correct: the obvious basis cannot be projected with past about six vectors, whatever it spans. Whether the subspace it spans was the right one is a separate question with a separate answer, and the two should not be run together.

That is the same discipline six routes to one spectrum insists on for the polynomial field — vary one thing, hold the rest, and refuse the comparison that varies two.

What survives from this into the rest of the field #

Two things.

Every method in this field is built out of shifted solves, and how they are spent is the decision. Eight at one point, eight at eight points, or — as interpolating at the model’s own poles does — eight at points chosen by the reduced model itself. The cost is the same in each case and the outcomes are not comparable in accuracy without care, but one of the three cannot be computed at all past six vectors and that narrows the field to two.

A rate that looks like a property of the problem can be a property of a parameter. The power-method rate in the Krylov field is fixed by the matrix. The same rate here is fixed by a number the user typed, and the repair is to change the number rather than to work harder on the arithmetic. It is worth checking, whenever a basis collapses, which of the two situations it is — because only one of them has a cheap way out.

The one number to look at #

If there is a single line to carry out of this essay it is that the quantity worth printing is κ₂ of the basis, and almost nothing prints it.

A reduced model is returned as three small matrices. Nothing about them says how they were computed or from what. The conditioning of the basis is available for the cost of one small SVD at the moment the basis is finished, it is the number that says whether the model means anything, and it is thrown away by every interface that returns only the model.

That is the same complaint orthogonal is a number makes about Q: the adjective is free and the measurement is nearly free, and only one of them can fail. Here the measurement costs r³ flops on an r that is single digits, against a construction that cost r sparse solves, so it is free by any accounting anybody would do.

The same failure, three fields apart #

It is worth putting the three instances of this collapse side by side, because the differences between them are more instructive than the similarity.

The Krylov sequence v, Av, A²v. The rate is |λ₂/λ₁|, a property of the matrix. Nothing inside the method changes it, and the repair — orthogonalise as the basis is built — is a repair to the arithmetic and not to the subspace. The method still converges; it just needs its basis maintained.

The moment sequence at one shift. The rate is a ratio of distances from σ to the poles, a property of a parameter. Orthogonalising helps in the same way, and there is a second repair available that the first case does not have: change the parameter. Spending the same solves at different points gives both a better-conditioned basis and a better approximation, which is why every implementation does it.

The second-order recurrence for a quadratic eigenvalue problem. Measured in a Krylov space for a problem that is not linear: ‖QᵀQ − I‖ runs 5.9·10⁻¹¹, 6.1·10⁻⁷, 1.6·10⁻³, 1.35 at m = 6, 10, 14, 20, while the linearised basis built by the same modified Gram–Schmidt sweep stays at 2.4·10⁻¹³. There the rate is a property of the recurrence rather than of the matrix or of a parameter, and neither repair is available: the recurrence is what buys the storage saving, and changing it gives back the saving.

Three mechanisms, one symptom, and three different amounts of recourse. What the comparison suggests as a habit is to ask, on meeting a collapsing basis, which of the three it is — because the answer decides whether the fix is free, cheap, or unavailable, and the symptom alone does not distinguish them.

The other thing all three share is that the answer often keeps improving while the basis collapses. That is the trap in each case: a run watching only its residual reports progress at every step, and the object it is building has stopped being a basis. Watching κ₂ costs nothing and is the only thing that says so.

What to print #

One number, at the moment the basis is finished: κ₂ of V, and κ₂ of W if the projection is two-sided. A construction that returns a reduced model without them has thrown away the only evidence that the model means what it claims. It costs an SVD of an n × r matrix on an r that is single digits, against a construction that already spent r sparse solves — which is to say it costs nothing, and it is the difference between a model and a plausible-looking set of small matrices.