Numerical Linear Algebra
where the algebra stops being the arithmetic
Two errors, and whose fault they are

Three errors and one number

This site's identity has two factors and a division of blame between them. Two fields have now added a third party and a fourth, and only one of the four is a property of anything — the others are decisions, made before the arithmetic, reported by nothing.
33 min read 29 figures Whose fault is itWhere the algebra stops being the arithmeticMeasured, not asserted

Worth reading first: The exact answer to a nearby problem · An index that is a pair · Elimination is a sequence of choices.

The exact answer to a nearby problem is the first essay of this field and it states the identity the whole site is built on:

forward error  ⪅  condition number × backward error.

A wrong answer has two possible authors and they can be measured separately. The problem was sensitive, or the algorithm was not stable, and only their product is what anyone sees. Eleven fields have since applied that division and it has held.

Two of the newest fields have now put something in front of it, twice, and the additions are of a kind the original two are not.

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. 1 The identity as this field first drew it: two factors and one visible product.

The third party

Six routes to one spectrum measures a quadratic eigenvalue problem across eight decades of an exact change of variable. Across that sweep:

  • the problem’s condition number is 4.98 at every stop, flat;
  • the algorithm’s backward error against the matrix it factorised is 7.6·10⁻¹³, flat;
  • the answer loses eleven orders.

Neither factor moved and the answer collapsed. What moved is a matrix nobody chose to look at: the linearised pencil’s own condition number, which runs from 49.7 to 7.8·10¹¹ over the same sweep. The substitution from a quadratic to a 2n × 2n matrix is exact in the algebra and is not neutral in the arithmetic, and the loss belongs to it. That is the reading a nearby problem of the wrong kind gives to a structured problem whose backward error is measured over the wrong set of perturbations: the number is right about a question nobody asked.

That is a third author, and naming it changes what a diagnosis has to include. The problem was fine and the algorithm was fine is no longer a complete sentence, because there is a step between them that has a conditioning of its own.

The fourth

The problem the solver was actually given measures a nonlinear eigenvalue problem where there is nothing to linearise at all — a finite pencil has finitely many eigenvalues and a square root has a branch cut. Every algorithm replaces the function first, then linearises the replacement.

For a computed eigenpair of the replacement:

  • the residual against the object the solver factorised is 8·10⁻¹⁵;
  • the residual against the problem that was asked is 3.0·10⁻⁵;
  • the forward error is 5.9·10⁻⁴.

The first number does not fall when the answer improves. It slowly rises, because a better approximation is a larger matrix. So the number a solver returns moves in the opposite direction to the accuracy of what it returned, over a sweep in which the answer improves by four orders.

The identity with a term in front:

forward error  ⪅  ‖g − r‖  +  condition number × backward error.

The new term is not amplified by anything, not produced by the arithmetic, and not a property of the problem.

What kind of thing each of the four is

Setting them out together is the point of this essay, because the differences matter more than the list.

author what it is who chose it visible in a residual
the problem’s conditioning a property of the question nobody no, but estimable
the algorithm’s backward error a property of the code the library author yes, directly
the reformulation a property of a substitution the library author, once no
the approximation a number somebody typed the caller, every time no

The first two are properties. Nobody decides a Hilbert matrix’s condition number, and a library’s backward stability is a fact about it that a user inherits.

The third is a decision made once, by whoever wrote the solver, and inherited by everyone who calls it. A user of a quadratic eigenvalue solver does not choose the linearisation and usually does not know which one is inside.

The fourth is different in kind. It is chosen every time, by the caller, usually without being told a choice is being made — the target set, the degree, the region. And it is the largest of the four on the problems where it appears: five orders from the region, against a factor of forty from the linearisation and nothing at all from the conditioning.

Why a residual cannot see two of them

The mechanism is the same in both cases and it is worth stating once.

A residual is a backward error for the object the arithmetic was performed on. That is exactly what makes it valuable: it is computable without knowing the answer, and it certifies that the computation did what it claimed on what it was given.

When the object handed to the arithmetic is the problem, up to an exact transformation, the residual is a statement about the problem. Every field before these two was in that situation, which is why a small residual is not a small error needed only one extra quantity — the condition number — to complete the account.

When the object is a reformulation, the residual is a statement about the reformulated object, and the map back to the original is not norm-preserving. When it is an approximation, the residual is a statement about a different function. In both cases the word backward stops pointing where it used to, and the number is correct, cheap, and about something else.

What to compute instead, in each case

Both have a repair and both repairs are cheap, which is the useful half of this essay.

For the reformulation: compute the backward error of the computed eigenvalue with respect to the original coefficients rather than to the linearised matrix. For a quadratic that is η = ‖Q(λ)x‖ / ((|λ|²‖M‖ + |λ|‖C‖ + ‖K‖)‖x‖), two matrix–vector products. It is the quantity that tracks the forward error in six routes to one spectrum: its ratio to the forward error stays within a factor of 5.8 across the whole sweep while the linearised residual’s moves by 1.5·10⁷.

For the approximation: compute the residual against the original function, which is one further evaluation of g per computed eigenvalue. It tracks the forward error to a factor of about twenty at every degree.

And in both cases: compute the term itself in advance. The linearisation’s conditioning is available from norms already in hand; ‖g − r‖ is a maximum of a scalar function over an interval. Both are free relative to the factorisation at the end of the chain, and neither is reported by anything — the same gap an estimate that can be fooled records for the one quantity of the four that does have an estimator.

The pattern, stated as a habit

Four authors, and the two that are hardest to see are the two that were chosen. That inversion is the thing worth carrying.

A property can be estimated, and this site has spent eleven fields building estimators for them: condition estimators, backward-error formulas, growth factors. A decision cannot be estimated — it can only be reported, by whoever made it — and nothing in the conventional output of a solver reports one.

So the habit the two newest fields suggest is a question rather than a measurement: what was replaced, between the problem that was asked about and the object that was factorised? If the answer is nothing, the two-factor identity is complete and the existing tools apply. If the answer is a linearisation, the reformulation’s conditioning is a third quantity to look at. If it is an approximation, there is a term in front of the whole identity that no amount of care in the arithmetic will reduce.

Asking that question costs nothing and it is not answered by any number a solver currently returns.

Where the two new authors do and do not appear

It is worth saying which problems have three or four authors and which have two, because the answer is not “the hard ones”.

Two authors is the ordinary case and covers every field before polynomial. A linear system, a least-squares fit, an eigenvalue problem for a matrix, a sparse factorisation: the object the arithmetic runs on is the object that was asked about, possibly reordered or rescaled by transformations that are exact. Everything this site built in eleven fields applies.

Three appears whenever a problem is reformulated into a different kind of object before being solved. A matrix polynomial linearised into a pencil is the instance measured here, and it is not the only one: forming the normal equations is a reformulation in exactly this sense, and the road that squares the problem is that essay written before the vocabulary existed. So is forming PQ to get Hankel singular values, in the product nobody had to form. The pattern in each case is the same: the residual against the reformulated object is small, the reformulated object’s conditioning is worse than the original’s, and the gap is the reformulation’s.

Four appears when the problem cannot be reformulated exactly at all — when there is no finite object with the same answer, and the only route is to replace the question. That is rarer and it is what the nonlinear eigenvalue problem is. It also covers a case this site has met without naming: truncating an infinite series to compute a matrix function, where the series that has to be squared back measures a truncation error that sits in exactly this position.

The reason to sort them is that the diagnostic differs. With two authors, measure both. With three, measure the reformulated object’s conditioning as well, and expect the answer to be a change of variable. With four, the dominant term is usually the one that was chosen, and the repair is to choose differently rather than to compute more carefully.

What has not changed

Worth closing on, because a reader could take this essay as an argument that residuals are unreliable and it is the opposite.

Every measurement above is possible because the backward-error framework works. The reason the reformulation’s contribution could be isolated is that the linearised solve’s backward error was measured and found flat; the reason the approximation’s could be is that the residual against the approximant was measured and found at rounding. Both new authors were found by the existing tools being applied carefully enough to notice that they had run out of explanation.

The identity is not wrong. It has always been a statement about a specific problem and a specific algorithm, and what the newest fields supply is a reason to check which problem and which algorithm those are.

A worked reading, on one number

It helps to run the account once on a single computed quantity, because the four authors are easier to keep apart when they are attached to one answer rather than to a table.

Take the quadratic eigenvalue problem of a damped chain of eight masses, in units chosen badly by a factor of 10⁸, and ask for its spectrum. A library returns twenty-eight numbers and, if asked, a residual of about 7.6·10⁻¹³. The forward error against the closed form is 1.3·10⁻³.

Four questions, in the order worth asking them.

Was the problem sensitive? The eigenvalue condition number of the quadratic is 4.98. No. Whatever went wrong was not the question’s fault, and this is the measurement that rules out the answer most people reach for first.

Was the algorithm unstable? Its backward error against the matrix it factorised is 7.6·10⁻¹³, which for a 16 × 16 real Schur factorisation is what stability looks like. No.

Was something reformulated? Yes — a quadratic became a 2n × 2n pencil — and the reformulated object’s condition number is 7.8·10¹¹ where the quadratic’s is 4.98. That is the loss, located, and it is eleven orders in the right place.

Was anything approximated? No. The linearisation is exact and there is no fourth term here.

Total account: eleven orders, all of it in the third author, none of it visible in either of the two numbers a library would return. And the repair is two lines of scaling computed from three norms, after which the same four questions give the same first two answers and 1.6·10⁻¹³ for the third.

That is the shape a diagnosis should have. Four questions, each with a number, and a sum that closes.

Why this took eleven fields to notice

A closing observation about how the two new authors were found, since the method is more portable than the finding.

Neither was found by a gate. Every check this site runs asks whether a computed quantity satisfies a relation, and both new authors produce computations in which every relation holds: the residual is small, the factorisation reconstructs, the assertions pass. A suite that measures whether numbers are right cannot find an author of error that makes no number wrong.

Both were found by sweeping a parameter that was supposed to be free. The change of units in the quadratic field is exact — the spectrum divides by γ and the closed form follows it — so it is a control: anything that moves across the sweep is manufactured rather than inherent. The target set in the nonlinear field is the same kind of thing, a decision that changes no eigenvalue and changes the answer.

That is a technique rather than a result and it is worth stating as one. To find an error source nobody is measuring, vary something that is not supposed to matter and see what moves. Both of this site’s newest authors were sitting in plain sight behind a parameter everybody treated as notation.

The one line

Before reading a residual as an accuracy estimate, ask what was replaced between the problem and the object that was factorised. If nothing was, the number means what it appears to. If a linearisation was, there is a third quantity to measure and it is the conditioning of the thing that was substituted. If an approximation was, there is a term in front of the whole identity, it was chosen rather than incurred, and it is usually the largest of the four.

At other settings

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 = H12b = A·(1, 2, …, 12)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 1.8·10⁻¹⁷the answer you wantedx = (1, 2, …, 12), exactlythe answer you gotx̂, wrong by 0.02 relativebackward error 1.8·10⁻¹⁷forward error 0.02κ = 1.8·10¹⁶κ · η = 0.33, and the measured forward error is 0.02.The algorithm is not at fault. The problem is.H12, LU with partial pivotingresidual and error differ
Fig. 2 The identity on a larger problem, where the amplification is greater and the account still closes.
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. 3 A small residual and a large error, which is what the two-factor identity explains.
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. 4 The amplifier, which touches neither of the two new terms.
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. 5 The standing every one of these measurements needs.
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. 6 Two backward errors for one computed eigenvalue, in the polynomial field.
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. 7 And the condition number that stayed at 4.98 while the answer collapsed.
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. 8 Six routes to one spectrum, differing by a factor the residual cannot see.
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. 9 And the two lines that remove the whole loss.
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. 10 Three residuals for one eigenvalue, and the free one is flat.
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.2 to 8.99, and its left end is 0.2 from the branch point — a choice made before anything is computed. The rational approximant with 8 poles clustered on the cut reaches 8.46·10⁻⁴; the polynomial of degree 9, which linearises to a matrix of the same size, reaches 0.00676. 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.06left end, from the cut0.2rational, worst8.5·10⁻⁴polynomial, worst0.0068linearisation size, both54committed before the solveand invisible to it
Fig. 11 The term the free residual cannot contain.
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. 12 And what else a replaced problem returns.
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. 13 A backward error measured against a structure rather than against any matrix.
The two perturbations that make one computed solution exact, on a 10×10 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, 1.89·10⁻¹⁷ relative, rank one, and constant along 0.97 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 2.65·10⁻¹³, larger by a factor of 1.41·10⁴. Both explain the same computed answer exactly. Only one of them is a problem anybody could have posed.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
Fig. 14 The shape a perturbation is allowed to have, which decides what backward means.
One 40×40 matrix described four ways, with its condition number and its backward error at eachThe matrix is ρ^|i−j| at ρ = 0.9, where κ = 284.2. Described as n² entries its condition number is that; as 79 constant diagonals it is 145.9; as the 40 a symmetric Toeplitz matrix has, 140.5; and as the one number ρ that the matrix actually holds, 1.65. The backward error of the same computed solution runs the other way — 1.84·10⁻¹⁷, 7.29·10⁻¹⁵, 8.74·10⁻¹⁴ — and at the last rung there is none: no ρ whatever has the computed answer as its exact solution, and the closest one leaves 100 per cent of the residual unexplained. The drop that matters is the last one, and it is not a linear-algebra structure at all.110¹10²10³10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²10²10⁶numbers that describe the matrixcondition number, and backward errordensetoeplitzsymmetricρ aloneno such problemcondition numberbackward errorone matrix, four descriptionsκ, all n² entries348as one number, ρ1.7backward error, unconstrained1.8·10⁻¹⁷as a symmetric Toeplitz matrix8.7·10⁻¹⁴fewer numbers, better conditionedand no nearby problem left
Fig. 15 Trading a description for a perturbation, in the structure field.
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. 16 Units, which are the third author’s usual disguise.
Four pivot rules on [[ε, 1], [1, 1]] with ε = 10⁻¹², before and after its first row is scaled by 1/εEight bars of forward error. On the matrix as given, no pivoting loses the answer entirely and all three pivot rules return it to rounding. Multiply the first row by 1/ε — which does not change the solution — and partial pivoting makes no interchange at all and returns the same wrong answer that no pivoting does. Scaled partial pivoting and complete pivoting are unmoved.forward error, relative to a solution of exactly (1, 1)no pivoting · as given1.3·10⁻⁴ 0 interchangesno pivoting · rows scaled1.2·10⁻⁴ 0 interchangespartial · as given2.2·10⁻¹⁶ 1 interchangepartial · rows scaled1.2·10⁻⁴ 0 interchangesscaled partial · as given2.2·10⁻¹⁶ 1 interchangescaled partial · rows scaled2.2·10⁻¹⁶ 1 interchangecomplete · as given1.1·10⁻¹⁶ 1 interchangecomplete · rows scaled2.2·10⁻¹⁶ 1 interchangethe same problem twicepartial, as given2.2·10⁻¹⁶partial, rows scaled1.2·10⁻⁴its relative residual1.2·10⁻¹⁶complete, rows scaled2.2·10⁻¹⁶the two systems have the same solutionand one pivot rule cannot see it
Fig. 17 And a decision that looks like an algorithm and is a scaling.
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. 18 Two residuals for one iteration, in the sequence field.
Distance from the truth and residual, for both methods, with 0% of the noise in the matrixFour bars, medians over 40 seeds. The upper pair is how far each answer is from the coefficients the problem was built from; the lower pair is ‖Ax − b‖ on the problem as given. Ordinary least squares minimises the lower quantity by definition, so its bar is the shorter of the two whatever happens above — and at this share it is the less accurate answer.the upper pair is distance from the truth; the lower pair is ‖Ax − b‖least squares · error0.04439total least squares · error0.1565least squares · ‖Ax − b‖4.734total least squares · ‖Ax − b‖4.971two orderingserror ratio (ls ÷ tls)0.28residual ratio (tls ÷ ls)1.1seeds40no vector makes the residual smallernot even the one the problem was built from
Fig. 19 A residual read the wrong way round.
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. 20 And one that cannot see what it is asked about.
Every column's 1-norm on a 12×12 matrix built to stop Hager's walk one column short12 bars, one per column, each the sum of the absolute values in that column — the quantity the estimator maximises, whose largest is the matrix's 1-norm. The walk visits 1 of them and returns 12.00; column 2 has 1-norm 114.00 and is never visited, because its entries alternate in sign and cancel against every sign vector the walk stands on.the estimator maximises this quantity over the columns it visitscolumn 1 ‹visited›12column 2 ‹the answer›114column 311.4column 411.4column 511.4column 611.4column 711.4column 811.4column 911.4column 1011.4column 1111.4column 1211.4estimate 12.0a walk that stopped earlythe estimate returned12the true 1-norm114columns visited1products with the matrix5the walk's own stopping test firedand every column it could see was smaller
Fig. 21 An estimator for the first author, which exists.
What Hager's estimator returns, as a share of the truth, against the size of the matrix built to defeat itThe estimate over the true 1-norm against n, both axes logarithmic. It falls from 0.526 at n = 4 to 0.0658 at n = 32, along the line 1/t where t is the construction's own multiplier — which is bounded only by the size of the matrix, so the ratio has no floor. A line at one marks a correct estimate.10¹10⁻²10⁻¹1size of the matrixestimate ÷ true 1-norma correct estimatewhat it returns1 / tno floorratio at n = 40.53ratio at n = 320.066products, either size5the estimate is always a lower boundwhich is the direction that flatters the matrix
Fig. 22 And the case built to fool it.
Iterative refinement from a 11-bit factorisation, κ = 10⁴A semi-logarithmic plot of forward error against refinement step. One curve falls steeply to the level of a double-precision solve; the other is nearly flat.012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹refinement step‖x − x*‖ / ‖x*‖a full double-precision solveresidual in11-bitresidual indoubleone argument apartκ·u of the factorisation4.9double residual, final0.0049same-precision, final1.430×30, κ = 10⁴, same factors in both runsidentical cost
Fig. 23 Buying back what the second author cost.
Two perturbation bounds and the error that was measured, on a 12×12 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 3.7·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)3.1·10⁸cond(A, x)8.5ratio of the bounds3.7·10⁷both curves above the data are boundsand only one of them is a measurement
Fig. 24 A backward error measured entry by entry rather than in a norm.
The error a reduced model has, and the bound computed before it existed‖H − Hᵣ‖∞ for balanced truncation of a 20-state model, against 2Σ_{k>r}σₖ — a number available from two Lyapunov solves before any reduced model is formed — and against σᵣ₊₁, which is a lower bound no reduced model of that order can beat. The measured error and the upper bound are the same curve: the ratio is 1.0000 to 1.0000 across the sweep, spread 1.0000. The bound is attained rather than approached, because on a spectrum falling this fast the tail sum is its own first term and removing one state costs exactly twice it. The measured error sits at 2.03 times the lower bound, so the whole question of what order to truncate at is answered by a curve nobody has to compute a model to see.123456710⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹order r kepterror and bounds2Σσ, and the errorσᵣ₊₁a bound that is an equalityorders7bound ÷ error, worst1spread over the sweep1error ÷ σᵣ₊₁, worst2.1computed before the modeland attained by it
Fig. 25 A field where the error is bounded before the computation instead.
Why a Gramian can be truncated at all: a rational approximation problemλₖ₊₁(P)/λ₁(P) for the controllability Gramian of a 30-state model, against Z_k², where Z_k is the smallest a degree-(k,k) rational function with poles on one side can be made on the other. The Gramian's spectrum falls from 1 to 2.48·10⁻⁸ in 10 steps, and the bound falls with it. Nothing about the model enters the bound except the two ends of the spectrum of −A, 9.86 and 3834, whose ratio is 388.8 — so the decay that makes model reduction possible is predictable from two numbers before the Gramian exists. That is the fact this whole field rests on and it is a fact about rational approximation.1234567891010⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹kλₖ₊₁ ÷ λ₁, and the boundZ_k²the Gramianpredicted from two numbersstates30κ of the spectrum389λ11 ÷ λ₁2.5·10⁻⁸the bound there5.2·10⁻⁴the cliff everything rests onand the reason for it
Fig. 26 And what makes that bound possible.
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 Half the digits lost for a reason that is nobody’s fault at all.
The condition number of one eigenvalue, at an off-diagonal entry of 1A rising curve against the size of an off-diagonal entry, on logarithmic axes, with a flat line at one below it and a vertical marker at the current value.110¹10²10³110¹10²10³off-diagonal entry ccondition number of the eigenvalue√(1 + c²)decoupled: 1measuredthree routes, one number‖A − ZTZᵀ‖/‖A‖1.4·10⁻¹⁵closed form1.4computed 1/|yᵀx|1.4worst measured movement0.65four eigenvalues, two conditioning numbersthe symmetric case has one, and it is 1
Fig. 28 The first author, measured for one eigenvalue.
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. 29 The hero again, beside the two terms this essay adds to it.

Where this goes next

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. 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 Six routes to one spectrum Three linearisations of one quadratic, each reduced to a standard eigenvalue problem two ways. All six have exactly the same eigenvalues in exact arithmetic. On a well-scaled problem they differ by noise; on a badly scaled one by a factor of forty; and two of the six are the same matrix. 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 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.

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 errorLinearisationResidualScalingUnit roundoff