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

A perturbation that moves every coefficient

The backward error of a polynomial eigenpair is measured against perturbations of all three coefficients at once. Restrict it to the one coefficient anybody is willing to move and the same computed answers are stable at one eigenvalue and unstable at another, by a factor that runs from 1.06 to 6,370 across a single spectrum.
33 min read 26 figures Whose fault is itWhere the algebra stops being the arithmeticMeasured, not assertedAssertions that reject

Worth reading first: The exact answer to a nearby problem · A matrix that depends on its own eigenvalue · The units the matrix is measured in.

The backward error of a computed eigenpair of λ²M + λC + K is the smallest perturbation of the polynomial that makes the pair exact. Every published formula computes the same one:

η  =  ‖P(λ)x‖ ÷ ( (|λ|²‖M‖ + |λ|‖C‖ + ‖K‖) · ‖x‖ )

with the three coefficient norms in the denominator, added together. That number answers a specific question — how far is this from exact, if all three coefficients are equally open to being wrong — and it is the right question when all three are data.

The same 16 answers, backward stable under one reading and not under anotherEvery real eigenpair of an overdamped chain of 8 masses, computed once through the first companion linearisation, and its backward error measured four ways — against perturbations of all three coefficients, of K alone, and of M alone. The unrestricted number is the one every published formula computes and it sits at the rounding level across the spectrum. Restricting the perturbation charges the same residual to less data, so it can only raise the quotient, and the factor has a closed form: (|λ|²‖M‖ + |λ|‖C‖ + ‖K‖) divided by the allowed part, checked here against the two computed numbers to 10⁻⁸. Because it carries |λ|², it is not a property of the polynomial: allowing only K to move costs between 1.06 and 45.69, and allowing only M costs between 1.99 and 6367. At a tolerance of 10⁻¹⁵ 15 of these pairs are stable against the whole polynomial and 7 against K alone.10⁻²10⁻¹110¹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹|λ|backward error of the computed pairone residual, four denominatorspairs measured16worst, all three1.4·10⁻¹⁵worst, K alone1.8·10⁻¹⁴worst, M alone9·10⁻¹²K-only factor, low1.1high46solid: every coefficient may movedashed: only one may
Fig. 1 Every real eigenpair of an overdamped chain, with its backward error measured against perturbations of all three coefficients, of K alone, and of M alone.

They are frequently not all data.

Where the coefficients come from

A quadratic eigenvalue problem from a damped structure has three matrices with three different histories.

The three coefficient norms are the same quantities the scaling essay moves apart, and their provenance differs.

M is a discretisation. It is known to whatever accuracy the mesh has, its entries are integrals of basis functions, and it is often exactly right in the sense that matters — a lumped mass matrix is a diagonal of known quantities.

K comes from measured material properties. Moduli, thicknesses, stiffnesses. It has real uncertainty in it, of a size somebody could state.

C is very often a model. Rayleigh damping sets C = αM + βK for two scalars chosen to match two observed decay rates. Modal damping assigns a ratio per mode. In neither case is C a measurement of anything; it is a device.

So asking “how far is the computed answer from exact under a perturbation of all three” prices a perturbation of C that nobody would accept as a description of the physics — the perturbed C would not be a Rayleigh damping matrix at all, and a nearby problem of the wrong kind is the essay about exactly that objection in the structured setting.

The restricted question is: how far from exact, if only the coefficient one is willing to move may move?

The ratio has a closed form and it is not a constant

Charging the same residual to less data can only raise the quotient, so a restricted backward error is never smaller. How much larger has a closed form:

η_K / η_all  =  (|λ|²‖M‖ + |λ|‖C‖ + ‖K‖) ÷ ‖K‖

which contains |λ|². So it is not a property of the polynomial — it is a different number at every eigenvalue.

Measured on an eight-mass overdamped chain, whose sixteen eigenvalues run from −0.01997 to −7.417 and whose coefficient norms are 2.83, 19.89 and 6.78:

restriction    factor at the smallest |λ|    at the largest
K only               1.059                       45.7
M only               1.99                     6,370

The K-only restriction costs six per cent at one end of the spectrum and a factor of forty-six at the other. The M-only restriction costs a factor of two at one end and three and a half orders at the other.

Every one of those factors is checked against the two backward errors computed separately, and the closed form matches to 10⁻⁸ at every eigenvalue. It is a formula being verified rather than quoted.

The perturbation is built, not bounded

A backward error formula is a ratio of two norms, and a ratio of two norms is a claim about a perturbation that exists. This field builds it.

For an allowed set of coefficients, the minimal perturbation puts

ΔX  =  −(α_X / d) · r xᵀ / ‖x‖²

with r the residual, α_X the coefficient’s weight and d the sum of the allowed weights. The perturbed polynomial then annihilates the pair exactly, and the library assembles it, evaluates it, and checks.

The check has a floor and the floor is instructive. Evaluating λ²Mx + λCx + Kx in binary64 commits a rounding of about u·(|λ|²‖M‖ + |λ|‖C‖ + ‖K‖)·‖x‖ — which is exactly the quantity the backward error divides by. So the residual of the perturbed polynomial cannot be driven below one unit roundoff’s worth of that denominator, and a check written as a ratio to the original residual reports 7.7·10⁻⁴ on a construction that is exact.

Measured against the floor instead, the worst case over thirty-six constructions — three restrictions on twelve pairs — is 0.34 floors. The construction is exact to the level at which its own correctness can be evaluated, which is the honest statement and is what the assertion checks.

That distinction is worth carrying. A verification whose tolerance is “zero” fails on every exact construction whose verification involves cancellation, and the useful tolerance is the arithmetic’s own resolution at that magnitude — which is computable rather than guessed.

Why the weights are the norms and not something else

The denominator of the standard formula is a weighted sum of coefficient norms rather than a plain one, and the weights — |λ|², |λ|, 1 — are worth deriving rather than accepting, because the whole essay is about what happens when the sum is truncated.

The perturbed polynomial is (λ²(M + ΔM) + λ(C + ΔC) + (K + ΔK))x = 0, so the perturbations enter multiplied by λ², λ and 1 respectively. A perturbation of M of size δ therefore contributes |λ|²δ to the residual it can absorb; the same-sized perturbation of K contributes δ. So at a large eigenvalue, M is the coefficient with leverage, and at a small one it is K.

That is why the M-only ratio spans three and a half orders where the K-only one spans one and a half. The M-only ratio is (|λ|²‖M‖ + |λ|‖C‖ + ‖K‖)/(|λ|²‖M‖), which at the small end of the spectrum is dominated by ‖K‖/(|λ|²‖M‖) — a ratio of 6,370 when |λ| is 0.02. Restricting to a coefficient with no leverage at that eigenvalue is charging the residual to something that could barely have caused it.

The general rule, and it is the operational form of the whole essay: restrict to the coefficient that has leverage at the eigenvalue in question, or accept a factor. For the low modes of a damped structure that is K; for the high ones it is M; and a single restriction applied across a whole spectrum is the wrong restriction at one of its ends.

The same answers, two verdicts

A backward error is read against a tolerance. So the practical consequence of the ratios above is that one set of computed answers is accepted and rejected by two readings of the same number.

At a change of units of 10⁴, on the same sixteen pairs:

tolerance    stable, all three    without C    K only    M only
10⁻¹³             0                  0           0         0
10⁻¹¹             8                  7           0         7
10⁻⁹             16                 16          16        10

At 10⁻¹¹ eight pairs are stable under the whole polynomial and none is stable under K alone. At 10⁻⁹ every pair is stable under all three and under K alone, and six are not under M alone.

Nothing about the computation differs between the columns. The same eigenpairs, from the same linearisation, from the same run. What differs is which perturbations the reader is willing to accept as a description of the problem — a modelling decision — and it moves the verdict on half the spectrum.

The same 24 answers, backward stable under one reading and not under anotherEvery real eigenpair of an overdamped chain of 12 masses, computed once through the first companion linearisation, and its backward error measured four ways — against perturbations of all three coefficients, of K alone, and of M alone. The unrestricted number is the one every published formula computes and it sits at the rounding level across the spectrum. Restricting the perturbation charges the same residual to less data, so it can only raise the quotient, and the factor has a closed form: (|λ|²‖M‖ + |λ|‖C‖ + ‖K‖) divided by the allowed part, checked here against the two computed numbers to 10⁻⁸. Because it carries |λ|², it is not a property of the polynomial: allowing only K to move costs between 1.03 and 45.59, and allowing only M costs between 1.99 and 2.664·10⁴. At a tolerance of 10⁻¹⁵ 24 of these pairs are stable against the whole polynomial and 14 against K alone.10⁻³10⁻²10⁻¹110¹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹|λ|backward error of the computed pairone residual, four denominatorspairs measured24worst, all three6.1·10⁻¹⁶worst, K alone2.6·10⁻¹⁴worst, M alone4.9·10⁻¹²K-only factor, low1high46solid: every coefficient may movedashed: only one may
Fig. 2 A longer chain, whose wider spectrum spreads the restricted and unrestricted backward errors further apart.

Why the change of units does not change the ratio

One measurement went the way it had to and is worth recording because it confirms the arithmetic rather than the modelling.

The scaling essay in this field is about a change of units γ that multiplies M by γ², C by γ and leaves K alone — which moves the coefficient norms apart by decades and is worth ten orders in the forward error of a badly scaled problem.

It does not move the ratio here at all. At γ = 1, 10² and 10⁴ the K-only factor spans 1.059 to 45.686 identically, to every digit. The reason is that the ratio is a quotient of weights, γ scales the weights and the eigenvalues in compensating ways, and the whole expression is invariant.

That is the right behaviour and it is a check on the measurement: a quantity that is supposed to be about the structure of the perturbation should not move when the units do. A condition number scaling cannot move is the essay about quantities that fail this test, and this one passes it.

What to report

The restricted backward error is not a new quantity to compute — it is the standard one with a different denominator, so it costs a division. What it needs is a decision about what goes in the denominator, and that decision belongs to whoever knows where the matrices came from.

State the allowed set. “Backward stable at 10⁻¹⁴” is incomplete. “Backward stable at 10⁻¹⁴ against perturbations of M, C and K” and “at 4·10⁻¹³ against perturbations of K alone” are two statements about one run, both true, and only the second answers a question about the structure.

Report the factor, not only the number. The ratio between the restricted and unrestricted errors is available in closed form from three norms and one eigenvalue, before either error is computed, and it says which end of the spectrum a restriction will bite at.

And be careful reading a tolerance. The table above is one set of answers passing at one reading and failing at another. A tolerance chosen against the unrestricted number and applied to a problem where only one coefficient is free is a tolerance that is too loose by up to three orders — at one end of the spectrum, and by six per cent at the other.

What a structured restriction is not

It is worth marking the boundary between this measurement and the structure field’s, because the two answer questions that sound the same.

A structured backward error asks for the smallest perturbation of a given form — Toeplitz, symmetric, positive — that makes the pair exact. Its denominator is the same; what changes is that the perturbation is restricted to a subspace of matrices, and computing it requires an optimisation over that subspace rather than a formula.

The restriction here is coarser and cheaper. It says nothing about the shape of ΔK, only that ΔM and ΔC are zero. So it is a restriction to a coordinate subspace of the three-coefficient space, which is why the minimal perturbation still has a closed form and why the ratio is a quotient of weights rather than the output of a minimisation.

Those two are composable and this field measures only one of them. A perturbation restricted to K and required to be symmetric is a further restriction, its backward error is at least this one, and the structure field’s measurement says the further factor is usually small — the structure never hurts by more than a modest constant. So the ordering of the three quantities is known even where only two are computed, which is the useful thing to take from having both fields.

Where this sits among the site’s backward errors

Three quantities are now in play for one computed pair, and it is worth separating them.

The linearisation’s backward error asks how far the computed pair is from exact for the 2n × 2n matrix the solver actually ran on. A backward-stable answer to a problem nobody asked is the essay about it, and the answer is that it is always tiny and always about the wrong object.

The quadratic’s backward error asks the same about the polynomial, unrestricted. That is the standard published quantity and it is what most of this field measures.

The restricted backward error asks it against a subset of the coefficients. It is never smaller, it is not a fixed multiple of the second, and the multiple is a function of the eigenvalue.

The same three-way separation appears in the approximation committed before the arithmetic, where the third quantity is not a backward error at all. Three numbers, one computation, and each answers a different question about what would have to be wrong for the answer to be right. The site’s whole vocabulary for error rests on the observation that a wrong answer has an author; this is the observation that who counts as an author is a choice, and that the choice changes the verdict on half a spectrum.

The refusals, and why there are three

Three assertions in this library are fed cases they must reject, and each one closes off a reading that a careless implementation would take.

A backward error with no coefficient allowed to move. The denominator is then zero, and the correct answer is that the question is malformed rather than that the error is infinite. A routine that returned infinity would be reporting a number where there is none.

The claim that a restriction can lower a backward error. Charging one residual to less data can only raise the quotient, so η_restricted ≥ η_all always. That is trivial to prove and easy to violate in code — a denominator assembled from the wrong weights, or a residual recomputed against the restricted polynomial rather than the original — and the assertion is fed a spectrum and required to find no counterexample.

And the claim that the restriction costs a fixed factor. This is the essay’s own finding turned into a refusal: a routine that computed the ratio once and applied it across a spectrum would be wrong by three orders at one end, and the assertion is handed a spectrum whose ratios span forty-three and required to reject the claim that they are all within one per cent of each other.

Each of those runs on every build. The second and third are the interesting ones, because both would otherwise be invisible: a restricted backward error that is too small looks like good news, and a constant factor looks like a simplification. Neither throws, neither is caught by a residual, and both are the kind of claim this site’s habit exists to make impossible to keep.

One line to add to a solver’s output

The cheapest form of everything above is a single extra line, and it needs three norms the solver has already computed and one eigenvalue it has just produced.

restriction factor to X  =  (|λ|²‖M‖ + |λ|‖C‖ + ‖K‖) ÷ (the allowed part)

Printed beside each eigenpair, it tells a reader immediately how much a restriction would cost at that eigenvalue — before either backward error is computed, and without any decision about which restriction is the right one. A reader who knows their damping matrix is a model reads the M-and-K column; one who trusts everything reads none of them.

That is the site’s habit applied to a modelling decision rather than to an arithmetic one: print the number that says how much the decision is worth, and let the person who knows the answer make it.

At other settings

The same 8 answers, backward stable under one reading and not under anotherEvery real eigenpair of an overdamped chain of 4 masses, computed once through the first companion linearisation, and its backward error measured four ways — against perturbations of all three coefficients, of K alone, and of M alone. The unrestricted number is the one every published formula computes and it sits at the rounding level across the spectrum. Restricting the perturbation charges the same residual to less data, so it can only raise the quotient, and the factor has a closed form: (|λ|²‖M‖ + |λ|‖C‖ + ‖K‖) divided by the allowed part, checked here against the two computed numbers to 10⁻⁸. Because it carries |λ|², it is not a property of the polynomial: allowing only K to move costs between 1.19 and 45.73, and allowing only M costs between 2 and 717.5. At a tolerance of 10⁻¹⁵ 8 of these pairs are stable against the whole polynomial and 3 against K alone.10⁻²10⁻¹110¹10⁻¹⁷10⁻¹⁵10⁻¹³|λ|backward error of the computed pairone residual, four denominatorspairs measured8worst, all three8.7·10⁻¹⁶worst, K alone1.8·10⁻¹⁴worst, M alone2.4·10⁻¹³K-only factor, low1.2high46solid: every coefficient may movedashed: only one may
Fig. 3 A four-mass chain, whose spectrum is narrow and whose restriction factors are correspondingly close together.
The same 12 answers, backward stable under one reading and not under anotherEvery real eigenpair of an overdamped chain of 6 masses, computed once through the first companion linearisation, and its backward error measured four ways — against perturbations of all three coefficients, of K alone, and of M alone. The unrestricted number is the one every published formula computes and it sits at the rounding level across the spectrum. Restricting the perturbation charges the same residual to less data, so it can only raise the quotient, and the factor has a closed form: (|λ|²‖M‖ + |λ|‖C‖ + ‖K‖) divided by the allowed part, checked here against the two computed numbers to 10⁻⁸. Because it carries |λ|², it is not a property of the polynomial: allowing only K to move costs between 1.1 and 45.73, and allowing only M costs between 2 and 2449. At a tolerance of 10⁻¹⁵ 12 of these pairs are stable against the whole polynomial and 5 against K alone.10⁻²10⁻¹110¹10⁻¹⁷10⁻¹⁵10⁻¹³|λ|backward error of the computed pairone residual, four denominatorspairs measured12worst, all three6.7·10⁻¹⁶worst, K alone2·10⁻¹⁴worst, M alone6.8·10⁻¹³K-only factor, low1.1high46solid: every coefficient may movedashed: only one may
Fig. 4 Six masses, where the spread has opened by an order.
The same 20 answers, backward stable under one reading and not under anotherEvery real eigenpair of an overdamped chain of 10 masses, computed once through the first companion linearisation, and its backward error measured four ways — against perturbations of all three coefficients, of K alone, and of M alone. The unrestricted number is the one every published formula computes and it sits at the rounding level across the spectrum. Restricting the perturbation charges the same residual to less data, so it can only raise the quotient, and the factor has a closed form: (|λ|²‖M‖ + |λ|‖C‖ + ‖K‖) divided by the allowed part, checked here against the two computed numbers to 10⁻⁸. Because it carries |λ|², it is not a property of the polynomial: allowing only K to move costs between 1.04 and 45.64, and allowing only M costs between 1.99 and 1.385·10⁴. At a tolerance of 10⁻¹⁵ 19 of these pairs are stable against the whole polynomial and 8 against K alone.10⁻²10⁻¹110¹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹|λ|backward error of the computed pairone residual, four denominatorspairs measured20worst, all three3·10⁻¹⁵worst, K alone2.5·10⁻¹⁴worst, M alone4.2·10⁻¹¹K-only factor, low1high46solid: every coefficient may movedashed: only one may
Fig. 5 Ten, where the M-only curve has separated from the other two at both ends.
The same 28 answers, backward stable under one reading and not under anotherEvery real eigenpair of an overdamped chain of 14 masses, computed once through the first companion linearisation, and its backward error measured four ways — against perturbations of all three coefficients, of K alone, and of M alone. The unrestricted number is the one every published formula computes and it sits at the rounding level across the spectrum. Restricting the perturbation charges the same residual to less data, so it can only raise the quotient, and the factor has a closed form: (|λ|²‖M‖ + |λ|‖C‖ + ‖K‖) divided by the allowed part, checked here against the two computed numbers to 10⁻⁸. Because it carries |λ|², it is not a property of the polynomial: allowing only K to move costs between 1.02 and 45.56, and allowing only M costs between 1.99 and 4.68·10⁴. At a tolerance of 10⁻¹⁵ 28 of these pairs are stable against the whole polynomial and 15 against K alone.10⁻³10⁻²10⁻¹110¹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹|λ|backward error of the computed pairone residual, four denominatorspairs measured28worst, all three6.7·10⁻¹⁶worst, K alone1.3·10⁻¹⁴worst, M alone1.9·10⁻¹¹K-only factor, low1high46solid: every coefficient may movedashed: only one may
Fig. 6 And fourteen, the largest chain the linearisation is afforded here.
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. 7 The two backward errors this essay adds a third to.
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. 8 The spectrum whose two ends the restriction factor differs by orders between.
The two perturbations that make one computed solution exact, on a 8×8 Toeplitz systemA Kac–Murdock–Szegő matrix at ρ = 0.95, solved by Levinson's recursion. Every backward-error claim on this site says the computed answer solves a nearby problem exactly, and the nearby problem is the top matrix: the smallest perturbation of any kind, 3.36·10⁻¹⁷ relative, rank one, and constant along 0.98 of the way to none of its diagonals. The bottom matrix is the smallest perturbation that is itself a symmetric Toeplitz matrix — the same kind of object the problem was posed with — and it is 1.7·10⁻¹⁴, larger by a factor of 505. Both explain the same computed answer exactly. Only one of them is a problem anybody could have posed.the smallest perturbation of any kind — 3.36·10⁻¹⁷the smallest Toeplitz one — 1.7·10⁻¹⁴both exact for the same x̂smallest of any kind3.4·10⁻¹⁷smallest Toeplitz one1.7·10⁻¹⁴the price of the constraint505diagonal defect, unconstrained0.98an exact answer to a nearby problemof a kind nobody posed
Fig. 9 A backward error measured against a structure rather than against a matrix, in the structure field.
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. 10 The same question asked of a Toeplitz solver.
The support points found a branch point nobody named, and cluster at it by 1.52Approximating γ√(λ + c) on an interval reaching to within 0.03 of its branch point, by a barycentric rational whose support points are chosen greedily: at each step the sample where the current approximation is worst becomes the next support point. The filled curve is the error against degree; the open one is the same barycentric form with support points spread evenly along the interval, which is the choice anybody makes who has no reason to make another. At degree 10 the adaptive rule is 702.3× better, and the advantage grows with the degree. Nothing in the algorithm is told that the function has a branch point or where it is: the points it selects lie between 0.0962 and 9.53 from it, a span of 99.06, with a mean ratio of 1.519 between consecutive distances — the geometric clustering the hand-built approximant in this field is given as a rule, arrived at from the residual.12345678910111210⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹degree of the approximantworst relative error on the target setchosen by the residualdegree reached12error there10⁻¹²even support, same degree1.1·10⁻⁷advantage702nearest support point0.096clustering ratio1.5filled: points chosen by the erroropen: points spread evenly
Fig. 11 The other error committed before the arithmetic, in the previous essay of this field.
The backward error of a hierarchical solve, against the error of the representation it usedSolve A_H x = b exactly and the residual against the matrix that was wanted is b − Ax = (A_H − A)x, so ‖b − Ax‖ ⁄ ‖A‖‖x‖ cannot be anything but the representation's own error. Measured across six accuracies spanning ten decades, the two track at a slope of 0.994 and sit a constant 6.7× apart, which is the difference between a norm of a matrix and a norm of that matrix applied to one vector. The reading is the one this field is built on: the accuracy knob is not an accuracy, it is a backward error chosen in advance. Everywhere else on this site a backward error is something an algorithm produced and somebody then measured; here it is a line in the program, and its size is known before the solve starts.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.99representation at 10⁻⁸1.2·10⁻⁹backward error there1.3·10⁻¹⁰representation at 10⁻¹²1.2·10⁻¹³backward error there2.5·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem
Fig. 12 A backward error chosen rather than incurred.
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: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 13 The two errors, of which this essay refines the first.
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. 14 The identity all three of these backward errors are read through.
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⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 15 The amplifier that turns any of them into a forward error.
Two perturbation bounds and the error that was measured, on a 8×8 matrix spread over 8 decades of unitsThree curves against the size of an entrywise relative perturbation. The normwise bound κ_∞·ε is a valid bound and sits 4·10⁷ times above the componentwise one cond(A, x)·ε, which is also a bound and is nearly attained by the worst of forty random perturbations at each size.10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ_∞ · εcond(A,x) · εmeasuredboth bounds holdκ_∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement
Fig. 16 An error measured entry by entry, which is a third kind of restriction.
Two condition numbers of one 8×8 system, as its rows are put into different unitsFour curves against the spread of the row units, in decades. κ_∞ of the scaled matrix rises from 9.83 to 1.9·10¹⁰ while the componentwise condition number stays at 6.98 throughout — the same system, the same solution, and one of the two numbers is a fact about the units. Hilbert's two numbers are drawn flat beside them at 3.4·10¹⁰ and 1.2·10¹⁰: a matrix whose sensitivity no scaling repairs.0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it
Fig. 17 Coefficient norms moved apart by a change of units, which this ratio is invariant to.
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. 18 A scaling recovering accuracy, in the field that measured what one costs.
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. 19 An answer known exactly, which the chain’s spectrum also is.
The accuracy asked for, against the accuracy obtained, for two kernels on one partitionε is applied to each block against that block's own largest singular value; the error is then reported against the whole matrix's norm. The two are not the same number and the dashed diagonal is where they would be. For 1/r the obtained error is 29 times smaller than the tolerance asked for; for log r on the identical partition it is 13 times smaller. Both curves are straight and parallel to the diagonal, so the knob does what a knob should — a decade in buys a decade out — and neither of them sits on it. The factor is how much of the matrix's mass lives off the diagonal, which is a property of the kernel; it is not the number of blocks, which is 112 and would put the curves on the other side of the diagonal.10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³ε asked for, per block‖A − A_H‖ ⁄ ‖A‖ obtainedasked = obtainedlog r1 ⁄ rtwo numbers, not one1/r, obtained at 10⁻⁸3.4·10⁻¹⁰log r, obtained at 10⁻⁸7.8·10⁻¹⁰1/r, obtained ⁄ asked0.034log r, obtained ⁄ asked0.078blocks in the partition112the tolerance is per blockand the error is per matrix
Fig. 20 A tolerance asked for against one obtained, which the verdict table is about reading.
How close Hager's estimate is to the true κ₁, over 200 seeded 8×8 matricesFive bars. The estimate is exactly the true condition number on 81% of the sample and inside ten per cent on 87%; the worst underestimate in the whole sample returns 38% of the truth. The last bar is the matrix built to defeat it, at 7.7%, well below anything the random sample reached.each bar is a percentage — of the sample, or of the true condition numberexactly right80.5%inside 10%87.0%inside a factor of 287.0%worst in the sample, ×10037.7%the constructed matrix, ×1007.7%usually exactexact share0.81worst of the sample0.38the constructed matrix0.077a routine that is right most of the timeand never wrong in the safe direction
Fig. 21 An estimate reported with its spread.
Past the cliff: what each method returnsLoss of orthogonality at four condition numbers past κ²u = 1. Cholesky QR returns a factor whose implied Q is 0.37 away from orthogonal at κ = 10⁸ and 1.49 at 10¹² — and at 3·10¹⁰ it refuses outright, on a pivot of -6.5·10⁻¹⁷. The sweep and the tree return a usable factorisation of every one of these matrices.‖QᵀQ − I‖ of the implied Qκ = 10⁸, Cholesky0.37κ = 10⁸, sweep8.5·10⁻⁹κ = 10¹⁰, Cholesky1.3κ = 10¹⁰, sweep3·10⁻⁷κ = 3·10¹⁰, Choleskyrefusedκ = 3·10¹⁰, sweep6.7·10⁻⁶κ = 10¹², Cholesky1.5κ = 10¹², sweep1.4·10⁻⁴the safe run is the one that failsrefusals in the range1wholly non-orthogonal returns2the pivot it refused on-6.5·10⁻¹⁷one of these outcomes is safeand it is the refusal
Fig. 22 A routine that returns something rather than refusing, which a backward error always does.
Forward error of a 6×6 Hilbert solve at eight precisionsA bar for each significand width from 12 to 53 bits showing the relative error in the computed solution, with the condition number times the unit roundoff marked as a prediction.κ = 1.5·10⁷ · the exact answer is (1, 2, …, 6)12 bits3.316 bits0.8620 bits4.324 bits0.09330 bits7.6·10⁻⁴36 bits6.2·10⁻⁶43 bits2·10⁻⁷53 bits4.5·10⁻¹¹dashed: κ · unit roundoffone matrix, eight arithmeticsmeasured against a known answer
Fig. 23 The precision knob, which moves all three backward errors together.
Growth factor under partial pivoting to n = 40: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8
Fig. 24 A bound never attained, in the field that measured it first.
Singular values of a rank-4 matrix with noise of relative size 10⁻⁸Ten singular values on a logarithmic axis. The first four sit near one; the rest sit at the noise level, and the vertical distance between the two groups is the evidence for the rank.1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 25 A threshold turning a number into a verdict, which a tolerance here also is.
Six is the degree, and the model reaches it from 24 numbersA reduced model of a forty-state system built out of 24 samples of its transfer function and nothing else — no state matrix, no eigenvalues, no factorisation of anything the model is of. The Loewner pencil of divided differences is formed from the samples, its singular values are read for the degree (6 here, with a gap of 2.01·10⁸ across the cliff), and the pencil is projected onto that many directions. The filled curve is the worst relative error at the samples the model was built from; the open one is the worst error at forty points it never saw. At the true degree of six the first is 6.99·10⁻¹² — interpolation is exact by construction, not by accuracy — and the second is 2.12·10⁻⁹.12345678910111210⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹order of the reduced modelworst relative errora model made of measurementssamples used24degree read6gap at the cut2·10⁸best model, order12its error2.7·10⁻¹¹states in the original40filled: at its own samplesopen: everywhere else
Fig. 26 A model built from measurements, in the next essay of this phase.

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 smallest perturbation of any kind — 1.89·10⁻¹⁷the smallest Toeplitz one — 2.65·10⁻¹³both exact for the same x̂smallest of any kind1.9·10⁻¹⁷smallest Toeplitz one2.7·10⁻¹³the price of the constraint1.4·10⁴diagonal defect, unconstrained0.97an exact answer to a nearby problemof a kind nobody posed A nearby problem of the wrong kind A good algorithm returns the exact answer to a nearby problem. A hundred and eighteen essays have measured the distance and not one has asked what the nearby problem looks like. On a Toeplitz system it is a rank-one matrix that is constant along none of its diagonals — and the smallest one that is Toeplitz is two and a half million times larger. 10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.99representation at 10⁻⁸1.2·10⁻⁹backward error there1.3·10⁻¹⁰representation at 10⁻¹²1.2·10⁻¹³backward error there2.5·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem An accuracy that is a backward error Every backward error on this site is something an algorithm produced and somebody then measured. This one is a line in the program. Solving with a compressed matrix gives a residual that is the compression's own error, at a slope of 1.000 over ten decades, so the knob that sets the storage sets the backward error directly.

What links here

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

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.

Backward errorCoefficient normsConditioningLinearisationModelling assumptionPolynomial eigenvalueStructured backward errorTolerance