Concept

Backward error — where it appears

The size of the smallest change to the problem for which the computed answer is exact. It is what an algorithm is answerable for, it is measurable without the true answer, and it is the smaller half of what a reader eventually sees.

Named by 68 essays across 17 fields — each of them below, with the objects they name alongside it.

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.

error · Backward error
nothing8.3%the right-hand side16.5%the matrix, slowly85.9%everything100.0%what can be reused: the answerwhat can be reused: the factorisationwhat can be reused: the factorisation, for a whilewhat can be reused: nothingshare of the cost of a sequence that shares nothingwhat 12 members cost in factorisationssame: factorisations1rhs: factorisations1drift: factorisations4independent: factorisations12what changes between the membersdecides what may be carried

The problem that arrives again

A hundred and thirty essays have solved a system once and measured how wrong the answer was. Almost no computation is shaped like that. A solve is one step of an outer loop, its answer is an input rather than a deliverable, and four quantities treated here as accuracy requirements turn out to be assets with a shelf life.

sequence · Sequence of solves
21-13-3-121-212-443-12as givenrows in the order 1 2 3 443-1201.251.252.502.51.5-30-0.5-0.52after step 1pivot 443-1202.51.5-3000.5400-0.21.4after step 2pivot 2.543-1202.51.5-3000.540003after step 3pivot 0.5‖PA − LU‖/‖A‖0largest multiplier0.75row order 4 3 2 1the pivot is chosen

Elimination is a sequence of choices

Gaussian elimination is taught as a procedure with no decisions in it. There is one decision at every step — which row to use — and every stability property the algorithm has comes from making it well.

elimination · Elimination
3456789101112110¹10²10³10⁴10⁵10⁶nwidest intermediate, in bitsrationals, not reducedfraction-free · rationals reduced · the answerone answer, three widthsn12answer40Hadamard bound52fraction-free40reduced rationals39unreduced1.4·10⁶the error is zero on every curvethe cost is the length of the numbers

An answer with no error in it

An integer matrix eliminated over the rationals rounds nothing, so the forward error is zero, the residual is the zero vector, and the identity this site is built on has no terms left. The cost does not vanish with the error. It moves into the length of the numbers, where three correct routes differ by four orders of magnitude.

exact · Exact cost
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

A small residual is not a small error

Substituting the answer back and finding that it fits is the most natural check there is, and it verifies the wrong thing. A residual of 10⁻¹⁷ is entirely compatible with an answer whose second digit is wrong.

error · Backward error
[ ε 1 ; 1 1 ] x = [ 1 ; 2 ], exact answer (1.000000, 1.000000)with partial pivoting1101U after elimination1.0000001.000000computed xbackward error 0forward error 0without10⁻¹⁷10-1·10¹⁷U after elimination0.0000001.000000computed xbackward error 0.25forward error 0.71no error is raisedgrowth 10¹⁷

The swap that is not optional

Run elimination without a row interchange on a matrix that needs one and nothing announces a failure. There is no division by zero, no warning, and an answer of the right shape. It is simply wrong, and how wrong depends on a number you did not look at.

elimination · Pivoting
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.

polynomial · Linearisation backward error
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

The scaling that buys ten orders

Two lines computed from three norms, a change of variable that is exact in both directions, and the whole of the loss the previous essay measured comes back — flat, at every stop, because after scaling every stop is the same problem.

polynomial · Polynomial scaling
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

The bound that is never attained

Partial pivoting's stability guarantee permits the entries to double at every step — a factor of 5.5·10¹¹ at n = 40. The measured growth on random matrices of that size is about three. The gap is eleven orders of magnitude, and the guarantee is still worth having.

elimination · Growth
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

An answer that is known

Almost every demonstration of numerical error estimates the error by computing the same thing more carefully. The Hilbert matrix does not need that: its inverse is a closed form in integers, so the true answer is available exactly and the error is measured rather than approximated.

error · Exact ground truth
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.

polynomial · Linearisation backward error
012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹refinement step‖x − x*‖ / ‖x*‖a full double-precision solveresidual in24-bitresidual indoubleone argument apartκ·u of the factorisation6·10⁻⁴double residual, final3.2·10⁻¹³same-precision, final1.3·10⁻⁴30×30, κ = 10⁴, same factors in both runsidentical cost

Buying the accuracy back

Factorise in single precision, then correct the answer using residuals computed in double, and the result is what a full double-precision solve would have given. Compute those residuals in single instead and the identical algorithm, at identical cost, recovers nothing.

arithmetic · Mixed-precision
forward error, relative to a solution of exactly (1, 1)no pivoting · as given1 0 interchangesno pivoting · rows scaled1 0 interchangespartial · as given0 1 interchangepartial · rows scaled1 0 interchangesscaled partial · as given0 1 interchangescaled partial · rows scaled0 1 interchangecomplete · as given0 1 interchangecomplete · rows scaled0 1 interchangethe same problem twicepartial, as given10⁻¹⁸partial, rows scaled1its relative residual10⁻¹⁷complete, rows scaled10⁻¹⁸the two systems have the same solutionand one pivot rule cannot see it

The pivot that reads the units

Partial pivoting compares the entries of a column and takes the largest. Those entries carry units, so the comparison depends on them — and there is a row scaling, on the standard two-by-two that pivoting exists to fix, which makes partial pivoting perform the identical catastrophic elimination it was introduced to prevent, with no interchange at all.

elimination · Pivoting
-14-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μ — the barrier parametercondition number, and relative errorκ₂, condensedκ₂, augmentederror, condensederror, augmentedagainst a BigInt answerκ₂ augmented, μ = 10⁻¹⁴3·10¹⁵its relative error10⁻¹⁵κ₂ condensed2.4·10¹⁶its relative error0.31the same step, written two waysand only one of them is solvable

A condition number sent to infinity

An interior-point method manufactures an ill-conditioned matrix on every iteration, deliberately, because the separating of a diagonal is how it discovers which constraints are active. Written one way the answer keeps fifteen digits at a condition number of 3·10¹⁵. Written the other way — the way almost every code writes it — it has none left.

constraint · Interior-point conditioning
10⁻³10⁻²10⁻¹1110¹10²10³10⁴pivot threshold τgrowth factor · entries in L+U, ÷ entries in Agrowthfillthe library default‖PA − LU‖/‖A‖ at τ = 0.12.9·10⁻¹⁶growth at τ = 0.138entries in L+U at τ = 0.1372one knob, two measurements, opposite directionsand the default is most of both

A threshold between fill and growth

One number decides how small a pivot an elimination will accept. At 0.001 the factor holds 172 entries and the matrix grows by 1,330; at 1 it holds 260 and grows by 1.2. The libraries ship 0.1, and the measurement says why.

sparsity · Sparse pivoting
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

A condition number scaling cannot move

Skeel's componentwise condition number is invariant under any row scaling — exactly, before any norm is taken, because two diagonal factors cancel entry by entry. It is never larger than the normwise one and can be arbitrarily smaller, and the ratio between them is a diagnostic for which kind of ill-conditioning a matrix has.

error · Scaling
10⁻²10⁻¹110¹10²10³10⁴10⁻⁴10⁻³10⁻²10⁻¹1backward error, in units of ushare of systems above ituCramereliminationCramer, control24-bit arithmeticworst Cramer, in u395worst elimination, in u1.2control, worst Cramer4.9κ of the worst system3.4·10⁶one derivation, two computationsand only one of them is stable

A rule that is correct and unusable

Cramer's rule gives every component of the solution in closed form, in terms of determinants, and it is a theorem. On two-by-two systems whose rows are nearly parallel it returns an answer with a backward error of 458 units of roundoff where elimination returns 1.3 — on a matrix whose condition number is 32,000 and which elimination solved perfectly.

elimination · Determinant
10²10³10⁴10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²steps taken‖RᵀR − AᵀA‖ ⁄ ‖AᵀA‖the bound, linear in the steps√k · uevery step safe, the chain notdrift after the run3.9·10⁻¹⁴the bound there3.3·10⁻¹³√k · u there6.1·10⁻¹⁵worst single amplification2.7worst leverage met0.69refreshes1backward stable onceand three thousand times is a different claim

Stable once, and three thousand times

A sliding window adds a row and removes one at every step and never looks at the data again. No single step of it amplifies by more than 2.72, no downdate fails, and after three thousand steps the triangular factor in memory is 3.9·10⁻¹⁴ from the matrix it is supposed to be a factor of — six hundred times growth from a per-step bound that says nothing about chains.

sequence · Sequence stability
02468101210⁻¹⁸10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸member of the sequencebackward error of the solvefresh order · equilibrated kept orderringed: another pivot replaced hereone order, two sets of unitsfresh order, last member9.5·10⁻¹⁷kept order, last member4.8·10⁻⁹equilibrated, last member6.5·10⁻¹⁷pivots replaced, kept2pivots replaced, equilibrated0the perturbation floor √u1.1·10⁻⁸a reused pivot order is safe or notdepending on what the rows are measured in

The order that was right last time

A pivot order computed once and reused across a sequence saves the symbolic phase, and the price is that a pivot which was large may now be small. Replacing it with √u·‖A‖ costs eight orders of backward error and iterative refinement recovers a factor of 8.8 of them. Divide each row by its largest entry first and the same reuse costs nothing at all.

sparsity · Sparse pivoting
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.

structure · Structured backward error
567891010⁴10⁵10⁶10⁷10⁸10⁹log₂ nmultiplicationsdense factorisation, n³⁄3the recursion, countedwhere the format starts payingratio at n = 641.5ratio at n = 5120.16exponent, first doubling2.1exponent, last doubling1.7backward error1.4·10⁻¹⁰cheaper is a sizenot a property

Where the format starts paying

A hierarchical solve costs 1.48 times a dense factorisation at 64 unknowns and 0.16 times it at 512. The crossover is between 64 and 128, it walks right when the accuracy is tightened, and the exponent between consecutive sizes is 2.13, 1.93, 1.74 — falling towards one and never arriving.

cost · Hierarchical solve
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.

polynomial · Approximation before linearisation
-9-7-5-3-1110⁵10⁶10⁷log₁₀ ε of the preconditionermultiplications for the whole solveno preconditioner: 41 stepsiteration count, rescaledmultiplications spenttwo curves, opposite waysκ21steps, no preconditioner41cheapest ε0.5its rank1tightest ⁄ cheapest2.4the count is what is printedand the cost is what is spent

The accuracy worth paying for

Used as a preconditioner, a hierarchical representation gets better at every accuracy — the iteration count falls monotonically all the way to the tightest tolerance. The total work does not. Its minimum sits at a rank-one preconditioner on an easy problem and six decades further along on a hard one.

hierarchy · Hierarchical solve
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² entries7.9·10⁴as one number, ρ62backward error, unconstrained2·10⁻¹⁷as a symmetric Toeplitz matrix5.6·10⁻¹²fewer numbers, better conditionedand no nearby problem left

The condition number of the model

Describe a 40×40 Toeplitz matrix by its 1,600 entries and its condition number is 78,800. Describe it by the one number it actually contains and the condition number is 61.9. The three decades in between are not an approximation or a bound — they are what κ has been over-stating, and the drop is not where the linear algebra is.

structure · Structured backward error

The inverse that is never formed

x = A⁻¹b is how the solution of a linear system is written and it is not how it is computed. The usual reason given is cost — three times the arithmetic. The real reason is that one of the two routes is backward stable and the other is not, and at κ = 10¹⁴ they differ by twelve orders of magnitude in the number that says whose fault a wrong answer is.

elimination · Inversion

The zero you are allowed to write

A deflation criterion sets a subdiagonal entry to zero because it is small. A drop tolerance discards an entry of a factor because it is small. A truncation discards a singular value because it is small. Three fields, three vocabularies, no shared arithmetic — and plotted as work saved against error accepted, one curve.

error · Deliberate zero

A perturbation that keeps the symmetry

The smallest perturbation that makes a computed answer exact is the backward error. Ask for the smallest one that also keeps the problem's structure and the number can only go up — and measured on a palindromic quadratic it goes up by 1.17, while the structure the computed spectrum has lost is not in either number.

structure · Structured backward error

The gap refinement can close

Multiplying by a computed inverse is not backward stable, and refinement at the working precision repairs it. That much is settled. The claim beside it — that the forward error does not move — was read at one conditioning and four corrections too late. Swept over ten, it moves at every one, and it lands on the LU route's own number after a single correction.

elimination · Inversion

An exact answer to a measured problem

The residual is the zero vector, nothing was rounded at any step, and the answer is wrong in its first digit. Data accurate to fourteen places, an exact solve of the system it defines, and an error of 10⁻⁵ — because conditioning was never a statement about arithmetic and removing the arithmetic error removes none of it.

exact · Exact cost

A knob calibrated in residuals

A formatted Cholesky has two numbers in it and only one of them is an accuracy. Across twelve trees — three sizes by four leaf sizes — the leaf moves the truncation count from 0 to 258 and moves the ranks of the blocks not at all, while the residual follows the tolerance at slopes between 1.022 and 1.046 and sits at about a tenth of it throughout.

hierarchy · Recompression

The repair the drift did not need

A sliding window's carried Cholesky factor drifts 3.9·10⁻¹⁴ from its data, and multiplying by κ(AᵀA) predicts eight lost digits in the coefficients, a stream conditioned at 10¹² losing the answer, and a periodic refresh of the factor as the default repair. Measured against coefficients computed exactly in rationals, all three come out differently. On a stream made ill-conditioned by scaling, the conditioning never reaches the coefficients. On a collinear stream, a freshly recomputed factor is as wrong as the drifted one. And one correction from the window's own rows reaches Householder's accuracy for a fraction of a refresh's cost.

sequence · Sequence stability

Deciding that a zero has arrived

The previous tolerances were offers — accept this much error, save this much work. A detection threshold is not an offer, because both directions are failures. One matrix here has three genuinely near-invariant subspaces, and the constant somebody typed decides which of them the recurrence stops at; at eight significand bits the same kind of constant produces a proof of something false.

error · Deliberate zero

A ranking whose order is not determined

Of the fourteen adjacent comparisons in a top-fifteen, fourteen survive perturbing the arithmetic at the rounding level, eleven survive moving the teleportation parameter across its usual range, and three survive removing one link. The computation is the strongest part of the answer.

graph · Perron frobenius

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.

error · Backward error

A worst case is as fragile as its margin

Wilkinson's matrix grows by 5.5·10¹¹ under partial pivoting, and adding Gaussian noise of 10⁻¹⁴ to every entry takes its median growth to exactly 2. The shooting matrix from a boundary-value problem grows by 1.1·10⁴, and noise ten million times larger leaves it untouched. The difference is what each worst case rests on. Wilkinson's rests on exact ties between candidate pivots, which any noise breaks. The shooting matrix's rests on a choice made by a margin of 5.6·10⁻³, and between 10⁻⁶ and 10⁻³ its median growth is that margin divided by the noise, times a constant between one half and four thirds.

elimination · Growth

Bracketing an error nobody can measure

The error of a low-rank Gramian factor is the one quantity a caller cannot compute, because computing it needs the Gramian the factor exists to avoid forming. Two numbers that can be computed sit either side of it — a rational factor known before the run, and a residual known after — and they stay a factor of four apart across a fourfold change of size.

reduction · Gramian decay

A tensor that cannot be decomposed

Every member of a certain sequence is exactly a sum of two rank-one terms, and both terms are written down in closed form. A three-hundred-sweep fit from a random start does not find them, and stalls at the same one per cent however far the sequence goes — while a fit started at the answer loses digits exactly as 2n² says it should.

error · Conditioning

The knob that moved two things

Decide how many digits the answer needs, divide by the condition number, and compress to that. It is the one rule licensed in advance here, and its two factors are not the independent inputs it reads as: the partition's leaf moves neither of them and moves the answer by nearly a factor of three, and the only knob here that raises κ halves the ranks while it does so.

hierarchy · Hierarchical solve

The growth a boundary-value problem supplies

Large growth under partial pivoting is usually said to need a matrix built for it. A two-point boundary-value problem solved by multiple shooting supplies one without being asked: its growth factor is e^(5T/6)/2 to four figures — 1.1·10⁴ at an interval of 12, 1.3·10⁵ at 15 — on a matrix whose condition number never exceeds 8.3, while rook and complete pivoting keep it below 2. And it is the finer shooting grid that grows: below a step of 0.3397 the choice partial pivoting makes turns on one entry against one, and above it there is no growth at all.

elimination · Growth

The offset that moved the slope

The accuracy is supposed to lift a storage curve and leave its growth alone. Measured at seven tolerances, the strong partition adds 10.9 numbers per unknown per doubling at two digits and 36.5 at twelve — the growth rate more than triples, so ten decades of accuracy cost 33 per cent more storage at 64 unknowns and 118 per cent at 512.

hierarchy · Storage growth

The number that cannot rank them

Levinson and Gaussian elimination are indistinguishable on the backward error a library reports — every one of ninety-six measurements between 1.16·10⁻¹⁷ and 5.73·10⁻¹⁷. The structured backward error separates them by up to a hundredfold, in whichever direction the point happens to give. Only the forward error ranks them, and only because this family's exact answer is known.

structure · Structured backward error

Two condition numbers of one matrix

κ₂ is a worst case over perturbations of a given norm, and a normwise perturbation may put its whole budget on the smallest entry. The componentwise number is a worst case over perturbations proportional to the entries, which is what a backward-stable factorisation actually makes. On one matrix they are 3·10¹³ and 13.3, and the error obeys the second.

error · Scaling

Which of the choices is doing the work

Elimination makes n − 1 decisions and they are not worth the same. On 8×8 standard normal matrices, removing the first pivot search and leaving the other six multiplies the median growth factor by 1.624; removing the last multiplies it by 1.000. The cost falls monotonically along the run, and the worst single matrix in the sweep grows by 2,366 when one early decision goes — so the median is the wrong statistic and the tail is where pivoting earns its reputation.

elimination · Elimination

One number that has to be right

Householder's orthogonality was called structural: a reflection is built from a unit vector, so rounding the vector names a different reflection rather than a broken one. Tested by breaking it, the claim is narrower and sharper. Perturb every component of the reflector by a relative 10⁻², and ‖QᵀQ − I‖ stays at 1.5·10⁻¹⁵ while the factorisation moves to 5·10⁻³. Perturb the one stored scalar by the same amount and ‖QᵀQ − I‖ is 6.5·10⁻². The structure is one degree of freedom, and the departure is four times its relative error.

orthogonality · Householder

The roots are not the coefficients

A polynomial whose roots are the integers one to twenty, expanded exactly, handed to the routine every library uses. The computed roots are wrong in the third digit, the computation is backward stable for the matrix it factorised, and above degree eighteen the coefficients are not double-precision numbers at all.

error · Conditioning

The order the greedy rule cannot choose

Wilkinson's matrix is the standard demonstration that partial pivoting's growth bound of 2^(n−1) is attained. It is attained by the row order the greedy rule picks, and not by the matrix: a single cyclic shift of the rows gives growth 2 at every size, with no multiplier above one. At n = 7 that is 64 against 2. And perturbing one entry by 10⁻¹² leaves the good order exactly where it was while putting every tie-break of the greedy rule back on 64.

elimination · Elimination

A triangle where the scalar was

Every level-3 QR assembles a block of reflectors into Q = I − Y T Yᵀ, and T is computed by a recurrence whose inputs are its own previous columns. A block of sixteen carries 136 computed numbers where sixteen separate reflections carry sixteen. The orthogonality it produces is 3.9·10⁻¹⁵ against the single reflector's 7.8·10⁻¹⁶ — a factor of five for a hundred and thirty-six times as many things that have to be right.

orthogonality · Householder

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.

error · Backward error

What a regression test can ask for

The machine's own variation on one solve is 3.2·10⁻¹², and the smallest defect whose answers clear it is one part in 10¹². The tolerance exists, it is bracketed on both sides by a factor of 1.42, and it is neither zero nor the 10⁻⁸ that usually gets typed.

machine · Regression tolerance

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.

polynomial · Linearisation backward error

Noise the growth amplifies

The shooting matrix's growth of 1.1·10⁴ fell under noise as its pivot margin divided by the noise, and the explanation offered was noise reversing partial pivoting's choices. But noise of 10⁻⁵ is five hundred times smaller than that margin. Add the same noise only to the entries that are not zero and every draw keeps the whole growth up to 10⁻³. The dense noise was not reversing the comparisons by itself: it sat in the zeros, the elimination multiplied it by the growth already made, and a comparison flips when that product reaches about four tenths of the margin — at every noise level from 10⁻⁶ to 10⁻³.

elimination · Growth

A count that comes out of a determinant

The number of spanning trees of a graph is the determinant of its grounded Laplacian, so it is a whole number known in advance. The elimination that computes it is backward stable at every size — and from sixteen vertices the answer is wrong, because the count has seventeen digits and a binary64 has sixteen.

graph · Graph elimination

The fifth author

Four authors of a wrong answer have been named on this site and each is a statement about one computation. The fifth is not: it is what separates two computations that are both correct, it is a backward error of measurable size, and no residual, bound or condition number contains it.

error · Backward error

A margin the factorisation records

Partial pivoting finds the second-largest candidate in every column it scans, and throws it away. Keep it, divide the gap by the growth reached at that step, and take the smallest over the steps before half the final growth has arrived. That one number, recorded by the factorisation that is already running, predicts the dense noise that halves the growth to within a factor of 0.92 to 2.4 on shooting matrices whose growth runs from 75 to 1.3·10⁵ and whose fragility spans three and a half decades — and on Wilkinson's matrix it is exactly zero.

elimination · Growth

The correction lost to its own two-by-two solve

Solving a band matrix through a circulant and a small correction was measured losing the answer to 10⁻⁴ where elimination kept it to 10⁻¹⁴, and a wrap that landed one sample on a zero was measured costing four orders more than one that landed a pair. Both measurements solved the two-by-two correction system by Cramer's rule. Solved with a row interchange, the corrected solve on the same matrix loses 2.9·10⁻¹⁰ — the cancellation, and nothing multiplied onto it — and the single landing costs what the pair costs. The loss was in the determinant.

structure · Circulant

A tolerance is priced by the problem

Three tolerances from three fields sit on one pair of axes and agree to within a factor of 5.74. That factor is the ratio of the two curves that cannot move. Change the only problem in the comparison and the third curve's fitted slope swings from 0.188 to 0.040 while the printed spread does not shift by a digit.

error · Deliberate zero

Three walks and one bound

A left-to-right sum, a chain of three thousand rotations and a conjugate gradient residual recurrence share no arithmetic and no vocabulary. Each has a standard bound that is linear in whatever it accumulates against. All three come out at a half — 0.486, 0.554 and 0.507 — and nothing is rescaled.

arithmetic · Summation

Two precisions guard the other edge

Run a lattice relation search at N digits and again at N/10, and accept its answer only if both runs return the same vector. Among six measured numbers, where the gap between the shortest and next vector vouches for one found relation in 35, the two runs agree on 34. Past a double's sixteen digits they never once agree on an accident, 0 of 561. They do agree on 62 accidents at fifteen digits or fewer — approximate relations that really are the shortest vector there — and the gap, which cannot see a relation among six numbers, can see those.

exact · Lattice reduction

The answer the last window left

A sliding window that corrects its least-squares answer at every step could start each correction from the previous step's corrected answer instead of from a fresh solve: the two windows share all but one row. On a stream with any noise in it, that start is three orders worse. The window's exact answer moves by 0.79 of itself in one step at κ(A) = 3·10⁶, a fresh seminormal solve is wrong by only 1.8·10⁻⁴, and one correction contracts either start by the same factor — so the fresh start ends at 2.2·10⁻⁸ and the carried one at 3.7·10⁻⁵. Only on data that agree exactly does carrying win.

sequence · Sequence stability

A threshold that holds the growth still

Partial pivoting's large growth on the shooting matrix rests on a margin of 5.6·10⁻³, and Wilkinson's rests on exact ties, which a perturbation at rounding breaks. A sparse code pivots with a threshold instead, taking the sparsest row among candidates within a factor τ of the largest. At every threshold from 1 to 0.1 both matrices grow by exactly as much as under partial pivoting. But the growth is now held by row counts, which no perturbation of the values can reverse: at τ = 0.1 it survives noise on the stored entries up to 0.56 on the shooting matrix, three hundred times more, and 0.18 on Wilkinson's, where partial pivoting's goes at 10⁻¹⁶. The margin that predicts it is the chosen candidate's height above the threshold line, and it needs no division by the growth.

elimination · Growth

A zero no twist can step around

Solved through its circulant wrap with the small capacitance system factorised stably, a banded matrix was found to lose only what the cancellation in the correction costs, and a twist of half a step kept that cancellation near one. That holds only while the cancellation is large. With it kept small, the corrected solve's error is a quarter of the wrap's condition number times u on every case measured — and on a symbol with a zero of fourth order, the square of the Laplacian, no twist can keep the wrap well conditioned, because the nearest sample any twist can reach sees (π/n)⁴. At n = 64 the corrected solve is 34 times less accurate than elimination, where the Laplacian's is four.

structure · Circulant

The number that moves when the problem does

Two quantities are offered as the condition number of one eigenvalue. One is unmoved to eight digits by a change of variable that is exact in both directions, and grows like the square root of the chain's length. The other is inflated by ten orders by that change of variable, and is ten times too large before anything has been done at all.

polynomial · Linearisation backward error

The licence is not the boundary

Cholesky QR is licensed by κ²u ≪ 1, which reaches equality at κ = 9.5·10⁷ in double precision. At 10⁸ the factor it returns is already 0.37 away from orthogonal, and it goes on returning factors as far as 10¹³ — refusing at scattered condition numbers in between, at different ones for eight columns and for six.

machine · Algorithm selection

Two matrices and one problem

Ax = λBx is what a finite element model, a structural vibration and a constrained optimisation actually produce, and it is not the one-matrix problem with a change of variables. Everybody is told not to form B⁻¹A because it is not symmetric. That is true, the departure from symmetry is about one, and it is not what decides the accuracy.

spectra · Pencil

A problem with no answer

If two matrices share a null vector then det(A − λB) is identically zero and every λ is an eigenvalue, which means none of them is. Perturb such a pencil by a ten-billionth and a solver returns six numbers with residuals below 10⁻⁹. Change the seed and it returns six different numbers, spread over forty-four, with residuals just as small.

spectra · Pencil

The residual the method reports

Conjugate gradients prints a relative residual of 6.9·10⁻²¹. The unit roundoff is 1.1·10⁻¹⁶, so that is not a small residual and not a large one — it is not a residual. The vector the method is holding at that step has ‖b − Ax‖/‖b‖ = 5.1·10⁻¹⁰, and nothing in the run says so.

iterative · Residual gap

The number that is re-derived

GMRES prints a residual it never computes from its answer either. On the matrix that sends a conjugate gradient recurrence 7.3·10¹⁰ wrong, and on two others chosen to be worse, its number is never more than a factor of 2.86 out — while the basis it is computed from has lost orthogonality entirely. The disease is not iterative methods, and it is not floating point.

iterative · Residual gap

A different equation on every grid

Upwinding is the exact discretisation of a convection–diffusion problem with diffusion ε + h/2, entry for entry, at a relative difference of between 0 and 1.26·10⁻¹⁶ on every mesh from 15 points to 511. The equation it is exact for is chosen by the mesh and not by ε — the added diffusion is 0.01563 on a 31-point grid whether ε is 0.2 or 0.001.

iterative · Convection

Named alongside it

The objects these essays reach for when they reach for this one.

Condition numberExact ground truthForward errorResidualGaussian eliminationGrowth factorPartial pivotingIterative refinementToleranceLinearisationUnit roundoffFlop count

All concepts