Numerical Linear Algebra
where the algebra stops being the arithmetic
Concept

Forward error — where it appears

The distance from a computed answer to the true one, which needs the true one and is therefore known only on constructed problems. It is the number a reader actually wants, and it equals the backward error times the condition number, which is why both factors are printed here.

Named by 36 essays across 16 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 01234510⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹log₁₀ κ(A)condition number, and relative errorκ(S)κ(ZᵀHZ)range-space errornull-space erroragainst a BigInt answerκ(S) at κ(A) = 10⁵4·10¹⁰κ(ZᵀHZ), all stops21range-space forward error5.3·10⁻⁶null-space forward error5.8·10⁻¹²both are the same algebraand only one squares

Two ways to remove a constraint

A constrained system can be reduced by eliminating the multipliers or by eliminating the constrained directions. Both give the same answer in exact arithmetic and inherit different condition numbers — one of them squares the constraint's, and the other does not contain it at all.

constraint · Saddle-point systems 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 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 21018263442500102030405060708090100110120130140150160170180190200210220230240250260270280290300310320330340significand bitsiterationspreconditioner roundedarithmetic roundedno preconditionersame bits, different casualtyerror, 3-bit preconditioner8.8·10⁻¹³error, 3-bit arithmetic0.16‖A − LLᵀ‖/‖A‖ of the factor0.083a direction may be roundeda measurement may not

The part of a solver that may be rounded

A preconditioner computed and applied with a three-bit significand still returns thirteen correct digits — it costs seventeen extra iterations and nothing else. Round the working arithmetic instead and the step count barely moves while the answer loses exactly the digits the format dropped.

combination · Mixed-precision 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 00.250.50.75110⁻¹110¹share of the noise placed in the matrixleast-squares error ÷ total least-squares errorequally accuratetotal leastsquares aheadordinary leastsquares aheadthe model, not the methodadvantage, all noise in b0.28advantage, all noise in A2.5seeds at each share40the same total noise at every pointand only where it sits changes

When the matrix is wrong too

Every least-squares problem here has assumed A is exact and b is not, and moved b onto the column space of A. Where both were measured, the smallest correction that makes the system consistent moves the matrix as well — and on the problems where that answer is more accurate, it has the larger residual, by construction rather than by luck.

leastsquares · Total least-squares 10¹10³10⁵10⁷10⁹10¹¹10¹³10¹⁵10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A), the matrix that was updatedrelative forward errorSherman–Morrisondirect solve of A + uvᵀone answer, two routesκ of the answer's matrix1κ of the matrix replaced10·10¹³update formula's error2.5·10⁻⁴direct solve's error1.1·10⁻¹⁶the question's condition number is 1at every point on this axis

A correction cheaper than the problem

Sherman and Morrison's formula updates a solved system for a rank-one change to the matrix, at 4n² operations instead of (2/3)n³. It is exact algebra. On a problem whose updated matrix is the identity — condition number one, the easiest system there is — it returns a forward error of 2.5·10⁻⁴ where a direct solve returns 10⁻¹⁶.

leastsquares · Low-rank update -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 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

The units the matrix is measured in

One linear system, written twice. The rows of the second are the rows of the first in different units, the solution is identical to the last bit, and the condition number has moved by eight orders of magnitude. One of those two numbers is a fact about the problem and the other is a fact about the notation.

error · Scaling 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 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 10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A)relative errorforward, A⁻¹bforward, A\bbackward, A⁻¹bbackward, A\bat κ = 10¹⁴η, LU solve2.2·10⁻¹⁷η, via the inverse4.5·10⁻⁵forward, LU solve2.8·10⁻⁴forward, via the inverse0.015one factorisation, two ways to use itand one of them forfeits the 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 01234510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1correction steprelative errorLU route: η = 2.2·10⁻¹⁷LU route: forward 2.8·10⁻⁴forward errorbackward errorwhat a correction buysη before refinement4.5·10⁻⁵η after four steps2.8·10⁻¹⁷forward, unchanged7.1·10⁻⁴cost of a step, flops1800the residual is repairableand the accuracy floor is the problem's

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 0246810121410⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹interior-point iterationrelative error, and μcertified from iterate 1μthe iterate's errorthe crossover's errorone solve, checkediterations15first certified iterate1iterate error there0.22crossover error there4·10⁻¹⁴the active set arrives long before the digitsand a crossover collects them at once

The active set before the digits

An interior-point method takes fifteen iterations on a quadratic programme with forty constraints, and its iterate has eight correct digits at the eleventh. Take the constraints its diagonal calls active at the first iterate, solve the equality problem they define once, and check the answer against the conditions for optimality. It passes, to thirteen digits. The step's matrix had a condition number of 43 at that iterate, and 7·10¹⁵ at the last.

constraint · Interior-point conditioning the upper pair is distance from the truth; the lower pair is ‖Ax − b‖least squares · error0.08523total least squares · error0.037least squares · ‖Ax − b‖3.965total least squares · ‖Ax − b‖4.078two orderingserror ratio (ls ÷ tls)2.3residual ratio (tls ÷ ls)1seeds40no vector makes the residual smallernot even the one the problem was built from

The two numbers a caller has

Choosing between the two least-squares methods is a statement about where the noise is, and the two quantities a caller can compute are both blind to it. The residual separates the answers by 0.14 per cent where their accuracies differ by 14, and κ(A) falls from 3.54 to 2.46 across a sweep in which the error rises by a factor of sixty-two.

leastsquares · Total least-squares 10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻⁷10⁻⁵10⁻³10⁻¹10¹relative error in the datarelative error in the answerthe answer's errorthe data's errornothing here is the arithmetic'sdata error10⁻¹²answer error10·10⁻⁴amplification10·10⁸rounding committed0residual0the residual is exactly zeroand the answer is wrong anyway

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 10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖x̂ − x‖ / ‖x‖modified, through Qᵀbclassical, either route‖QᵀQ − I‖, modifiedHouseholdermodified, on [A b]κ²uκuforward error at κ = 10⁸Householder4·10⁻⁹modified, through Qᵀb0.13modified, on [A b]2.7·10⁻¹⁰classical840×8, three seeds, against an exact rational solveQ is not orthogonal; x is right

The right-hand side as one more column

Modified Gram–Schmidt's Q is 4.3·10⁻⁹ from orthogonal at κ = 10⁸, and a least-squares solve that multiplies b by it is wrong by 0.13. Hand the same routine b as an extra column instead and the answer is right to 2.7·10⁻¹⁰ — closer than Householder's 4.0·10⁻⁹. Classical Gram–Schmidt gains nothing from the same trick, to the last bit.

orthogonality · Gram–Schmidt -1-0.500.511.520eigenvalue of P⁻¹K1 − φ1φS never formedγ1smallest ν0.056furthest from φ, 1 − φ0.57MINRES steps11an exact Schur approximation, bought by changing Hthe golden ratio without S

Where the augmentation puts the cost

Add γAᵀA to the objective block of a saddle-point system and its Schur complement tends to I/γ, so the cheapest possible approximation becomes the right one and the golden-ratio spectrum arrives — within 7.6·10⁻⁶ at γ = 10⁶. MINRES falls from 21 steps to 6. The inner solve with the augmented block rises from 14 conjugate gradient steps to 43, their product does not fall at all, and the answer loses seven and a half digits on the way.

constraint · Block preconditioning 10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³relative compression error, 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.

error · Backward error 234567891010⁻⁴10⁻²1ADI steps krelative error, and max|r|²max|r|²the measured errora bound that is an identitysteps drawn9ratio, worst0.77spread of the ratio1error at k = 103.9·10⁻⁴its residual1.1·10⁻⁴four orders of errorand one constant

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 110¹10²10⁻¹²10⁻⁸10⁻⁴110⁴ncondition number, and residual reachedmarks above: κ of the step · dashes: its closed formmiddle: a fit from a random startbelow: the same fit started at the answera decomposition that is ill-conditionedκ at n = 1283.3·10⁴its closed form3.3·10⁴cosine of the terms1from a random start0.012from the answer1.4·10⁻¹⁰the answer existsand cannot be found

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 10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹relative compression error, chosen‖x − x*‖ ⁄ ‖x*‖κ × the chosen backward errorthe error in the answerboth factors known firstκ21chosen at 10⁻⁸1.4·10⁻⁹predicted forward2.9·10⁻⁸measured forward1.6·10⁻⁹bound ⁄ measured26the amplifier, used forwardsfor once

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 10⁻³10⁻²10⁻¹110⁻¹⁹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹1 − ρbackward errorLevinson, Toeplitzelimination, ToeplitzLevinson, any kindelimination, any kindat ρ = 0.999Levinson, any perturbation2.9·10⁻¹⁷Levinson, Toeplitz only6·10⁻¹²the ratio between them2·10⁵diagonal defect of the first0.97the number that is reportedand the number that is asked about

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

One line that buys a quarter of the run

The adaptive forcing rule has a floor on it that no published statement of the rule carries: do not solve a step to an accuracy the outer loop will not use. Removing it costs 9 to 27 per cent of the whole inner run. Keeping it costs between 23 and 2,600 times the forward error — accuracy the residual test never asked for and both runs satisfy the test either way. The line is a trade between a residual and an error, and which of the two the caller meant decides whether it is a saving.

sequence · Inexact newton

A test with no tolerance in it

An interior-point method's own stopping test is a tolerance on μ, and at the tightest it can be set to it stops after 15 iterations with 1.6·10⁻¹². A crossover from the iterate at 1.5 returns a point whose error is 1.8·10⁻¹⁴ — ten times sooner and a hundred times better, from an iterate carrying one correct digit. One attempt costs an eighth of a step, and the guess's own margin says which iterate to spend it on.

constraint · Interior-point conditioning

The circulant the problem did not contain

A matrix that differs from a circulant in two corner entries can be solved through the circulant, by a transform and a two-by-two correction, and the cost claim is exact. The accuracy claim is not. On tridiag(−1, 2 + σ, −1), whose condition number stops at 1,712, the correction is wrong by 1.2·10⁻⁴ at σ = 10⁻⁸ while elimination is right to 1.1·10⁻¹⁴ — because the periodic neighbour is singular at σ = 0 and the two-by-two system inherits that. Solved by Cramer's rule, as here, the two amplifications multiply; a later measurement found that a pivoted solve of the same two-by-two system removes the second.

structure · Circulant

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

Two near-zeros cost less than one

Solve a well-conditioned tridiagonal matrix through a nearly singular wrap and the correction's accuracy is not set by how singular the wrap is. At κ = 4·10⁷ one wrap returns the answer to 8.8·10⁻¹¹ — better than κ·u — and another, at κ = 3.8·10⁷, returns it to 3.3·10⁻⁶. The difference is how many of its samples sit near the symbol's zeros. A real wrap lands on a conjugate pair, a rank-two correction absorbs the pair exactly, and its two-by-two system has condition number 1.00 — which mattered because that system was solved by Cramer's rule; solved with pivoting, the single landing costs what the pair costs.

structure · Circulant

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

One sequence and two recurrences

CGLS and LSQR compute the same iterates — the minimiser over a space is unique, so there is nothing to choose between them in the algebra. At κ = 10⁶ they cost 42 steps and 47. At κ = 10¹⁰ they cost 110 and 209, across four seeds, and the quantity that separates them is the orthogonality of a basis neither of them keeps.

iterative · Krylov

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

Nine steps of pessimism

A proved bound is 8 to 26 times the error it bounds, at every precision from 16 to 40 significand bits. A carried interval is (√2)ᵐ times too wide after m re-enclosures. The two cross between eight and nine, so the method everybody warns against is the tighter of the two for a short computation.

arithmetic · Interval

The reading that never moves

Across thirty runs — five grids from 36 to 196 unknowns, six working precisions from 53 significand bits down to 8 — the residual conjugate gradients stops on stays between 1.10·10⁻¹³ and 9.95·10⁻¹³. Over the same thirty runs the error of the answer spans a factor of 2.39·10¹¹, and the step count more than doubles. The one number the run publishes is the only one that responds to neither axis.

machine · Stopping test

The error the method already knows

Summing the exponential's Taylor series throws away a known number of digits, and the number is on the machine while the sum is being formed. The largest term divided by the answer, times the unit roundoff, tracks the relative error that comes out — to within a factor of nine, across fourteen orders of magnitude of it — and nothing reports it.

spectra · Matrix function

Named alongside it

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

Condition numberBackward errorResidualExact ground truthUnit roundoffHilbert matrixIterative refinementStopping criterionLinearisationSaddle-point systemsBackward stabilityComponentwise condition number

All concepts