Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: 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.
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

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.

 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

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.

 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

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.

 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

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⁻¹⁵.

 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

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.

 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

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⁻¹⁶.

 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

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.

 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

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.

 -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

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.

 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

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.

 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

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.

 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

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.

 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

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.

 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

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.

 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

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.

 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

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.

 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

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.

 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

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.

 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

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.

 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

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.

 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

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.

 110¹10⁴10⁵block size bwords movedthe count: √(M/3) = 5measured best: b = 8words movedat the best block1.6·10⁴at b = 13.9·10⁴at b = 243.9·10⁴derived from M with no measurement, and scannedthe two agree Where the flop count stopped predicting the time

A block size is a property of the machine

Three lines of counting say the best block size is √(M/3). Scanned over every integer at five fast memories, the measured optimum is √M − 2 — exactly, at all five. The count has the right scaling and the wrong constant, low by a factor of 1.56, and the wrong form: the answer is affine in √M rather than proportional to it.

 10¹10²10²size ncondition numberlimit 81measureda limit, as a fraction of itselfreached at n = 1280.99still to go0.011κ at n = 8, as a fraction0.52every point is below the line and none of them is on itthe limit is not a value Structure, and the solver that cannot see it

A limit the matrix never reaches

Szegő's theorem gives a Toeplitz family's condition number in closed form — ((1+ρ)/(1−ρ))², which is 81 at ρ = 0.8. The 8×8 section reaches 52% of it, the 128×128 reaches 98.9%, and none of them ever arrives. A statement about a family is not a statement about the matrix in front of you.

 00.250.50.75110⁻¹110¹10²10³fraction of the method's own rangerelative errorfloor 0.141truncation KTikhonov λCGLS steprandomised rankfour methods, one floortruncation K0.14Tikhonov λ0.14CGLS step0.14randomised rank0.14four knobs from four fieldsand one obstruction underneath them Methods that were designed apart

Four knobs and one floor

A truncation, a Tikhonov parameter, a step count and a randomised rank, on one problem with an answer that is known. Their best errors are 0.1445, 0.1406, 0.1426 and 0.1449 — a spread of 3% across four methods that share no arithmetic.

All themes