Numerical Linear Algebra
where the algebra stops being the arithmetic

The algebra is exact. The arithmetic is not. Everything interesting happens in the gap, and the gap can be measured.

Solve a thirteen-by-thirteen Hilbert system whose answer is exactly the integers one to thirteen, and double-precision elimination returns 0.05 where 8 belongs. Nothing went wrong: the backward error of that solve is 10⁻¹⁷, which is to say the algorithm answered a neighbouring problem perfectly. The fault is the problem's, not the algorithm's, and those two things are separately measurable — which is what this site is about. The geometric half of linear algebra is explained superbly elsewhere. This is the other half: where the algebra stops being the arithmetic, and whose fault it is.

Backward and forward error against the condition number, on 8 × 8 systemsA log–log plot over twelve decades of condition number, at 8 × 8, twenty seeds a point. The backward error is flat — median 6.12·10⁻¹⁷ at κ = 10 and 2.4·10⁻¹⁷ at κ = 10¹³, worst 1.45·10⁻¹⁶ anywhere on the sweep — while the forward error climbs from 4.25·10⁻¹⁶ to 5.72·10⁻⁵. At the right-hand end the two are a factor of 2.38·10¹² apart, and nothing about the computation that produced them differs.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 1 Two errors, for matrices built to a prescribed condition number. The flat line along the bottom is how well the algorithm did — unchanged across twelve orders of magnitude of difficulty. The rising line is how wrong the answer is. Both are measured, twenty matrices at each condition number, and the distance between them is the condition number itself.

Start anywhere

one essay from each of the 22 fields

[½, 1)[1, 2)[2, 4)0.5124gap 0.125gap 0.25 — twice as wide8 values per octavespacing doubles at each power of two The arithmetic underneath · 19 essays

What a float can hold

The representable numbers are not a fine fuzz spread evenly over the line. They are evenly spaced inside each power-of-two interval and twice as far apart in the next one up, and almost everything else in this subject is a consequence of that one fact.

8 figures 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 Two errors, and whose fault they are · 18 essays

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.

6 figures 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, and the swap · 16 essays

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.

8 figures A = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle Orthogonality, measured · 16 essays

Orthogonal is a number

"Q is orthogonal" is a claim about a measurable quantity, ‖QᵀQ − I‖, and on the eight-by-eight Hilbert matrix two standard algorithms return 10⁻¹⁵ and 1 for it. The one that returns 1 still reconstructs the matrix perfectly, which is why nothing warns you.

8 figures everything Ax can reachb = (1.1, 0.4, 1.5)Ax, the closest reachable pointr = b − Ax‖Aᵀr‖ / (‖A‖‖r‖)1.7·10⁻¹⁶‖b‖² − ‖Ax‖² − ‖r‖²1.3·10⁻¹⁵‖r‖1.3200 random nearby points of the plane were tried; none is closer.a 3×2 system, Householder QRperpendicularity is checked Least squares, and the road not to take · 16 essays

The projection and the right angle

The least-squares solution is the one whose residual is perpendicular to everything the columns can reach. That is not a mnemonic — it is an equation, Aᵀr = 0, and the computed answer satisfies it to 10⁻¹⁶.

7 figures 081624324048566400.250.50.751index kfilter factor fₖno regularisation: fₖ = 1truncationTikhonovthe same sum, three weightsTikhonov, relative error0.11truncation, relative error0.11no filter at all5.5·10⁸both filters are one expression with a different weightfₖ = 1 is the catastrophe Regularisation, and the answer that is chosen · 15 essays

When the answer is a choice

A backward-stable least-squares solve of this problem returns an answer whose relative error is 5.5·10⁸. Nothing went wrong. The singular values decay exponentially with no gap anywhere in them, the data does not determine the answer, and something outside the data has to choose — which is the computation rather than a preliminary to it.

8 figures 10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹size of the perturbation ‖δA‖how far the eigenvalues moveJordan block, ε^(1/8)symmetric, ≤ ‖δA‖rounding error alone moves it to 10⁻²six seeds per symmetric point; Jordan is closed formsymmetry beats precision Eigenvalues, singular values, rank · 32 essays

Symmetry is worth more than precision

A symmetric matrix gives up its eigenvalues to full accuracy however ill-conditioned it is. An unsymmetric one can move them by the eighth root of a perturbation, so the rounding involved in merely storing the matrix shifts the spectrum by a hundredth.

6 figures 00.0939569-2-1.24498-0.4899610.2650581.020081.7751real partimaginary parttwo routes, one spectrumeigenvalues16rows8complex16against the closed form2.1·10⁻¹⁵n rows and 2n eigenvaluesso the eigenvectors are not a basis The eigenvalue problem that is not linear · 20 essays

A matrix that depends on its own eigenvalue

A damped structure does not produce Ax = λx. It produces (λ²M + λC + K)x = 0, where the matrix whose null vector is wanted is a function of the number being solved for — so there is nothing to factorise, an n × n problem has 2n answers, and the eigenvectors cannot be a basis.

7 figures 10⁻¹110¹10²10³10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹frequency ω|H(iω)|one function, two routesstates24relative, ω ≤ 10010⁻¹⁴against the terms3.9·10⁻¹⁵‖H‖∞0.0924 statesone function of one variable Reduction, and what a model is for · 16 essays

A model that is a rational function

A state matrix has a hundred thousand rows and the thing anyone wants from it is a function of one complex variable. The number that says how much of that size was ever the complexity is a rank — and the rank a derivation writes down cannot be computed, while one built from samples alone can.

7 figures 0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖ ⁄ ‖e₀‖ in the A-normmeasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417 Iterating, instead of factorising · 39 essays

The rate the condition number predicts

Conjugate gradients converge at a rate governed by the square root of the condition number. That is a bound rather than an estimate, it is provable, and it is loose enough that provisioning iterations from it wastes nine out of ten.

6 figures the matrix, lower triangle408 entriesits Cholesky factor1739 entries · 1331 created→‖A − LLᵀ‖/‖A‖1.4·10⁻¹⁶fill, symbolic1331fill, numeric1331n = 144 · density 3.2% · bandwidth 12same matrix, renumberedthe answer is identical to rounding Sparsity, and what elimination costs · 17 essays

The factor is not sparse

A sparse matrix has a factor that is not sparse, and the gap between them is the entire reason iterative methods exist. The entries elimination creates can be counted before any arithmetic runs, from the graph alone.

5 figures 024681012141610⁻¹1index kλthe transformthe eigensolverC = F* Λ F is a factorisationworst relative disagreement2.6·10⁻¹⁵‖Cx − b‖/‖b‖ from the transform solve4.7·10⁻¹⁶imaginary part of a real spectrum1.2·10⁻¹⁶n = 16, and the whole matrix is 16 numberseigenvectors known in advance Structure, and the solver that cannot see it · 15 essays

The matrix that is one row

A circulant of size 16 is sixteen numbers, has no zero entry anywhere, and hands over its entire spectrum in closed form — the discrete Fourier transform of its first column, exactly. An eigensolver spends a sweep of Jacobi rotations over 256 entries arriving at the same answer, and agrees to 1.2·10⁻¹⁵.

5 figures 051015202530354010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, in orderσ ⁄ σ₁eight digitsan independent draw per entry1 ⁄ ra cliff, and a control1/r rank at 10⁻⁸5log r rank at 10⁻⁸5noise rank at 10⁻⁸96σ₂ ⁄ σ₁0.024σ₆ ⁄ σ₁2.8·10⁻⁹the block has full rankand five useful columns Neither sparse nor dense · 18 essays

A block nobody can call sparse

A 96 × 96 block of a kernel matrix has ninety-six nonzero singular values and five that matter. It has no zero entries, it is not described by fewer numbers than it contains, and neither of the two ways this collection already knows to make a large matrix affordable applies to it.

7 figures 03672108144180216024681012eigenvalues in orderλthe closed formmarks: the assembled matrix, decomposeda spectrum nobody computedrows of the matrix216numbers that describe it108λ smallest0.59λ largest11worst |computed − exact|7.1·10⁻¹³the matrix is never neededand neither is its decomposition When the index is a tuple · 15 essays

An index that is a pair

A discretisation on a two-dimensional grid of n points a side has n² unknowns and a matrix with n⁴ entries — 10⁸ at n = 100. What that matrix is instead is two Kronecker products of an n × n matrix, which is 2n² numbers, and nothing has been approximated: assembling it was the mistake.

7 figures the matrix, measuredvertices40edges223‖L·1‖∞0zero eigenvalues1components, by search1λ₂1.5laid out at its own eigenvectorsand the row sums are exactly zero The matrix that is a graph · 19 essays

A matrix with no numbers in it

A graph arrives as vertices and edges. Two different matrices can be built from it, they answer different questions, and one of them has a null vector that is exact — the only object on this site whose kernel is known before anything runs.

7 figures 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 Exact arithmetic, and what it costs instead · 15 essays

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.

7 figures 04812162010⁻¹10⁻⁰.⁵1target rank k‖A − Aₖ‖₂published boundrandomisedσₖ₊₁, optimalhow far apart the three areworst seed spread1.6bound / median at k = 125.9median / optimum at k = 121.960×60, 6 seeds, oversampling p = 5band is best to worst Randomised, and the guarantee that changes kind · 16 essays

A bound that holds with probability

Every other guarantee in this collection is deterministic. The randomised low-rank approximation offers one that holds with a probability, the seed changes the answer, and the honest figure is a band rather than a line.

7 figures 015304560759010512010⁻²10⁻¹110¹steprelative sizeleast error: 20discrepancy stop: 7errorresidualthe knob is an integerleast error, at step20error there0.14error at step 1206the residual falls at every stepthe error turns and keeps rising Methods that were designed apart · 14 essays

A parameter that counts steps

The regularisation field's knob is a positive real number chosen by one of three rules. The iterative field's is an integer nobody called a knob — where to stop. On the same problem the best step is 20 and the best λ is 0.025, and they reach 0.1426 and 0.1406.

10 figures 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 When the problem arrives again · 16 essays

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.

5 figures 10²10³10⁴10⁵matrix size ncountoperations, bothwords, unblockedwords, blocked (b = 6)the answer does not move‖PA − LU‖/‖A‖, unblocked2.8·10⁻¹⁶‖PA − LU‖/‖A‖, blocked2.8·10⁻¹⁶difference between them0the dashed curve is both orderings' operation countthe solid pair is what they cost Where the flop count stopped predicting the time · 16 essays

The same arithmetic at a different price

A blocked and an unblocked elimination perform 72,568 operations each — the same operations, associated differently — choose the same pivots, and return a factorisation identical to the last bit: ‖PA − LU‖/‖A‖ = 4.487946226420872·10⁻¹⁶ in both. One of them moves 41,332 words between fast and slow memory and the other moves 19,476.

6 figures 110¹10⁻¹¹10⁻⁸10⁻⁵pieces the vector was divided intodistance from the exact sum, relativethe published boundκ · uone vector, one algorithmdistinct answers21runs26spread, in ulps2.3·10⁷κ of the sum10⁸bound ÷ worst error2.6·10⁴nobody chose pand no answer is the answer The answer that depends on the machine · 20 essays

The same program, twice

One vector of 4,096 numbers, one summation algorithm, one precision, twenty-six runs — and twenty-one different answers. Nothing in the program chose between them, every one of them satisfies the textbook bound, and the exactly rounded answer is not among them.

7 figures -1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue4 negative10 positivecounted before it was formedpositive10negative4at zero0innermost ratio39the zero block is a theoremand so is the count either side of it The matrix a constraint makes · 16 essays

The zero that is not a missing entry

A constrained minimisation produces a matrix with a zero block, and the zero is a theorem rather than a sparsity pattern. No pivot order makes it positive definite, no precision changes that, and Cholesky does not fail somewhere on it — it fails at the first constraint row, on a number the problem already contained.

7 figures

All 404 essays, by field — the 22 above are where each field starts; the index carries the rest in reading order.

Other ways in

the same essays, indexed differently

All 22 fields Depth ladders What is refused here Concepts Every generator What's new Search

Threads running through

themes, not chapters

Where the algebra stops being the arithmetic

Every identity in linear algebra is exact, and every implementation of it is not. The gap is not a defect to be apologised for — it is the subject. QᵀQ is the identity in the algebra and is measurably not the identity in the arithmetic, and the difference between those two sentences is what this site is about.

241 essays

Whose fault is it

A good algorithm returns the exact answer to a nearby problem. So when the answer is wrong there are two possible authors, and they are separately measurable: the backward error is what the algorithm did, the condition number is what the problem did to it. Almost every essay here reports both.

179 essays

Identical algebra, different arithmetic

The subject supplies pairs of algorithms that a textbook derivation cannot tell apart and a computer can: classical against modified Gram–Schmidt, the normal equations against QR, elimination with and without a row swap. Each pair is one derivation and two behaviours, and the second is only visible if you run it.

81 essays

Measured, not assumed

No decomposition is drawn here without its residual printed beside it. ‖A − QR‖/‖A‖, ‖QᵀQ − I‖, ‖PA − LU‖/‖A‖ — the numbers are on the figures, and a figure that shows a factorisation without one is refused rather than drawn.

382 essays

Two routes to a number

A determinant from the pivots and from cofactor expansion. A singular value from one-sided Jacobi and from the eigenvalues of AᵀA. A Hilbert solution in floating point and in exact rational arithmetic. Neither route can confirm itself, and where they disagree the disagreement is the finding.

185 essays

Precision as the knob

The decompositions here are written against an arithmetic context rather than against plus and times, so the same code runs at any significand width from sixteen bits to fifty-three. Where a figure has a slider, the slider is usually the mantissa — and every frame it can reach was generated, and asserted, at build time.

33 essays

Exact ground truth

Most numerical demonstrations estimate their own error by computing the same thing more carefully, which is circular. The Hilbert matrix has a closed-form rational inverse, so for that whole class of problems the answer is known in BigInt integers and the error is measured rather than approximated.

116 essays

Claims that can fail

An assertion is evidence only if it can fail. Every check on this site is fed a counterexample somewhere in the build — the orthogonality test is shown a shear, the reconstruction test is shown the wrong matrix, the growth measurement is shown Wilkinson's matrix — and the build fails if any of them accepts.

176 essays