Field

Elimination, and the swap

Gaussian elimination is the first algorithm anybody learns and the row interchange is the part nobody is given a reason for. Here it is run without one, on a matrix where that fails — quietly, returning an answer of the right shape — and the growth factor that governs the whole story is plotted against a bound it never attains.
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.

[ ε 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¹⁷

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.

0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8

The bound that is never attained

Partial pivoting's stability guarantee permits the entries to double at every step — a factor of 5.5·10¹¹ at n = 40. The measured growth on random matrices of that size is about three. The gap is eleven orders of magnitude, and the guarantee is still worth having.

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

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¹¹110¹10²10³10⁴condition number of the matrixgrowth factorbound 2^11partial pivotingCholeskyno pivot to gain fromCholesky growth, every κ1Cholesky interchanges0partial pivoting, worst9the bound, 2^112048both eliminations reach the same growthand only one of them had to swap to get there

A factorisation with nothing to pivot for

Cholesky's growth factor is not bounded by one. It is equal to one, at every size and every condition number, and the two-line reason is why the algorithm needs no pivoting at all — not "usually gets away without it". Its only failure is the square root of a non-positive number, which is exactly the test for definiteness, and in floating point that test moves with the precision.

D from PAPᵀ = LDLᵀ — the shaded pairs are 2×2 pivots10⁻⁶0.749······0.749·········1.2·10⁻⁶1.4······1.4·········2.1·10⁻⁶0.549······0.549·········3.6·10⁻⁶0.614······0.614·three rules, one matrix‖PAPᵀ − LDLᵀ‖, blocks5.8·10⁻¹⁷‖PAPᵀ − LDLᵀ‖, diagonal3.1·10⁻¹¹growth, blocks1.3growth, diagonal5·10⁵the zero block is what the problem saysand one rule does not need it to be nonzero

When symmetry is not enough

The matrix [[0, 1], [1, 0]] is symmetric, nonsingular and perfectly conditioned, and there is no diagonal entry to pivot on. Every factorisation restricted to symmetric interchanges and one-by-one pivots fails on it, at any depth of searching, because every entry it could search is zero. The repair is to take two variables at once.

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.

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.

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.

0510152025303540012345matrix size ngrowth factorpartial pivotingrook pivotingcomplete pivotingsolid: median of 30 · dashed: worst of themat n = 40partial pivoting, median3.3rook pivoting, median2complete pivoting, median1.6rook, worst of the draw2.8the same matrices under every rulerook 3.3× partial's search

A pivot that searches one row and one column

Rook pivoting looks down a column for its largest entry, along that entry's row for a larger one, and back down that entry's column, until it finds an entry largest in both. On Gaussian matrices of size 64 it keeps the median growth factor at 2.53 against partial pivoting's 4.06 and complete pivoting's 1.88, and it compares 6,987 entries against 2,080 and 89,440. On Wilkinson's matrix it holds the growth at exactly 2 where partial pivoting reaches 9.2·10¹⁸. And on a matrix built to make it walk it compares 113,376 entries — more than complete pivoting.

10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹110²10⁴standard deviation of the noise σmedian growth factormargin 0.00564shootingWilkinsonwhat the growth rests onWilkinson, no noise5.5·10¹¹Wilkinson, σ = 10⁻¹⁴2shooting, σ = 10⁻⁸1.1·10⁴pivot margin0.0056thirty draws at every pointa tie breaks at any noise; a margin needs its own size

A worst case is as fragile as its margin

Wilkinson's matrix grows by 5.5·10¹¹ under partial pivoting, and adding Gaussian noise of 10⁻¹⁴ to every entry takes its median growth to exactly 2. The shooting matrix from a boundary-value problem grows by 1.1·10⁴, and noise ten million times larger leaves it untouched. The difference is what each worst case rests on. Wilkinson's rests on exact ties between candidate pivots, which any noise breaks. The shooting matrix's rests on a choice made by a margin of 5.6·10⁻³, and between 10⁻⁶ and 10⁻³ its median growth is that margin divided by the noise, times a constant between one half and four thirds.

0246810121416110¹10²10³10⁴10⁵10⁶length of the interval Tgrowth factorpartialrookcompleteh = 0.3, 102 unknownspartial, T = 15.01.3·10⁵e^(5T/6)/21.3·10⁵complete2largest κ8.3a boundary-value problem, not a constructionκ single-digit throughout

The growth a boundary-value problem supplies

Large growth under partial pivoting is usually said to need a matrix built for it. A two-point boundary-value problem solved by multiple shooting supplies one without being asked: its growth factor is e^(5T/6)/2 to four figures — 1.1·10⁴ at an interval of 12, 1.3·10⁵ at 15 — on a matrix whose condition number never exceeds 8.3, while rook and complete pivoting keep it below 2. And it is the finer shooting grid that grows: below a step of 0.3397 the choice partial pivoting makes turns on one entry against one, and above it there is no growth at all.

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.

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.

10⁻⁹10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²110¹10²10³10⁴standard deviation of the noise σmedian growth factorthe marginevery entrynonzeros onlyrelative, per entryone error in Estep 0.3, margin 0.00564, thirty drawsno noise1.1·10⁴every entry, σ = 10⁻⁵536nonzeros only, σ = 10⁻³1.1·10⁴one error in E, σ = 10⁻², draws kept21every curve is the same matrixonly what the noise touches changes

Noise the growth amplifies

The shooting matrix's growth of 1.1·10⁴ fell under noise as its pivot margin divided by the noise, and the explanation offered was noise reversing partial pivoting's choices. But noise of 10⁻⁵ is five hundred times smaller than that margin. Add the same noise only to the entries that are not zero and every draw keeps the whole growth up to 10⁻³. The dense noise was not reversing the comparisons by itself: it sat in the zeros, the elimination multiplied it by the growth already made, and a comparison flips when that product reaches about four tenths of the margin — at every noise level from 10⁻⁶ to 10⁻³.

10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³r½ recorded during one factorisationnoise that halves the growthT 6, h 0.3: 75T 9, h 0.3: 910T 12, h 0.3: 11000T 15, h 0.3: 1.3·10⁵T 12, h 0.25: 11000T 12, h 0.2: 11000T 9, h 0.15: 910measured ÷ predicted, shooting matricessmallest ratio0.92largest ratio2.4Wilkinson: recorded r½0Wilkinson: noise that halves it10⁻¹⁶one factorisation, one extra comparison a stepdashed: a factor of three either way

A margin the factorisation records

Partial pivoting finds the second-largest candidate in every column it scans, and throws it away. Keep it, divide the gap by the growth reached at that step, and take the smallest over the steps before half the final growth has arrived. That one number, recorded by the factorisation that is already running, predicts the dense noise that halves the growth to within a factor of 0.92 to 2.4 on shooting matrices whose growth runs from 75 to 1.3·10⁵ and whose fragility spans three and a half decades — and on Wilkinson's matrix it is exactly zero.

noise that halves the growthshooting, partial pivoting0.0018shooting, τ = 0.10.56Wilkinson, τ = 0.10.1800.10.20.30.40.50.60.70.80.9110⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1one minus the thresholdstored noise that halves the growthshooting, T = 12, h = 0.3Wilkinson, n = 24the same growth at every thresholdheld by a choice noise must push further to reverse

A threshold that holds the growth still

Partial pivoting's large growth on the shooting matrix rests on a margin of 5.6·10⁻³, and Wilkinson's rests on exact ties, which a perturbation at rounding breaks. A sparse code pivots with a threshold instead, taking the sparsest row among candidates within a factor τ of the largest. At every threshold from 1 to 0.1 both matrices grow by exactly as much as under partial pivoting. But the growth is now held by row counts, which no perturbation of the values can reverse: at τ = 0.1 it survives noise on the stored entries up to 0.56 on the shooting matrix, three hundred times more, and 0.18 on Wilkinson's, where partial pivoting's goes at 10⁻¹⁶. The margin that predicts it is the chosen candidate's height above the threshold line, and it needs no division by the growth.

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.

step 0.3rate per interval1.3largest r that catches it1.3110¹10²10³10⁴trigger sensitivity rgrowth factor1.021.324dashed: e to the five-sixths of the stepthe trigger must be finer than the growth

A trigger finer than the growth

Threshold pivoting cannot remove the shooting matrix's growth of 1.1·10⁴, because at τ = 1 it is partial pivoting. The proposal was a factorisation that notices growth on a step and searches as rook pivoting does after it, paying for the stricter search only where growth is made. But no step makes much: the largest one-step rise is exactly the growing mode's rate over one shooting interval, 1.284 at a step of 0.3, and 1.133 at 0.15. A trigger finer than that rate holds the growth at 2.00 for about two thirds of rook's extra search; one coarser never fires. And a trigger fine enough searches as rook on a quarter of an ordinary random matrix's steps.

at ε = 10⁻¹⁰Bunch–Kaufman, |L||D||Lᵀ| ÷ A6.7rook (bounded), |L||D||Lᵀ| ÷ A9.310⁻¹⁰10⁻⁹10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹110²10⁴10⁶10⁸10¹⁰coupling εsizeBunch–Kaufman: multiplierBunch–Kaufman: productrook (bounded): multiplierrook (bounded): productdashed: ‖|L||D||Lᵀ|‖ ÷ ‖A‖the multipliers grow and the product does not

Where the multipliers go

Bunch–Kaufman bounds the growth in D and not the entries of L, and the warning attached to that is that everything which later uses the factors inherits the size of L. On a matrix built to make those entries 1.2 over ε, they reach 1.2·10¹⁰ while the solve's backward error stays at 1.6·10⁻¹⁶, |L||D||Lᵀ| stays at 6.7 times ‖A‖, and a step of refinement changes nothing. The large multipliers are where the rule has put the matrix's ill-conditioning. A direction of negative curvature read from those factors finds 7·10⁻¹⁶ of the curvature that is there; the bounded rule's finds 13%.

fractions of the curvaturebounded rule's direction, every ε0.14nearest zero at ε = 10⁻⁸, ÷ λ10⁻¹⁶10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹coupling ε|quotient| ÷ |smallest eigenvalue|Bunch–Kaufman's directionthe bounded rule'seigenvalue nearest zeroinverse iteration goes where the eigenvalue nearest zero isa solve amplifies the smallest, not the most negative

A curvature direction the factors cannot refine

A direction of negative curvature read from Bunch–Kaufman's factors held 7·10⁻¹⁶ of the curvature that was there, and the bounded rule's held 14 per cent. The prediction was that a few steps of inverse iteration with the same factors would recover it from either. One step leaves both under a thousandth, and two put both on positive curvature, at the eigenvalue nearest zero — because a solve amplifies the smallest eigenvalue in magnitude, not the most negative. What recovers the curvature is the matrix, not its factors: Lanczos from either direction reaches ninety-nine per cent in six to nine products at every coupling, and power iteration from Bunch–Kaufman's direction has not reached a tenth after sixty.

thirty matricesBunch–Kaufman: unpushed, usable on12bounded (rook): unpushed, usable on511.251.51.7522.252.5051015202530push: the pivot multiplied bymatrices on which it is nearest the smallest eigenvalueBunch–Kaufmanbounded (rook)dashed: all thirtyno small push is reliable

The push a pivot needs belongs to the matrix

An indefinite factorisation's most negative pivot looks like an estimate of the most negative eigenvalue, and on one planted family Bunch–Kaufman's, pushed a tenth further, was a usable shift for inverse iteration. On thirty random symmetric matrices neither rule's pivot is: unpushed it is nearest the smallest eigenvalue on 12 and 5 of them, pushed a tenth on 13 and 9, and only a push of two serves 29. The push each matrix needs runs from 0.59 to 2.27 for Bunch–Kaufman — whose pivot is more negative than every eigenvalue on eight matrices and under half the smallest on others — and from 0.85 to 2.01 for the bounded rule. What does work needs no push at all: multiply the shift by a quarter until the factorisation of A − σI has no negative pivot, which Sylvester's law says is exactly when σ is below the spectrum. It lands every time, in a median of two or three factorisations.

468101214160246810nlog₂ of the growth factorW ⊗ I₁W ⊗ I₂W ⊗ I₃look 1 step aheadlook 2 steps aheadlook 3 steps aheadlook 4 steps aheaddashed: 2^(n/d − 1); each family defeats every depth up to its owna fixed depth stretches the exponential

A doubling that arrives two steps late

A pivot rule that looks one elimination ahead chooses Wilkinson's growth row for row, because the doubling it sets up arrives one step later; looking two ahead escapes it. The argument said a matrix whose damage arrives two steps late would defeat the two-step rule the same way. It exists, and it is two Wilkinson matrices interleaved: W ⊗ I₂. Every look-ahead of depth one or two reproduces the greedy rule's growth on it exactly, 4 to 64 from n = 6 to 14, and depth three gives 2. Three copies defeat depth three and yield to four. The growth each depth is held to is 2^(n/d − 1): a fixed depth stretches the exponential and never removes it, at a cost that rises by a factor of about n for every step it looks.

10¹10²10³10⁴size n = m·dentries compared by the pivot search4102040complete pivotingrook pivotinghalf rookpartial pivotingone dot per matrix and rulethe rook's price does not see the deferral

The doubling is in the row already

Interleave d copies of Wilkinson's matrix and every pivot rule that looks up to d eliminations ahead walks into growth 2^(n/d − 1); escaping costs a look one step deeper, and a deeper look multiplies its trial eliminations. Rook pivoting looks no further than the present step and escapes every one of them — growth exactly 2 at every deferral and every size — for 3.2 times partial pivoting's comparisons whatever the deferral. The doubling a look-ahead hunts for in the future is already visible in the row in front of it. And the leak that was supposed to let the two-step rule escape was placed where no pivot is ever chosen: at every size up to one, it changes nothing.

1251020shifted symmetricsmallest pivot, plainsmallest pivot, diagonal firstharmonic mean of the spectrumGram, 48 rowssmallest pivot, plainsmallest pivot, diagonal firstharmonic mean of the spectrumGram, 30 rowssmallest pivot, plainsmallest pivot, diagonal firstharmonic mean of the spectrumover the smallest eigenvaluethirty matrices a family, logarithmic axisthe pivot follows the mean, not the bottom

A definite pivot sits at the mean

On a positive-definite matrix every Cholesky pivot is at or above the smallest eigenvalue, so a shift read off the smallest pivot always has to be pushed down, and it was predicted to need a median of four pushes and to be usable unpushed on none. The second half holds on all ninety matrices measured. The first does not: shrinking by a quarter at a time takes a median of six to ten steps, because the smallest pivot sits within a factor of two of the spectrum's harmonic mean, and the harmonic mean sits four to nine times above the smallest eigenvalue. The push was a device for the indefinite case, which has no floor to aim at. A definite matrix does, and halving towards zero gets below the spectrum in a median of two to four.

11.2251020one weak edgesmallest pivot, plainsmallest pivot, diagonal firstbest vertex left to the enda vertex on a weak tethersmallest pivot, plainsmallest pivot, diagonal firstbest vertex left to the enda well three widesmallest pivot, plainsmallest pivot, diagonal firstbest vertex left to the endover the smallest eigenvalueweak weight one thousandth, logarithmic axisa vertex, not an order

A vertex, not an order

A definite matrix's smallest Cholesky pivot sits at the harmonic mean of its spectrum, and the prediction was that a localised bottom would change that for the pivoted factorisation alone: on a path with one weak edge, taking the largest diagonal first would land within a factor of 1.2 of the smallest eigenvalue while the unpivoted pivot stayed at the mean. On the prediction's own family neither rule comes within a factor of three, because a weak edge does not localise the bottom on a vertex — it splits it over the shorter piece, and the best pivot any order can reach is the smallest eigenvalue times that piece's length. Where the bottom really does sit on one vertex, both rules land within two per cent, the unpivoted one wherever in the order that vertex happens to fall. The pivot counts the vertices the bottom lives on; the order chooses nothing.

All essays