Numerical Linear Algebra
where the algebra stops being the arithmetic
The eigenvalue problem that is not linear

The problem the solver was actually given

A linearisation is exact — it has the polynomial's eigenvalues, with their multiplicities, and the whole loss is arithmetic. A nonlinear eigenvalue problem does not offer that. Every algorithm replaces the function first, and the term that replacement contributes is committed before any number is rounded and appears in no residual.
36 min read 31 figures Whose fault is itWhere the algebra stops being the arithmeticMeasured, not asserted

Worth reading first: A problem with infinitely many eigenvalues · A matrix that depends on its own eigenvalue · The exact answer to a nearby problem · Elimination is a sequence of choices.

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.

Three residuals, and the one a solver returns is the one about nothingFor a rational approximant with m poles, on a target set reaching to 0.3: ‖T̃(λ)x‖ — the residual against the problem the eigensolver was handed — stays at 1.2·10⁻¹⁵ to 1.5·10⁻¹² — it does not fall with m, and past a point it slowly rises, because each added pole makes the fit's own basis worse conditioned; ‖T(λ)x‖, which costs one further evaluation of γ√(λ + c), falls with the approximation; and the forward error against the closed form falls with it, staying 20.5 to 24.8 times larger. The free number rises by 1306 across the sweep while the answer improves by 885. So the free number says nothing and the nearly-free number says almost everything. Evaluate the residual against the function you asked about, not against the one you handed over is the whole practical content of this field.23456710⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²poles in the approximantresiduals and errorforward error‖T(λ)x‖‖T̃(λ)x‖one extra evaluationagainst the approximant1.5·10⁻¹²against the problem asked1.9·10⁻⁶forward error4.8·10⁻⁵‖g − r‖ there8.5·10⁻⁵the free residual is flatand the answer is not
Fig. 1 Three residuals for the same computed eigenvalues. The flat one is the one a solver returns.

The three numbers

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

quantity at m = 2 at m = 7 what it is
‖T̃(λ)x‖ 1.2·10⁻¹⁵ 1.5·10⁻¹² the residual against the problem handed over
‖T(λ)x‖ 1.8·10⁻³ 1.9·10⁻⁶ the residual against the problem asked
|λ − λ*| 4.2·10⁻² 4.8·10⁻⁵ the forward error

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:

poles m 2 3 4 5 6 7
forward ÷ ‖g − r‖ 0.662 0.582 0.508 0.598 0.554 0.562

A spread of 1.30 while both fall by two orders. The constant is 1/|d/dλ (μ − λ + g(λ))| at the root — how much a perturbation of the function moves the eigenvalue — and it does not depend on the approximant, which is why it is flat.

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.

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^{−λ}, 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.

At other settings

Three residuals, and the one a solver returns is the one about nothingFor a rational approximant with m poles, on a target set reaching to 1: ‖T̃(λ)x‖ — the residual against the problem the eigensolver was handed — stays at 1.8·10⁻¹⁵ to 2·10⁻⁸ — it does not fall with m, and past a point it slowly rises, because each added pole makes the fit's own basis worse conditioned; ‖T(λ)x‖, which costs one further evaluation of γ√(λ + c), falls with the approximation; and the forward error against the closed form falls with it, staying 19.5 to 21.2 times larger. The free number rises by 1.14·10⁷ across the sweep while the answer improves by 2.79·10⁴. So the free number says nothing and the nearly-free number says almost everything. Evaluate the residual against the function you asked about, not against the one you handed over is the whole practical content of this field.23456710⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²poles in the approximantresiduals and errorforward error‖T(λ)x‖‖T̃(λ)x‖one extra evaluationagainst the approximant2·10⁻⁸against the problem asked2.5·10⁻⁸forward error5.4·10⁻⁷‖g − r‖ there3.4·10⁻⁷the free residual is flatand the answer is not
Fig. 2 A target set that stands clear of the branch point, where the approximation term is smallest.
Three residuals, and the one a solver returns is the one about nothingFor a rational approximant with m poles, on a target set reaching to 0.7: ‖T̃(λ)x‖ — the residual against the problem the eigensolver was handed — stays at 1.1·10⁻¹⁵ to 10·10⁻¹⁰ — it does not fall with m, and past a point it slowly rises, because each added pole makes the fit's own basis worse conditioned; ‖T(λ)x‖, which costs one further evaluation of γ√(λ + c), falls with the approximation; and the forward error against the closed form falls with it, staying 19.6 to 24.8 times larger. The free number rises by 8.72·10⁵ across the sweep while the answer improves by 8634. So the free number says nothing and the nearly-free number says almost everything. Evaluate the residual against the function you asked about, not against the one you handed over is the whole practical content of this field.23456710⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²poles in the approximantresiduals and errorforward error‖T(λ)x‖‖T̃(λ)x‖one extra evaluationagainst the approximant10·10⁻¹⁰against the problem asked6.8·10⁻⁸forward error1.7·10⁻⁶‖g − r‖ there3.1·10⁻⁶the free residual is flatand the answer is not
Fig. 3 Reaching a little further, and every curve but the flat one moves.
Three residuals, and the one a solver returns is the one about nothingFor a rational approximant with m poles, on a target set reaching to 0.5: ‖T̃(λ)x‖ — the residual against the problem the eigensolver was handed — stays at 1.2·10⁻¹⁵ to 7.1·10⁻¹¹ — it does not fall with m, and past a point it slowly rises, because each added pole makes the fit's own basis worse conditioned; ‖T(λ)x‖, which costs one further evaluation of γ√(λ + c), falls with the approximation; and the forward error against the closed form falls with it, staying 21.2 to 24.8 times larger. The free number rises by 6.05·10⁴ across the sweep while the answer improves by 2396. So the free number says nothing and the nearly-free number says almost everything. Evaluate the residual against the function you asked about, not against the one you handed over is the whole practical content of this field.23456710⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²poles in the approximantresiduals and errorforward error‖T(λ)x‖‖T̃(λ)x‖one extra evaluationagainst the approximant7.1·10⁻¹¹against the problem asked3.7·10⁻⁷forward error9.1·10⁻⁶‖g − r‖ there1.5·10⁻⁵the free residual is flatand the answer is not
Fig. 4 Halfway, where the forward error at seven poles is three orders worse than at reach one.
Three residuals, and the one a solver returns is the one about nothingFor a rational approximant with m poles, on a target set reaching to 0.4: ‖T̃(λ)x‖ — the residual against the problem the eigensolver was handed — stays at 2.4·10⁻¹⁵ to 3.3·10⁻¹¹ — it does not fall with m, and past a point it slowly rises, because each added pole makes the fit's own basis worse conditioned; ‖T(λ)x‖, which costs one further evaluation of γ√(λ + c), falls with the approximation; and the forward error against the closed form falls with it, staying 21.5 to 24.8 times larger. The free number rises by 1.4·10⁴ across the sweep while the answer improves by 1456. So the free number says nothing and the nearly-free number says almost everything. Evaluate the residual against the function you asked about, not against the one you handed over is the whole practical content of this field.23456710⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²poles in the approximantresiduals and errorforward error‖T(λ)x‖‖T̃(λ)x‖one extra evaluationagainst the approximant3.3·10⁻¹¹against the problem asked8.4·10⁻⁷forward error2.1·10⁻⁵‖g − r‖ there3.4·10⁻⁵the free residual is flatand the answer is not
Fig. 5 And further still.
The error that happens before any arithmetic: |g − r| on the target setγ√(λ + c) minus its approximant, plotted against the distance from the branch point at −0.4. The target set runs from 0.602 to 8.99, and its left end is 1 from the branch point — a choice made before anything is computed. The rational approximant with 8 poles clustered on the cut reaches 2.22·10⁻⁵; the polynomial of degree 9, which linearises to a matrix of the same size, reaches 1.71·10⁻⁴. Both errors are committed at this stage. No residual computed later contains either of them, because a solver's residual is taken against whichever of these two functions it was handed.110¹10⁻⁸10⁻⁶10⁻⁴distance from the branch point, λ + c|g − r|the eigenvaluesdegree 98 poleschosen, not computedreach0.3left end, from the cut1rational, worst2.2·10⁻⁵polynomial, worst1.7·10⁻⁴linearisation size, both54committed before the solveand invisible to it
Fig. 6 The term itself, drawn on the interval where it was committed.
The error that happens before any arithmetic: |g − r| on the target setγ√(λ + c) minus its approximant, plotted against the distance from the branch point at −0.4. The target set runs from 2.94 to 8.99, and its left end is 3.34 from the branch point — a choice made before anything is computed. The rational approximant with 8 poles clustered on the cut reaches 4.3·10⁻⁸; the polynomial of degree 9, which linearises to a matrix of the same size, reaches 1.91·10⁻⁷. Both errors are committed at this stage. No residual computed later contains either of them, because a solver's residual is taken against whichever of these two functions it was handed.10¹10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶distance from the branch point, λ + c|g − r|the eigenvaluesdegree 98 poleschosen, not computedreach1left end, from the cut3.3rational, worst4.3·10⁻⁸polynomial, worst1.9·10⁻⁷linearisation size, both54committed before the solveand invisible to it
Fig. 7 And on an easy target set, where it is fourteen orders smaller.
36 eigenvalues come back, 6 of them mean somethingThe real eigenvalues of the n(m+1) = 36 that a rational approximant with 5 poles produces, on the real line, with the branch point at −0.4 marked. The 6 discs to the right are the answers, agreeing with the closed form to 7.74·10⁻⁴. The 6 crosses to the left are exact eigenvalues of the approximant lying where γ√(λ + c) is not a real number at all — so the residual that says nothing about accuracy is decisive as a filter: against T̃ all 12 look perfect, and against T the spurious ones do not have a residual. The remaining 24 are complex and are discarded the same way.-10123456789-101λbranch point at −0.4the closed formwhich of these is an answereigenvalues returned36wanted6past the branch point6complex24worst against the closed form7.7·10⁻⁴all of them exactfor a problem nobody asked
Fig. 8 What else the linearisation returns.
Two approximants at equal linearisation size, at a reach of 0.06‖g − r‖ against the size of the matrix the eigensolver is then handed: n(m+1) for a rational approximant with m poles and nd for a polynomial of degree d, so d = m + 1 is the fair comparison. At this target set the ratio of the polynomial's error to the rational's runs 1.2, 2.9, 8.0, 23.7, 73.2 across the five sizes. On a target set that stands clear of the branch point the two are the same measurement and the ordering flips with the noise; on one that reaches it the rational pulls away, because a polynomial converges algebraically against an algebraic singularity. Neither method is better — the question is different, and which question it is was decided by whoever said where to look.263850627410⁻⁵10⁻³10⁻¹rows in the linearisation‖g − r‖ on the target setpolynomialrationalequal cost, two basesreach0.06largest size drawn78rational there2·10⁻⁵polynomial there0.0015the ratio73no winner on an easy targetand two orders on a hard one
Fig. 9 Two approximants at equal linearisation size.
One quadratic eigenvalue problem in 9 systems of units: what the solver reports and what the answer is worthAn overdamped chain of 6 masses, with λ replaced by γμ so that the coefficients become (γ²M, γC, K). That substitution is exact in both directions and divides the spectrum by γ exactly, so the closed form is still available and every error here is measured against it. The backward error of the eigenpair for the LINEARISED MATRIX — the residual a solver's own error analysis is about — is 7.1·10⁻¹⁶ at γ = 1 and 2.57·10⁻¹² at γ = 108 — it moves by a factor of 7022 while the other two move by 7.37·10¹¹. The backward error for the QUADRATIC, which is what the person who posed the problem is entitled to, grows by 7.03·10¹¹ across the same sweep, and the forward error follows it: 9.14·10⁻¹⁵ to 0.00674. Nothing went wrong with the solver at any stop.0246810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ γ, the change of unitsrelative errorforward errorη, the quadraticη, the linearisationagainst a closed formη(linearisation), worst2.6·10⁻¹²η(quadratic), worst4.6·10⁻⁴forward error, worst0.0067coefficient spread4.2·10¹⁵the solver is right at every stopabout a problem nobody asked
Fig. 10 The backward error of a linearisation, which is the term this essay puts something in front of.
The condition number the problem has, and the one the solver's error analysis is written againstThe same eigenvalue of the same overdamped chain of 10 masses, in seven systems of units. Its condition number as an eigenvalue of the QUADRATIC — Tisseur's, with the three coefficient norms in the numerator and yᵀQ′(λ)x in the denominator — is 5.554 at γ = 1 and 5.554 at γ = 10⁶, a spread of 1 over six decades: it cannot move, because a change of units is not a change of problem. Its condition number as an eigenvalue of the LINEARISED MATRIX runs 55.22 to 8.685·10¹¹, a factor of 1.57·10¹⁰. The forward error follows the second one, and the first one is the honest description of the problem — so the substitution has manufactured an ill conditioning that belongs to the algorithm rather than to the question.012345610⁻¹10²10⁵10⁸10¹¹log₁₀ γ, the change of unitscondition numberthe linearisationthe quadraticone problem, two amplifiersκ(quadratic), first5.6κ(quadratic), last5.6κ(linearisation), last8.7·10¹¹how far the first moved1the problem is as well conditioned as everand the method is not
Fig. 11 And the condition number that multiplies it.
Six ways of computing one spectrum, at a change of units of 106Three linearisations — the first companion form, the second, and the symmetric member of the DL(Q) family — each reduced to a standard eigenvalue problem in two ways: by inverting the leading coefficient, whose eigenvalues are λ, and by inverting the trailing one, whose eigenvalues are 1/λ. All six have exactly the eigenvalues of the quadratic in exact arithmetic. Measured against the closed form at γ = 106, the best is first/leading at 5.473·10⁻⁶ and the worst second/trailing at 2.234·10⁻⁴, a spread of 40.81. Under the leading reduction the symmetric form and the first companion form are THE SAME MATRIX — −A₁⁻¹A₀ is [[−M⁻¹C, −M⁻¹K], [I, 0]] for both — so the distinction between them is one the inversion discards, and it survives only under the trailing one.first · leading5.47·10⁻⁶first · trailing4.62·10⁻⁵second · leading4.27·10⁻⁵second · trailing2.23·10⁻⁴symmetric · leading5.47·10⁻⁶symmetric · trailing4.05·10⁻⁵all six are the same algebrabest route5.5·10⁻⁶worst route2.2·10⁻⁴spread across the six41condition of the linearisation8.3·10¹²the spectra agreeand the arithmetic does not
Fig. 12 Six routes to one spectrum, where the third party was the reformulation.
What two lines of scaling are worth: the same quadratic in nine systems of unitsThe forward error against the closed form for an overdamped chain of 8 masses, before and after Fan–Lin–Van Dooren scaling. The unscaled curve runs 7.49·10⁻¹⁴ to 0.001258 — every digit gone by the far end — and the scaled one runs 1.26·10⁻¹³ to 8.23·10⁻¹⁴, flat to within a factor of 2.01. The scaling is γ = √(‖K‖/‖M‖) and δ = 2/(‖K‖ + γ‖C‖), computed from three norms and nothing else, and the map back is λ = γμ with no rounding in the statement. Flatness is the half that matters: after scaling every stop of the sweep IS the same problem, so there is nothing left for the change of units to do.0246810⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ γ, the change of unitsforward erroras writtenafter scalingone change of variableunscaled, worst0.0013scaled, worst1.7·10⁻¹³orders recovered10scaled coefficient spread4.5the answer was never the problemthe units were
Fig. 13 And the two lines that removed it.
Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u30×30, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 14 A small residual and a large error, in the field that named the distinction.
Backward error, forward error and the condition number, with measured valuesTwo boxes at the top — the problem posed and the nearby problem the algorithm answered exactly — and two answers below them, with the distances between all four labelled by numbers from a Hilbert solve.the problem you posedA = H10b = A·(1, 2, …, 10)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 2.3·10⁻¹⁷the answer you wantedx = (1, 2, …, 10), exactlythe answer you gotx̂, wrong by 2.7·10⁻⁴ relativebackward error 2.3·10⁻¹⁷forward error 2.7·10⁻⁴κ = 1.6·10¹³κ · η = 3.6·10⁻⁴, and the measured forward error is 2.7·10⁻⁴.The algorithm is not at fault. The problem is.H10, LU with partial pivotingresidual and error differ
Fig. 15 The identity this essay adds a term to.
How much a perturbation of the right-hand side is amplified, κ = 10⁴The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10000worst found 76206×6, 200 directionsmedian reaches 0.29 of κ
Fig. 16 The amplifier, which does not touch the new term.
The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 17 The standing every number here is measured against.
A spectrum with no last eigenvalue: 12 of the infinitely many solutions of a delay problemT(λ) = A − λI + 0.5e^−λI for the model matrix at n = 4. The problem decouples in A's eigenvector basis, so each of its 4 eigenvalues μ contributes the scalar equation μ − λ + γe^−λ = 0, whose solutions are λ = μ + W_k(γe^−μ) — one for every branch of the Lambert W function, and therefore countably infinitely many, running off to the left along two curves. The marks are the closed form, each checked by substitution to 1.2·10⁻¹⁵. There is no computation that returns this set: the only well-posed question is which eigenvalues lie inside a region, and the circle of radius 4 contains 4 of them.-3-2-101234-5-3-1135real partimaginary part4 insidea countable spectruminside the contour4drawn12existing∞worst branch residual1.2·10⁻¹⁵there is no last eigenvalueso the question has to change
Fig. 18 The other nonlinear problem in this field, whose function has no branch point.
Counting the eigenvalues inside a circle: an integral that is an integer once it has convergedThe number of eigenvalues of the delay problem inside |z| = 1.5, computed as (1/2πi)∮tr(T(z)⁻¹T′(z))dz by the trapezoidal rule. The true count is 1, from the closed form. The error runs 0.473, 0.39, 0.289, 0.142, 0.0268, 7.58·10⁻⁴ at 4 to 128 points — a straight line on this axis is a constant factor a doubling, and this falls faster than that, which is what exponential convergence looks like when the integrand is analytic on the contour. The answer is then ROUNDED, and 0 of the 6 quadratures round to the wrong integer while looking exactly as confident as the ones that do not. The only evidence available is the distance from the nearest integer, which is the quantity plotted.11.31.61.92.210⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10²log₁₀ quadrature pointsdistance from the true counthalf an eigenvaluean integer, eventuallytrue count1at 4 points1at 128 points1finest error7.6·10⁻⁴the integral is an integerand a rounding hides how far it was
Fig. 19 And the only finite question that can be asked of it.
The 16 eigenvalues of an 8 × 8 quadratic eigenvalue problem, computed and in closed formλ²M + λC + K for a chain of 8 masses with C = 0.3M + 0.1K. The crosses are the closed form — one scalar quadratic per eigenvalue of K, whose roots are known exactly — and the discs are the eigenvalues a real Schur factorisation returns from the 16 × 16 first companion linearisation. There are 16 of them for a matrix with 8 rows, of which 16 are complex and arrive in conjugate pairs, so the eigenvectors cannot be independent: 9 vectors in 8 dimensions never are. The worst disagreement between the two routes is 2.12·10⁻¹⁵, and the routes share nothing but the three coefficient matrices.00.0939569-2-1.24498-0.4899610.2650581.020081.7751real partimaginary parttwo routes, one spectrumeigenvalues16rows8complex16against the closed form2.1·10⁻¹⁵n rows and 2n eigenvaluesso the eigenvectors are not a basis
Fig. 20 The field’s opening picture, where the substitution was exact.
Where the missing eigenvalues went: an 6 × 6 quadratic with 1 of its masses removedA chain of 6 masses with the first 1 of them set to zero, so M is singular of rank 5. det Q(λ) is a polynomial of degree at most 12; interpolated exactly in BigInt rationals at 13 nodes it has degree 11, so 1 of the 12 eigenvalues are at infinity — the same object a descriptor pencil has, arriving here because a degree of freedom with no inertia is an algebraic constraint. The float route counts the singular values of M judged to be zero and reaches 1, backed by a gap of ∞ between consecutive singular values. One integer, two routes, and only the second of them is a decision.finite eigenvalues (degree of det Q)11at infinity (2n − degree)1at infinity, by the rank of M12n, if M were nonsingular12a degree, not a decisiondegree of det Q11at infinity1by the rank of M1singular-value gap∞the count is a degreeand the other route is a judgement
Fig. 21 And the count that made it so.
The roots of ∏(x − i) up to degree 10, computed from the coefficients that define themA polynomial whose roots are the integers 1 to n, expanded exactly in BigInt and then handed to the route every library takes: the eigenvalues of the companion matrix. The measured relative error runs 1.711·10⁻¹³ at degree 6 to 1.127·10⁻¹⁰ at degree 10 — the third digit of an answer whose truth is an integer. The other curve is the root condition number times the unit roundoff, which is a prediction made from the coefficients alone, and it tracks the measurement to within two orders at every degree. Every coefficient here is still exactly representable, which is what makes the failure a conditioning one rather than a storage one.5791110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²degreerelative error of the worst rootmeasuredκ × ua polynomial is its coefficientsworst root, measured1.1·10⁻¹⁰predicted, κ × u1.4·10⁻¹⁰root condition number6.1·10⁵largest coefficient1.3·10⁷the roots are integersand the coefficients are not the roots
Fig. 22 Where a polynomial’s coefficients stop determining its roots.
Fitting the same degree-10 polynomial in two basesOn the left, the data and two fitted curves that lie on top of each other. On the right, the condition numbers of the two design matrices, ten orders of magnitude apart.00.10.20.30.40.50.60.70.80.910.20.40.60.81xyboth fits, drawn on top of one anothercondition number of the design matrixmonomial2.1·10⁷Chebyshev2.3largest fitted coefficientmonomial176Chebyshev0.51rms residual: 4.2·10⁻⁶ and 3.1·10⁻⁷ — the data is fitted either way.30 points, degree 10, single precisionthe basis is part of the problem
Fig. 23 And the basis those coefficients are written in.
Two routes to e^A as a 8×8 bidiagonal matrix walks up to a defective oneA_δ is upper bidiagonal with diagonal λ, λ+δ, …, so its eigenvalues are exact and its eigenvectors have a closed form. As δ falls the eigendecomposition route's relative error rises to 2.92·10⁶⁵ while scaling and squaring's falls to 4.09·10⁻¹², which is the distance ‖A_δ − A₀‖ = 7·10⁻¹² to the matrix whose exponential is known in closed form. The two agree at the right of the axis, where the matrix really is far from the limit; the divergence is κ(V), which reaches 3.3·10⁸².10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹⁴10⁻¹10¹²10²⁵10³⁸10⁵¹10⁶⁴δ, the gap between consecutive eigenvaluesrelative error in e^Aan answer with no correct digitsV f(Λ) V⁻¹scaling and squaring‖A_δ − A₀‖exact eigenvalues throughoutκ(V) at the smallest δ3.3·10⁸²eigen route2.9·10⁶⁵scaling and squaring4.1·10⁻¹²distance to the limit7·10⁻¹²the eigenvalues are the diagonaland they are exact at every stop
Fig. 24 Approximating a function of a matrix, in the spectral field.
‖Aᵏ/k!‖ against ‖e^A‖ for Moler and Van Loan's 2×2, scaled by 0.25The terms of the Taylor series rise to 81.6 at k = 4 and fall back; the sum they produce has norm 5.37, drawn as the flat line. The ratio is 15.2, so that many digits are cancelled away before truncation is a question, and the computed exponential has a relative error of 1.87·10⁻¹⁵ after 36 terms. Scaling and squaring on the same matrix returns 1.65·10⁻¹⁴.0510152025303510⁻⁶10⁻⁴10⁻²110²k‖Aᵏ/k!‖‖e^A‖ = 5.4largest term 82what the series throws awaylargest term82‖e^A‖5.4digits cancelled away15error after the sum1.9·10⁻¹⁵every term is computed correctlyand the sum has lost seven digits
Fig. 25 And what a truncated series costs there.
Where a quadratic stops being hyperbolic, located by a Cholesky and by a sineThe overdamping margin min over modes of (βκ)² − 4κ, for a chain of 10 masses, against β. It reaches zero at β* = 1/sin(π/2(n+1)) = 7.02667418333, which is the closed form. Bisecting on a completely different question — does a Cholesky of −Q(μ) complete for some μ — gives 7.02667418333, agreeing to 14 digits. Neither route computes an eigenvalue. The marks below the axis are the largest imaginary part in the computed spectrum, which is zero to the rounding level above β* and not below it, so a third route agrees with the other two about where the boundary is.56.171117.342228.513349.6844510.855610⁻¹110¹10²stiffness damping βoverdamping marginβ* = 7.02667two routes to a boundaryclosed form β*7by certificate7difference1.9·10⁻¹³bisection steps44a factorisation that completesand a sine, agreeing to twelve digits
Fig. 26 A certificate anybody can re-check, which is what a residual against T is.
A double root, approached: the pair separates like √ε and the accuracy fails like √uA chain of 12 masses at β = β*(1 + ε), where β* is the critical damping and the smallest mode's two eigenvalues coincide at ε = 0. The separation of that pair is 0.6909√ε at every one of twelve decades — a spread of 1.013 in the constant — and the computed spectrum's worst relative error against the closed form rises as the pair closes, reaching 7.465·10⁻⁸ at the boundary itself. That is √u times a small constant: half the digits, on a problem where κ(K) is 67.83 at every stop and the coefficients are integers. Nothing here is ill conditioned in any sense this site has used before; what is unbounded is the derivative of the map from a coefficient to a double root.-15-13-11-9-7-5-3-110⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ ε, distance past the critical dampingseparation, and relative errorat ε = 0: 7.47·10⁻⁸separationcomputed errornothing is ill conditionedthe √ε constant0.69spread of it, twelve decades1error at the boundary7.5·10⁻⁸κ(K), unchanged throughout68half the digitsand no condition number to blame
Fig. 27 Where half the digits go for a reason that is nobody’s fault.
The residual conjugate gradients reports and the residual of the vector it holds, along one runA 50×50 diagonal matrix with one eigenvalue at 10⁻¹⁰ and the rest between 1 and 50, so every product with it is exact to a rounding and nothing below can be blamed on the operator. The lower curve is r ← r − αAp, updated by the recurrence, which is what a stopping test reads. The upper curve is ‖b − Ax‖/‖b‖, recomputed from the iterate. They start as the same vector and end 8.6·10⁸ apart, with the reported one at 6.89·10⁻²¹ — below the unit roundoff of 1.11·10⁻¹⁶, which is the shortest proof available that it is not the residual of anything.0122436486072849610⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best6.9·10⁻²¹the answer's, at its best5.9·10⁻¹²unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10⁹iterations drawn96the recurrence remembers every roundingand the stopping test is written in it
Fig. 28 Two residuals for one iteration, in the field that first drew the pair.
What the solver reported, and what it achievedFinal relative residual and final relative error against the working precision, on a logarithmic vertical axis. The residual sits at the stopping tolerance at every precision — the method converged, by every test available to it. The error runs from 1.4·10⁻¹³ at 53 bits to 0.00597 at 8, tracking the unit roundoff, which is drawn beside it.614223038465410⁻¹³10⁻¹⁰10⁻⁷10⁻⁴working significand bitsrelative sizeerrorresidualuthe gap the solver cannot seeerror ÷ residual at 24 bits1.5·10⁶error ÷ residual at 16 bits2.3·10⁷error ÷ residual at 8 bits2.4·10¹⁰every convergence test passesand the answer is wrong in proportion to u
Fig. 29 A residual that cannot see what it is asked about.
Backward error with and without the structure, for Levinson and for elimination, n = 12Both solvers are backward stable on this family by the usual measure: the smallest perturbation of any kind that explains either answer stays below 3.1·10⁻¹⁷ at every ρ, and the two curves are indistinguishable. Insisting the perturbation be a symmetric Toeplitz matrix — the kind of object the problem was posed with — lifts both by orders, and at ρ = 0.999 Levinson's rises by a factor of 2.48·10⁶. Neither number is alarming in absolute terms; what is worth carrying is that the reassuring one is the one that is reported, and the two are not measuring the same thing.10⁻³10⁻²10⁻¹110⁻¹⁹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹1 − ρbackward errorLevinson, Toeplitzelimination, ToeplitzLevinson, any kindelimination, any kindat ρ = 0.999Levinson, any perturbation3.1·10⁻¹⁷Levinson, Toeplitz only7.6·10⁻¹¹the ratio between them2.5·10⁶diagonal defect of the first0.97the number that is reportedand the number that is asked about
Fig. 30 A backward error measured against the structure of the problem rather than against any matrix.
Three residuals, and the one a solver returns is the one about nothingFor a rational approximant with m poles, on a target set reaching to 0.3: ‖T̃(λ)x‖ — the residual against the problem the eigensolver was handed — stays at 1.2·10⁻¹⁵ to 1.5·10⁻¹² — it does not fall with m, and past a point it slowly rises, because each added pole makes the fit's own basis worse conditioned; ‖T(λ)x‖, which costs one further evaluation of γ√(λ + c), falls with the approximation; and the forward error against the closed form falls with it, staying 20.5 to 24.8 times larger. The free number rises by 1306 across the sweep while the answer improves by 885. So the free number says nothing and the nearly-free number says almost everything. Evaluate the residual against the function you asked about, not against the one you handed over is the whole practical content of this field.23456710⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²poles in the approximantresiduals and errorforward error‖T(λ)x‖‖T̃(λ)x‖one extra evaluationagainst the approximant1.5·10⁻¹²against the problem asked1.9·10⁻⁶forward error4.8·10⁻⁵‖g − r‖ there8.5·10⁻⁵the free residual is flatand the answer is not
Fig. 31 The hero again, beside the approximants that produce it.

Where this goes next

0246810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ γ, the change of unitsrelative errorforward errorη, the quadraticη, the linearisationagainst a closed formη(linearisation), worst7.6·10⁻¹³η(quadratic), worst1.2·10⁻⁴forward error, worst0.0013coefficient spread4.2·10¹⁵the solver is right at every stopabout a problem nobody asked A backward-stable answer to a problem nobody asked One quadratic eigenvalue problem, in nine systems of units, with a change of variable that is exact in both directions. The residual the solver prints stays at the rounding level at every stop. The answer loses eleven orders of magnitude, and the two facts are consistent. the problem you posedA = H10b = A·(1, 2, …, 10)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 2.3·10⁻¹⁷the answer you wantedx = (1, 2, …, 10), exactlythe answer you gotx̂, wrong by 2.7·10⁻⁴ relativebackward error 2.3·10⁻¹⁷forward error 2.7·10⁻⁴κ = 1.6·10¹³κ · η = 3.6·10⁻⁴, and the measured forward error is 2.7·10⁻⁴.The algorithm is not at fault. The problem is.H10, LU with partial pivotingresidual and error differ The exact answer to a nearby problem A good algorithm does not give an approximate answer to your problem. It gives the exact answer to a problem very close to yours — and once that is the definition, a wrong result has two possible authors and they can be measured apart. -2-101234-5-3-1135real partimaginary part4 insidea countable spectruminside the contour4drawn12existing∞worst branch residual1.6·10⁻¹⁵there is no last eigenvalueso the question has to change A problem with infinitely many eigenvalues Let the matrix depend on λ through something that is not a polynomial and three things stop being true at once. There is no linearisation, there is no characteristic polynomial, and "compute the spectrum" is not a request that can be granted — the only finite question is how many eigenvalues are inside this circle.

What links here

Computed from the collection, not written here: the essays that point at this one.

Reads more easily once this is understood

Essays that name this one as worth reading first.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Approximation before linearisationBackward errorCondition numberExact ground truthForward errorLinearisationNonlinear eigenvalue problemRational approximationResidual