A backward-stable answer to a problem nobody asked

The spine of this site is one inequality:

forward error ⪅ condition number × backward error.

A good algorithm returns the exact answer to a nearby problem, the condition number says how much the problem magnifies a perturbation, and only their product is what anybody sees. Every essay here that reports an error reports which of the two factors it came from, and the exact answer to a nearby problem is where the habit was set.

There is a phrase in that sentence doing more work than it looks: a nearby problem. Nearby in what, and a problem of which kind? For a linear solve those questions have obvious answers and nobody has to ask them. For a quadratic eigenvalue problem they do not, and this essay is the measurement of what happens when the obvious answer is the wrong one.

The family is one problem in different units #

The sweep is built out of a substitution rather than a perturbation, and that is what makes it airtight. Replace λ by γμ in Q(λ) = λ²M + λC + K and the coefficients become

(γ²M, γC, K),

with the spectrum divided by γ exactly. Nothing is approximated, nothing is chosen, and the map back is a multiplication. If time is measured in milliseconds instead of seconds, γ is a thousand; if a stiffness is given in newtons per metre and a mass in tonnes, γ is whatever the unit conversion says. Every stop of the sweep is the same physical problem written down in a different unit system, and the closed form of the damped chain still gives its exact spectrum at every one of them.

So there is no wriggle room about what the right answer is, at any stop.

Two backward errors, and they are not the same number #

The linearised matrix A is 2n × 2n, and the real Schur factorisation applied to it is backward stable in the sense the form a real matrix can reach establishes: the computed eigenpair is exact for A + ΔA with ‖ΔA‖/‖A‖ at the rounding level. That theorem is true here and nothing in this essay disputes it.

The backward error of the same pair for the quadratic is a different quantity:

η(λ, x) = ‖Q(λ)x‖ / ((|λ|²‖M‖ + |λ|‖C‖ + ‖K‖)‖x‖),

the smallest relative perturbation of M, C and K severally for which the pair is exact. The denominator is the whole content. A residual divided by ‖Q(λ)‖ would be the backward error of a perturbation of Q as a single object — which is not a quadratic problem, and so not an answer to anything anybody asked.

And the two numbers can be far apart, for a reason that is visible in the linearisation’s own shape. C₁ has identity blocks in it. Their norm is one whatever the coefficients are, so a perturbation that is small relative to ‖A‖ can be enormous relative to the smallest of ‖M‖, ‖C‖ and ‖K‖ — and there is no reason for the perturbed matrix to be a linearisation of any quadratic at all.

What is measured #

At n = 8, on an overdamped chain, with γ running from 1 to 10⁸:

The middle column moves by a factor of about a thousand across the sweep. The outer two move by eleven orders of magnitude.

The solver did nothing wrong at any stop. It was handed a matrix and it returned eigenpairs that are exact for a matrix a rounding away. What it was not handed is the problem, and the residual it prints — the middle column — is a statement about the object it was given rather than about the object the caller has.

That the forward error tracks the quadratic’s backward error and not the linearisation’s is the half that makes this more than a definition. Measured, the ratio of forward error to η(quadratic) stays inside a factor of 5.8 across the whole sweep, while the ratio to η(linearisation) moves by a factor of 1.5·10⁷. One of those two numbers predicts the answer’s accuracy and the other one does not, and the one that does not is the one a solver reports.

The condition number that does not move #

Take the spine’s inequality seriously and there should be a condition number in it. There is, and it separates in the same way.

The eigenvalue condition number of the quadratic at a simple λ is

κ(λ) = (|λ|²‖M‖ + |λ|‖C‖ + ‖K‖)‖x‖‖y‖ / (|λ| |yᵀQ′(λ)x|),

with x and y the right and left null vectors of Q(λ). It is the direct analogue of a condition number for one eigenvalue — the 1/|yᵀx| that governs how far a single eigenvalue of a matrix moves — with the two changes that make it a statement about a polynomial: the numerator measures the three coefficients separately, and the denominator carries Q′(λ) rather than the identity.

Measured: κ(quadratic) is 4.98 at every stop, over six decades of γ — a spread of 1.0002. It cannot move, because a change of units is not a change of problem and the quantity is built to be invariant under one. κ(linearisation) goes from 49.66 to 7.8·10¹¹.

So the accounting closes. The problem is as well conditioned at γ = 10⁶ as it was at γ = 1; the answer has lost six digits; and the entire loss is attributable to a substitution the caller made before the solver was reached. That is whose fault is it with a third party in it — not the algorithm and not the problem, but the reformulation — and it is the first time this site has had to name one.

What a perturbation of the linearisation would have to be #

It is worth being concrete about why the two backward errors can differ at all, because the argument is short and it explains the size of the gap rather than merely its existence.

Suppose a solver returns a pair that is exact for A + ΔA with ‖ΔA‖ = u‖A‖. For that to say anything about the quadratic, ΔA would have to be a perturbation that keeps A a linearisation of some quadratic — that is, it would have to leave the identity blocks alone and disturb only the blocks built from M, C and K. Nothing in the error analysis of a Schur factorisation promises that. The perturbation it delivers is an arbitrary 2n × 2n matrix of the right size, and an arbitrary matrix of that size is not the linearisation of anything. That is the same shape of objection a nearby problem of the wrong kind makes about a structured problem, and it bites harder here, because the structure being violated is not a modelling assumption but the definition of the object.

So the question becomes: given the ΔA the solver actually produced, what is the smallest perturbation of the three coefficients that reproduces the same eigenpair? That is what η for the quadratic computes, by construction, and the answer involves dividing the residual by |λ|²‖M‖ + |λ|‖C‖ + ‖K‖ rather than by ‖A‖. When the three norms are wildly different those two denominators are wildly different, and the ratio between them is what the sweep is walking along.

The bound the literature gives for C₁ has the same shape: the quadratic’s backward error is at most the linearisation’s times a factor built from the coefficient norms and |λ|, and the factor is one when everything is of comparable size. Measured here, the ratio grows by about one order per decade of γ, so the bound is not merely satisfied — it is nearly attained, which is the only thing that makes a bound worth quoting.

What a caller can actually see #

The uncomfortable part of this measurement is that every quantity a caller normally has is flat. The residual for the linearised matrix is at the rounding level. The computed eigenvalues look plausible: they are negative reals in an overdamped problem, spread over the range the physics suggests, in conjugate pairs where they should be. The iteration converged. Nothing warns.

Two things are visible without the closed form, and both are worth having.

The coefficient norms. ‖M‖, ‖C‖ and ‖K‖ are three cheap numbers, and their spread is what predicts the trouble — the same diagnostic role the units the matrix is measured in gives to a row scaling, and available for the same reason: it is a property of the data rather than of the run. At γ = 10⁸ they are 2.8·10¹⁶, 2.0·10⁹ and 6.8 — a spread of 4·10¹⁵ — and nothing else in the computation is that far out of balance.

The quadratic’s own backward error. η(λ, x) costs one matrix-vector product with each coefficient and three norms. It is the number that should be reported and, on the evidence of the table above, the only one worth reporting — and computing it is the same act orthogonal is a number recommends for a factorisation: measure the property being claimed rather than the one the routine happens to return. A library that prints the linearised residual is answering a question about its own internals; a library that prints η is answering the caller’s.

Which half of the eigenvector, and why it is worth saying #

C₁’s eigenvector is [λx; x]. So the quadratic’s eigenvector can be read off either half of the computed 2n-vector, and the two halves differ by a factor of λ. They therefore have different backward errors, and on a badly scaled problem the difference is not small.

The measurement above reports the better of the two at every stop, deliberately: a comparison that took the worse half would be measuring a choice this essay is not arguing about, and the finding survives the most favourable reading of the method. It is worth knowing that the choice exists, because a library that takes the top half unconditionally is leaving a factor of |λ| on the table, and |λ| is 10⁻⁸ at the far end of this sweep.

What this does not say #

Three things, because the finding is easy to over-read.

It is not that linearising is wrong. There is nothing else to do. A matrix polynomial has no QR algorithm of its own, and the linearisation is what makes fifty years of eigenvalue machinery applicable to a problem it was not written for. The finding is that the substitution has a price, that the price is measurable, and that it is invisible in every number a solver returns.

It is not that the sweep is exotic. γ = 10⁸ is a unit conversion. The measured examples that motivated this analysis in the literature are structural models where mass is in kilograms and stiffness in newtons per metre, which is exactly a mismatch of this size, and nobody involved thought they were doing anything unusual.

And it is not unfixable. The repair is two lines of scaling computed from three norms, it recovers all ten orders, and it is the scaling that buys ten orders. What this essay establishes is why anyone should bother: without the measurement, a code that scales looks like a code with a superstition in it, since every residual it prints was already at the rounding level.

Where the eigenvector comes from, and what it costs to get one #

A detail that matters for the measurement’s honesty. The real Schur factorisation returns eigenvalues; the eigenvectors here come from inverse iteration on the linearised matrix, four steps from a fixed deterministic start with the shift moved off the eigenvalue by a relative 10⁻¹⁰.

That is the standard device and it looks alarming the first time: the matrix being solved with is deliberately nearly singular. It works for exactly that reason — the ill conditioning aims the solve at the eigenvector, since the component along it is the one amplified — and it is the same argument the inverse that is never formed makes about why a solve is not an inverse. What matters here is that the resulting pairs are good ones for the linearisation: their η(linearisation) is at the rounding level at every stop of the sweep, which is the column the table’s middle holds. If they were not, this essay would be measuring a bad eigenvector rather than a bad substitution.

The sweep in the other direction #

The table above raises γ, which makes ‖M‖ large. Lowering it makes ‖K‖ large instead, and the same thing happens with the roles exchanged: the eigenvalues become enormous, |λ|²‖M‖ dominates the scale in the denominator, and the small-modulus eigenvalues are the ones that lose their digits. There is nothing special about the direction, which is worth saying because it kills the obvious defence — that the problem is large numbers and could be avoided by working in sensible SI. What matters is the ratio of the three coefficient norms, and no choice of units makes all three of ‖M‖, ‖C‖ and ‖K‖ comparable for every problem, because they carry different physical dimensions and their ratios have units of time.

That is why the repair is a change of variable rather than a convention: the right γ is a property of the problem’s own numbers and has to be computed from them.

The thing to carry out of this #

A backward error is a promise about which problem the answer is exact for, and the promise is only as good as the reader’s assumption about what the data is. This site has met that once before, from a different direction: a nearby problem of the wrong kind shows a backward-stable computation whose nearby problem is not of the same structure as the one posed, so a perturbation the analysis permits is one the model forbids.

Here the mismatch is more elementary and more common. The nearby problem is not a quadratic eigenvalue problem at all. It is a nearby 2n × 2n matrix, and the set of quadratics that linearise to it is, in general, empty.

At other settings #