The problem the solver was actually given

The polynomial field opened on an object that is not a matrix, and a backward-stable answer to a problem nobody asked measured what turning it into one costs. Both of those substitutions are exact: a linearisation has the eigenvalues of the polynomial, with their multiplicities and their Jordan structure, and the entire loss is in the arithmetic that follows.

A nonlinear eigenvalue problem does not offer that. Take

T(λ) = A − λI + γ√(λ + c) I,

the radiation term a wave problem on a truncated domain produces. It has no linearisation, and the reason is not that nobody has found one: a finite pencil has finitely many eigenvalues and a square root has a branch cut. So every algorithm does the same thing first. It replaces the function by one it can linearise, and only then computes.

T(λ)  →  T̃(λ) = A − λI + r(λ)I  →  a matrix of size n(m+1)  →  eigenvalues
         ^^^^^^^^ approximation    ^^^^^^^^^^^^ exact         ^^^^^ rounding

Three places an error can come from, and the middle one is free. The first happens before any number is rounded. It is a choice, it is computable in advance, and no residual contains it.

The three numbers #

For a computed eigenpair (λ, x) of the approximated problem, at m poles:

The first is what any eigensolver returns, and it is at rounding whatever the approximation was. It does not fall as m rises. Past a point it slowly rises, because each added pole makes the fit’s own basis worse conditioned — so the number a caller reads gets worse while the answer gets four orders better.

The second costs one further evaluation of γ√(λ + c). It falls with the approximation and tracks the forward error to a factor of about twenty at every m.

So the free number says nothing and the nearly-free number says almost everything. Evaluate the residual against the function that was asked about, not against the one handed over is the whole practical content of this essay, and it is one line of code.

The identity, with a term in front #

This site’s spine is

forward error  ⪅  condition number × backward error,

and every field so far has been about separating the two factors on the right. The polynomial field added a third party in six routes to one spectrum: the reformulation, which manufactures a backward error out of nothing when the units are wrong, while the problem’s own conditioning does not move.

This is a fourth, and it sits in front rather than inside:

forward error  ⪅  ‖g − r‖  +  κ × backward error.

The first term is the approximation. It is not amplified by the condition number, it is not produced by the arithmetic, and it is not a property of the problem. It is a property of a decision, made by the caller, before the solve began.

That makes it the cleanest instance of whose fault is it this site has. The problem is as well conditioned as it ever was. The algorithm is backward stable, at 10⁻¹⁵. The answer has four digits. Nobody did anything wrong, and the missing eleven digits were spent by whoever chose the approximant.

The measurement that makes it an identity #

The claim above would be a description if the forward error merely correlated with ‖g − r‖. It does not correlate — it equals it, times one number:

A spread of 1.30 while both fall by two orders. The obvious candidate for the constant is 1/|d/dλ (μ − λ + g(λ))| at the root — how much a perturbation of the function moves the eigenvalue — which does not depend on the approximant and would explain why the column is flat. It is the right mechanism and it is not that number, and the discrepancy is worth a section because the repaired statement is stronger than the one it replaces.

That is the same shape as the ADI identity in why a Gramian can be truncated at all: a bound that is a fixed fraction of its quantity over orders is not a bound, it is an identity with a constant in front, and knowing which of the two is in hand decides whether the number can be used to predict rather than merely to reassure.

A maximum divided by a maximum #

The derivative can be evaluated. At the six roots — 2.939, 3.843, 5.060, 6.339, 7.450, 8.201 — the quantity 1/|f′| with f(λ) = μ − λ + γ√(λ + c) comes out as

1.696 1.573 1.473 1.406 1.366 1.344

and the column above never rises above 0.662. The proposed constant is larger than the measured one at every root, by a factor of between two and three, and it is not even a single number: it varies by 26% across the spectrum while the measured ratio was supposed to be flat.

The reason is in how the ratio was formed. Its numerator is the largest of the six eigenvalue errors; its denominator is the largest of |g − r| over the whole target interval. Those are maxima over two different sets that have no point in common — the peak of the approximation error sits wherever the fit’s equioscillation puts it, and the eigenvalues sit wherever the spectrum of A puts them. No eigenvalue is at the peak, so every one of them sees less error than the maximum, and the quotient comes out below one for a reason that has nothing to do with the derivative.

Taken per eigenvalue instead of as a maximum, the six ratios behind the single value 0.662 at m = 2 are

0.005 0.424 0.599 0.338 0.184 0.662

— a factor of a hundred and thirty between the smallest and the largest. That is the actual spread of the quantity. The 1.30 in the row above is the spread of the maximum of six such numbers across six degrees, and a maximum of a scattered set is a stable statistic while the set is not. This collection has now named that failure four times from four fields, and it is the same shape every time: a sample maximum quoted where a property was meant.

The identity, at the resolution it holds at #

What survives, and survives much better than the original claim, is the pointwise version. Take one eigenvalue, evaluate the approximation error at its own root rather than at the peak, and divide by the derivative there:

|λ − λ*| = |g(λ*) − r(λ*)| ⁄ |f′(λ*)| .

Measured over all six eigenvalues at all six degrees, the two sides agree to between 0.955 and 1.037 of each other at m = 2, and from m = 5 onwards to 1.000 at four digits. Not a bound, not a correlation, not a bound with a constant in front — an identity, holding to the precision that first-order perturbation theory is capable of, and getting tighter as the approximation improves because that is when the second-order term stops mattering.

And the ratio the maximum version reports is stable across the whole range, which is what makes it a usable rule of thumb even though it is not the identity.

The approximant residual falls by six hundred and the ratio does not move. Across reaches of 1, 0.7, 0.5 and 0.4 the approximant residual runs 1.99·10⁻⁸, 9.98·10⁻¹⁰, 7.06·10⁻¹¹ and 3.31·10⁻¹¹ — a factor of 600 — while the forward error divided by the true residual stays between 19.5 and 24.5 at every one of them.

That is the conditioning factor being measured directly and shown not to depend on the approximant. The whole argument of the section is that the forward error is the approximation error carried by a factor belonging to the problem; if the factor were a property of the approximation it would move as the approximation improved six-hundredfold, and it moves by 25%. The range does drift slightly upward — 19.5–21.2 at reach 1 against 21.5–24.6 at 0.4 — which is the second-order term the identity above absorbs, and it is small enough that a caller multiplying a residual by twenty would be right to within a quarter across the range.

The essay’s argument is unchanged and better supported. The forward error is the approximation error, carried to the eigenvalue by a factor that belongs to the problem and not to the approximant. What had to go is the number 0.55, which was never the factor: it was one maximum divided by another, and the identity it was standing in for is exact.

The practical consequence is a sharper version of the advice further down. Computing ‖g − r‖ on the target set gives an upper bound over the whole interval, which is what a caller can do before any eigenvalue exists and is worth doing. Once an eigenvalue does exist, evaluating |g − r| at that λ costs one further evaluation and gives the forward error rather than a bound on it — a certificate rather than a reassurance, from the same two function calls the essay already asks for. assertTheIdentityIsPointwiseAndNotAMaximum holds all four parts: that the quoted ratio lies below 1/|f′| at every root, that the per-eigenvalue ratios behind it span orders, that the maximum of them barely moves anyway, and that the pointwise statement is an identity to four digits.

The exact answer, and where it comes from #

Everything above is measured against a closed form. The nonlinearity here is scalar, so substituting u = √(λ + c) turns the eigen-equation μ − λ + γu = 0 into a scalar quadratic u² − γu − (μ + c) = 0, one per eigenvalue of A. The positive root gives λ = u² − c, and the whole spectrum is written down rather than computed.

Checked by substitution into the original equation at 2.2·10⁻¹⁶.

That standing matters more here than anywhere in this field. Every number above is a distance from a true eigenvalue, and a distance from a better computation would be circular — the better computation would itself be an approximation, and the term this essay is about would be present in both and cancel. It is the habit an answer that is known established for the site and the reason a square root rather than an exponential was chosen for the family.

Why a square root #

Worth one paragraph, because the choice of test problem is doing work.

The delay problem in a problem with infinitely many eigenvalues has an entire nonlinearity, $e^{-\lambda}$, which a polynomial approximates geometrically. Everything in this essay would still be true there and the numbers would be less interesting: the approximation error falls so fast that it stops being the dominant term after a handful of degrees, and the three-way separation collapses.

A branch point does not go away. A polynomial converges to √ algebraically and a rational with poles on the cut converges geometrically, so the choice of approximant is visible, the approximation error is a real term over a useful range of degrees, and there is something for two approximants and one matrix size to compare. The choice of approximant is only visible on a function that has something to be wrong about.

What the free residual is good for #

It is not useless, and saying so precisely matters because “the residual says nothing” is too strong.

‖T̃(λ)x‖ is a correct backward-error statement about the linearised eigenproblem, and it does its job: it certifies that the eigensolver did what it claimed on the matrix it was handed. When it rises — as it does past m = 6 here, to 1.5·10⁻¹² — that is real information, and it is information about the arithmetic, which is the other half of the trade two approximants and one matrix size is about.

What it cannot do is say anything about the distance to a true eigenvalue, because the problem it is a residual for is not the problem that was posed. A caller who reads it as an accuracy estimate has made the mistake this essay exists to name, and the mistake is invisible: the number is small, it is correctly computed, and it means what it says. It just says something else.

That is the same shape as a small residual is not a small error, one level up. There the gap between the two is a condition number, and the condition number can be estimated. Here the gap is an approximation error, and it can be computed exactly — which makes this the easier case, once anybody thinks to look.

How the approximant becomes a matrix #

The middle step is the one this essay calls free, and it is worth showing that it really is, because “exact” is a strong word.

Fit r(λ) = c₀ + Σ αⱼ/(λ − ξⱼ) to g on the target set, with the poles ξⱼ chosen in advance. Then the approximated problem (A + c₀I − λI + Σ αⱼ/(λ − ξⱼ) I)x = 0 becomes a standard eigenvalue problem of size n(m+1) by introducing yⱼ = x/(λ − ξⱼ):

λx  = (A + c₀I)x + Σ αⱼ yⱼ,
λyⱼ = x + ξⱼ yⱼ.

Substitute and the two lines are identities. Nothing is dropped, nothing is truncated, and the eigenvalues of the n(m+1) × n(m+1) matrix are exactly the eigenvalues of T̃ — which is what makes the error decomposition three-term rather than four.

It is not even a pencil, which is worth a sentence given how much of this field is about pencils. The substitution produces a plain matrix because the rational function’s poles are simple and its constant term is finite. A rational approximant with a repeated pole, or one that grows at infinity, gives a pencil instead, and then everything six routes to one spectrum says about which reduction to take applies again on top.

What the step costs is size. Every pole added is n more rows, so a better approximation is a larger matrix and the eigensolve is O(n³m³). That is the trade two approximants and one matrix size prices, and it is the reason the approximation term cannot simply be driven to zero.

What a caller should do #

Three things, in order of cost.

Compute ‖g − r‖ on the target set when the approximant is built. It is a maximum of a scalar function over an interval, it costs a few hundred evaluations, and it is an upper bound on the approximation term before any eigenvalue exists. Nobody reports it and everybody could.

Return the residual against T rather than against T̃. One evaluation of g per computed eigenvalue. It is the number that tracks the forward error, and it is the only one of the three a solver can compute without a closed form.

Do not read the linearisation’s residual as accuracy. It is a check on the eigensolver and it is worth having as one.

The first two together turn an error with three sources into an error with three reported components, which is the arrangement this site has been arguing for since its first field. What is unusual here is how cheap it is: the term nobody reports is the one that is easiest to compute.

The reading a residual invites, and why it is so hard to resist #

There is a reason this mistake is common enough to be worth an essay rather than a footnote, and it is not carelessness.

Every habit this site has built says that a residual is the honest quantity. It is computable without knowing the answer, it is what backward stability is a statement about, and half the collection is an argument for reporting it rather than a forward error nobody can obtain. Then a solver returns one at 10⁻¹⁵ and the habit fires correctly on the wrong object.

What has changed is that there are now two problems, and the residual is a perfectly good answer about the wrong one. That situation has not arisen before in this collection: everywhere else the matrix handed to the solver is the problem, possibly badly scaled or badly linearised, and a residual against it is a residual against the thing asked about up to a transformation that is exact. Here the transformation is not exact, and the word backward stops pointing where it used to.

The repair is not to distrust residuals. It is to be explicit about which problem a residual is against, which is a discipline rather than a technique, and which becomes automatic as soon as somebody writes both numbers down side by side. The measurement above is the argument for writing both down: they differ by ten orders and only one of them moves when the answer improves.

The fourth party, named #

This field has now accumulated a list of things that can be blamed for a wrong answer, and it is worth setting it out because the list is the field’s contribution to the site.

The problem. Its condition number amplifies whatever perturbation reaches it. Nothing can be done about this from inside an algorithm, and measuring it is how a user learns that the answer they want does not exist at the precision they are working in.

The algorithm. Its backward error is the perturbation it manufactures out of rounding. A good algorithm keeps this at the unit roundoff and a bad one does not, and the two are distinguishable by measurement rather than by inspection.

The reformulation. New with this field. A linearisation is exact in the algebra and changes the conditioning of the object being factorised, so a badly scaled quadratic loses ten orders through a step that added no error of its own. The repair is two lines and the diagnosis is that the quadratic’s condition number does not move while the linearised matrix’s runs to 10¹¹.

The approximation. New here. It is committed before the solve, it is not amplified by anything, it appears in no residual against the object the solver factorised, and it is the only one of the four that is a decision. The other three are properties — of a problem, of an algorithm, of a reformulation — and this one is a number somebody typed.

That last point is what makes it worth separating rather than folding into the reformulation. A reformulation is chosen once, by whoever wrote the solver, and a user of the solver inherits it. An approximation is chosen every time, by the user, and usually without being told that a choice is being made at all.

One sentence #

A solver that reports a residual of 10⁻¹⁵ against a problem it constructed, on an answer that is wrong in the fourth digit, has not made a mistake and has not told anybody anything. The number that would have told them costs one evaluation of the function they asked about, and nobody computes it.