Numerical Linear Algebra
where the algebra stops being the arithmetic
Concept

LU factorisation — where it appears

Writing a matrix as a permutation times a lower triangular times an upper triangular, which is what a general solve computes and reuses. It is what a general solve computes and reuses across right-hand sides, and its residual against the matrix is what says whether the pivoting worked.

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

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

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.

cost · Blocking 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 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 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 the estimator maximises this quantity over the columns it visitscolumn 1 ‹visited›12column 2 ‹the answer›114column 311.4column 411.4column 511.4column 611.4column 711.4column 811.4column 911.4column 1011.4column 1111.4column 1211.4estimate 12.0a walk that stopped earlythe estimate returned12the true 1-norm114columns visited1products with the matrix5the walk's own stopping test firedand every column it could see was smaller

An estimate that can be fooled

Nobody computes a condition number, because forming an inverse costs more than the solve did. Every library estimates it instead, from four or five products with a factorisation already in hand. The estimate is exactly right on four random matrices out of five — and there is a matrix, three distinct entries wide, on which it returns a twentieth of the truth.

error · Condition-estimation 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 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 does this matrix look nearly singular?green: the test agrees with the truth · red: it does not · the bar under each number is its magnitude, over sixty-two decades|det A||det A|^(1/n)σₘᵢₙ1/κ = σₘᵢₙ/σₘₐₓ0.1·I at n = 40perfectly conditioned10⁻⁴⁰0.10.11κ = 10¹⁰, |det| = 1nearly singular1110·10⁻⁶10·10⁻¹¹Hilbert at n = 8nearly singular2.7·10⁻³³8.5·10⁻⁵1.1·10⁻¹⁰6.6·10⁻¹¹the two counterexamplesκ of the scaled identity1its determinant10⁻⁴⁰κ of the normalised matrix10¹⁰its determinant1det(cA) = cⁿ det(A)so a determinant carries the units n times over

The number that decides nothing

The determinant is the first scalar anybody attaches to a matrix and the last one worth consulting. A tenth of the identity has a determinant of 10⁻⁶⁰ and a condition number of exactly one. The Hilbert matrix's determinant stops being right at n = 13 and stops being a number at n = 29, and nothing in between reports either.

error · Determinant 110¹10²10⁵block size bwords movedbest block: b = 10the recursion: no block size126,742 wordsbest block, scanned1.1·10⁵the recursion1.3·10⁵recursion ÷ best1.2M = 144 words; the recursion never reads it1.16× the best of 24 blocks

The recursion that was never told the memory

A blocked elimination has to be tuned to its fast memory, and tuned to one memory it costs up to 2.9 times the best at another. A recursive elimination splits the columns in half down to one and reads no memory size at all. On eight fast memories from 36 to 576 words it moves between 0.94 and 1.28 times the words of the best tuned block, with the same 585,200 operations and the same pivots — and on a machine with two caches it beats the block tuned to either cache on six machines of seven.

cost · Blocking 110¹123base case, columnswords ÷ the best block's√M − 2 = 10M = 144M = 144base of one1.2panels ≤ 10, worst1.2panels of 123the base case is a block size, rounded to a halvingand it has the block's cliff at √M − 2

The block size a recursion still has

A recursive elimination is sold as having no block size, and every real one switches to plain loops below some width. Swept over that width, the traffic is a staircase with its steps at the halvings of n, and its cliff sits where the blocked elimination's does — the first panel wider than √M − 2 moves 1.53 to 4.47 times the words, on six memories of six. On three caches the innermost decides, and a third cache costs every tuned block up to 14 per cent and the recursion nothing.

cost · Blocking 8×8 standard normalmatrices60median growth, fully pivoted1.3costliest single decision1.6step 1 — median cost1.624×step 1 — moves a row90%step 2 — median cost1.545×step 2 — moves a row87%step 3 — median cost1.419×step 3 — moves a row88%step 4 — median cost1.333×step 4 — moves a row80%step 5 — median cost1.058×step 5 — moves a row75%step 6 — median cost1.021×step 6 — moves a row80%step 7 — median cost1.000×step 7 — moves a row43%red: growth when that one search is skippedblue: how often that search moves a row

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 an ordering that is theregreedy at n = 764best ordering at n = 72its largest multiplier1456710¹10²matrix size ngrowth factorgreedy, ties forwardsgreedy, ties backwards, ε = 10⁻¹²greedy, ties backwards, ε = 0the best orderingthe best order is 2 3 4 5 6 1 at n = 6one cyclic shift of the rows

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 -12-10-8-6-4-2024681011.522.533.5leaf width minus (√M − 2)words ÷ pure recursionM = 64M = 144M = 256the edge is the same column at every memoryand the cheapest leaf sits on it

The leaf that sits on the edge

A recursive elimination's base case was found to be a block size in disguise, with a cliff where the blocked elimination's is, and the choice read as a trade: the processor wants wide leaves, the cache wants narrow ones. Measured at every width rather than at the halvings of 96, there is no trade inside the edge. Leaves exactly √M − 2 wide are the cheapest the recursion can have in words as well as calls — 0.81, 0.80 and 0.84 of the pure recursion's traffic at 64, 144 and 256 words — and one column wider moves 1.77 to 3.15 times it. And a matrix of 100 columns, which halves unevenly, meets the cliff in two steps rather than one.

cost · Blocking ties broken by 10⁻¹²one step ahead, n = 248.4·10⁶two steps ahead, n = 2424812162024110¹10²10³10⁴10⁵10⁶10⁷matrix size ngrowth factorgreedyone step aheadtwo steps aheadone step ahead lies exactly on greedythe damage is done a step before it shows

One step ahead is one step short

Partial pivoting takes the largest entry in the column and, on Wilkinson's matrix, walks into growth of 2^(n−1) that a cyclic shift of the rows avoids entirely. A rule that chose instead the pivot whose elimination leaves the smallest trailing submatrix was expected to see the good order at the first step. It sees nothing there: every first pivot leaves a largest entry of exactly 2, and with the ties broken by 10⁻¹² it prefers the greedy row by 10⁻¹². It attains 2^(n−1) at every size. Looking two eliminations ahead, the greedy row scores 4 and every other row 2, and the growth is 2 at every size up to 24. On random matrices one step of look-ahead helps below n = 16 and is worse than greedy on more than half of them by n = 32.

elimination · Elimination M = 144, edge 10halving, worst size1.1aligned, worst size0.86961041121201280.70.80.911.1columns nwords ÷ pure recursionhalving to the edgesplit at the edgedashed: the pure recursion's wordsone base case, two ways to reach it

Leaves cut to the edge on purpose

A recursive elimination's cheapest leaf is exactly the square root of M, less 2, columns wide, and the rule drawn from it was to set the base case there and let the halvings put the leaves at or below it. On thirty-two sizes from 96 to 127 columns, halving to that base case moves more words than the pure recursion on sixteen of them at 144 words of fast memory, because the halvings stop at six and seven columns, not ten. Cut every dimension at a multiple of the edge instead and every size keeps the saving: 0.79 to 0.86 of the pure recursion's words, against halving's 0.91 to 1.07, with the same arithmetic and the same pivots. What it cannot make full is the one leftover leaf, and a leftover of one column is where it loses.

cost · Blocking

Named alongside it

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

Flop countPartial pivotingBackward errorData movementBlocked algorithmCondition numberMemory hierarchyExact ground truthGaussian eliminationGrowth factorPermutationPivoting

All concepts