Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

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

 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.

 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.

 natural1739reverse Cuthill–McKee1354minimum degree1026nested dissection1413matrix: 408 entries · dense factor: 10440bandwidth 12 · 4.26× the matrixbandwidth 12 · 3.32× the matrixbandwidth 123 · 2.51× the matrixbandwidth 108 · 3.46× the matrixn = 144, five-point stencilevery ordering fills in; none avoids it Sparsity, and what elimination costs

The order decides the memory

Four elimination orderings on one matrix give factors of 1,739, 1,354, 1,413 and 1,026 entries. All four factorisations are exact, all four return the same answer, and the one with the better asymptotics is not the one that wins.

 10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴ Least squares, and the road not to take

The road that squares the problem

The normal equations are the first method every course teaches and the method no library uses. Forming AᵀA squares the condition number, and below ε = √u it does not degrade — it produces a matrix that is exactly singular, from data that was perfectly usable.

 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 The matrix a constraint makes

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.

 10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at allbinary64 throughoutone function, two spellings · zero below 1.5·10⁻⁸ The arithmetic underneath

Cancellation takes the answer, not a digit

Subtracting two nearly equal numbers is exact. That is what makes it dangerous — the subtraction introduces no error at all, it exposes error the operands were already carrying, and the exposure can consume every significant figure at once.

 [ ε 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¹⁷ Elimination, and the swap

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.

 for each previous column i, subtract the projection of column j onto qᵢclassicalr[i][j] = qᵢ · a[j] ↑ the ORIGINAL columnv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 7 and 8: 1columns 6 and 8: 0.13columns 6 and 7: 0.13modifiedr[i][j] = qᵢ · v ↑ what is LEFT of itv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 1 and 8: 4.4·10⁻⁷columns 2 and 8: 2.7·10⁻⁷columns 3 and 8: 2.4·10⁻⁸The two R factors agree to 1.2·10⁻⁶ relative. The two Q factors do not.the 8×8 Hilbert matrixone word, eight orders Orthogonality, measured

Two Gram–Schmidts

One argument changes. Classical Gram–Schmidt projects the original column onto each previous direction; modified projects what is left of it. In exact arithmetic the coefficients are identical. In floating point they differ by eight orders of magnitude in the thing that matters.

 10⁻³10⁻²10⁻¹110¹12345678910conductance, and the two bounds on itpathcyclegridbarbelltwo blockshypercubepreferentialstarcompletethe bar is the inequalitythe dot is the graph The matrix that is a graph

Two Laplacians of one graph

The combinatorial Laplacian D − A and the normalised one, which conjugates it by the inverse square roots of the degrees, are built from the same object, are not similar to each other, and answer different questions. On a graph whose degrees are equal they coincide. On one whose degrees span an order of magnitude their second eigenvalues are sixteen times apart.

 05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here Iterating, instead of factorising

An orthogonalisation nobody calls one

Conjugate gradients are derived as a minimisation and behave as an orthogonalisation, which is why the finite-termination property in every textbook is not a property the method has in floating point.

 the mirrorx, length 4.000Hx = (-4.000, 0)v = x − αe₁safe sign: α = −‖x‖, so v is formed from a sum and nothing cancelsunsafe sign: α = +‖x‖ gives ‖v‖ only 35.1% of ‖x‖ + ‖x‖ — the digits go‖HᵀH − I‖5·10⁻¹⁶‖Hx‖ − ‖x‖8.9·10⁻¹⁶second component2.2·10⁻¹⁶built from a unit vectororthogonality is structural Orthogonality, measured

A reflection cannot stop being one

Householder QR holds orthogonality at 10⁻¹⁵ whatever the condition number of the matrix, and Gram–Schmidt does not. The reason is not that it is more careful. It is that its Q is built from unit vectors, and rounding a unit vector gives a different reflection rather than a broken one.

 0123456789101110⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110³Newton steptolerance asked for, and iterations paiditerations paidtolerance asked forouter residualthe adaptive policy, step by stepNewton steps10inner iterations, total1009first step's cost1last step's cost271final outer residual3.4·10⁻¹¹the rule reads the last two residualsand asks for nothing it cannot use When the problem arrives again

A tolerance that reads its own residual

The cheapest constant forcing term costs 980 inner iterations and arrives with a hundred times the forward error of the dearest, which costs 9,358. A rule that sets each step's tolerance from the ratio of the last two residuals costs 1,009 and arrives with neither problem — and it is not a constant, so it does not appear on the curve the constants are compared on.

 the matrix43 entriestip eliminated first253 entriestip eliminated last43 entries‖A − LLᵀ‖/‖A‖, tip first1.4·10⁻¹⁶‖A − LLᵀ‖/‖A‖, tip last0dense factor is n(n+1)/2 = 253 · sparse factor is 2n − 1 = 43one row swapped to the endnothing numerical chose between them Sparsity, and what elimination costs

Two ends of the same arrow

One matrix, one row moved from the front of the elimination order to the back, and the factor goes from completely dense to no fill at all. Both factorisations are exact to rounding, and nothing numerical chose between them.

 1357911131510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1index kσₖ ÷ σ₁σ₁√utwo factorstheir productthe normal equations, againroutes agree to k =6product floors at9.7·10⁻¹⁰σ₁√u2.3·10⁻⁹κ(P)κ(Q)1.8·10³⁶√ of it1.3·10¹⁸do not form the productthe σ below the line are the bound Reduction, and what a model is for

The product nobody had to form

The Hankel singular values are the square roots of the eigenvalues of PQ. Form that product and half of them stop existing, at a floor this site can predict from one number — and the fix is the one the least-squares field has had since its first essay, arriving in a place with no least-squares problem in it.

 051210241536204810⁻¹⁷10⁻¹⁴terms consumederror accumulated so farone accumulator8 piecesthe steps are the partial sumspeak partial sum32the answer1.6·10⁻⁵error, one piece7.6·10⁻¹⁴error, 8 pieces2.4·10⁻¹⁴mean at p = 648.1·10⁻¹⁵the walk sets the sizeand nothing sets the value The answer that depends on the machine

Where the disagreement comes from

The error of a reduction is a walk whose step length is the spacing of the running total, not of the answer. That one sentence predicts the size of the disagreement to a factor of two, explains why dividing the work makes it smaller, and explains why the value cannot be predicted at all.

 051015202530354010⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹iteration‖r‖ / ‖b‖plain CGIC(0) CGwhat the preconditioner didκ(A)48κ(L⁻¹AL⁻ᵀ)5.1‖A − LLᵀ‖/‖A‖0.0832D Laplacian, n = 100√κ ratio predicts 3.07× Iterating, instead of factorising

Changing the condition number on purpose

Preconditioning is usually introduced as a trick that makes an iteration converge faster. It is not a trick. It is solving a different system with the same solution and a condition number chosen rather than inherited, and the new condition number is computable.

 10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸ The arithmetic underneath

The order they are added in

Addition is associative in the algebra and is not associative in the arithmetic. The same million numbers, added in a different order, give answers that differ in the third significant figure — and the fix is not a wider float, it is a different order.

 classical Gram–Schmidt4.62·10⁻¹⁰modified Gram–Schmidt1.49·10⁻¹²Householder, one sweep2.03·10⁻¹⁴reduction tree, 16 leaves1.48·10⁻¹⁵departure from orthogonality, logarithmicthe tree, at four depths‖AᵀA − RᵀR‖/‖AᵀA‖, depth 13.4·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 24.3·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 31.7·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 41.5·10⁻¹⁵the same algebra, four timestwo of them are products of reflections Where the flop count stopped predicting the time

A reduction that changes the order

A tall-skinny QR computed as a tree of independent block factorisations touches a 512×12 matrix once instead of twelve times, computes a completely different sequence of roundings from the sweep it replaces, and returns ‖AᵀA − RᵀR‖/‖AᵀA‖ = 1.65·10⁻¹⁵ against the sweep's 9.95·10⁻¹⁵. On the same matrix classical Gram–Schmidt returns 4.6·10⁻¹⁰.

 the matrix105 entriescorner first — sparsest227 entries, growth 1.9·10¹¹largest first — safe242 entries, growth 1.19the middle factor is the smaller one, and its answer has no correct digitsboth factorisations reproduce the matrix‖PA − LU‖/‖A‖, sparsest3.8·10⁻¹⁷‖PA − LU‖/‖A‖, pivoted5.4·10⁻¹⁷forward error, sparsest3·10⁻⁵forward error, pivoted4.8·10⁻¹⁶red marks are entries elimination createdthe fill argument and the stability argument disagree Sparsity, and what elimination costs

Structure and stability stop being separable

The sparsest variable to eliminate on this matrix has a diagonal entry of 10⁻¹². Eliminating it produces the smaller factor, reproduces the matrix to 3.8·10⁻¹⁷ — better than pivoting does — and returns an answer wrong in the fifth digit.

 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 The eigenvalue problem that is not linear

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.

 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 The arithmetic underneath

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.

 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 Elimination, and the swap

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.

 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 Least squares, and the road not to take

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

All themes