Numerical Linear Algebra
where the algebra stops being the arithmetic
Concept

Growth factor — where it appears

The largest entry produced anywhere in an elimination over the largest entry of the matrix, which is exactly one for Cholesky and 2ⁿ⁻¹ at worst for partial pivoting. It is the quantity every stability bound for elimination is written in, its worst case is attained on one constructed matrix, and in practice it is never large.

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

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 [ ε 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.

elimination · Pivoting 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.

elimination · Growth 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

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.

sparsity · Sparse pivoting 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.

elimination · Pivoting 10⁻³10⁻²10⁻¹1110¹10²10³10⁴pivot threshold τgrowth factor · entries in L+U, ÷ entries in Agrowthfillthe library default‖PA − LU‖/‖A‖ at τ = 0.12.9·10⁻¹⁶growth at τ = 0.138entries in L+U at τ = 0.1372one knob, two measurements, opposite directionsand the default is most of both

A threshold between fill and growth

One number decides how small a pivot an elimination will accept. At 0.001 the factor holds 172 entries and the matrix grows by 1,330; at 1 it holds 260 and grows by 1.2. The libraries ship 0.1, and the measurement says why.

sparsity · Sparse pivoting 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.

elimination · Cholesky -18-15-12-9-6-300log₁₀ ‖PKPᵀ − LDLᵀ‖ / ‖K‖share of orderingsbestworstexistence and stabilityfactorise1unregularised0.69worst growth6.4·10⁵growth × γ0.64the ordering is free to chooseand not free of consequence

The regularisation that legalises every order

Perturb a saddle-point matrix's two blocks in opposite directions and it acquires a factorisation with a diagonal D under every symmetric permutation — not under a good one, under all of them. Five hundred random orderings, five hundred successes, and a growth factor that spans six orders across them.

constraint · Quasi-definite 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.

elimination · Cholesky 02468101210⁻¹⁸10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸member of the sequencebackward error of the solvefresh order · equilibrated kept orderringed: another pivot replaced hereone order, two sets of unitsfresh order, last member9.5·10⁻¹⁷kept order, last member4.8·10⁻⁹equilibrated, last member6.5·10⁻¹⁷pivots replaced, kept2pivots replaced, equilibrated0the perturbation floor √u1.1·10⁻⁸a reused pivot order is safe or notdepending on what the rows are measured in

The order that was right last time

A pivot order computed once and reused across a sequence saves the symbolic phase, and the price is that a pivot which was large may now be small. Replacing it with √u·‖A‖ costs eight orders of backward error and iterative refinement recovers a factor of 8.8 of them. Divide each row by its largest entry first and the same reuse costs nothing at all.

sparsity · Sparse pivoting natural113 predicted · 113 countedminimum-degree63 predicted · 63 countedreverse Cuthill–McKee63 predicted · 63 countedthe shaded entries are fill: zeros of K that the factorisation makes nonzeroallocated before the numbersnatural113minimum-degree63reverse-cuthill-mckee63predicted minus counted0the symbolic phase decides the memoryand nothing later is allowed to argue

An ordering that does not wait for the numbers

A sparse factorisation's memory is decided by an ordering computed from the graph, and its stability by pivots computed from the values, and the two decisions fight. On one family of matrices they do not — the ordering can be chosen for fill alone, and the fill the symbolic phase predicts is the fill the factorisation produces — exactly, not as a bound.

sparsity · Sparse pivoting 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.

elimination · 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.

elimination · Growth 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.

elimination · Growth -18-15-12-9-6-300log₁₀ ‖PKPᵀ − LDLᵀ‖ / ‖K‖share of orderingsbestworstexistence and stabilityfactorise1unregularised0.67worst growth5·10⁷growth × γ0.5the ordering is free to chooseand not free of consequence

The perturbation that does the work

A saddle-point matrix made quasi-definite is perturbed in both blocks, and the laws measured for it moved both together. Moved apart, the laws all belong to one block. The zero block's perturbation γ decides whether every ordering factorises, sets the worst ordering's growth at 0.51/γ, and costs the answer 1,451 per unit — the reciprocal of the smallest eigenvalue of AH⁻¹Aᵀ to three figures. The perturbation of H moves none of the first two and costs 19 per unit. Refinement removes each block's perturbation at the rate its own Schur complement sets, so γ's limit sits fifty times nearer than δ's.

constraint · Quasi-definite 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 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⁻³.

elimination · Growth 10⁻³10⁻²10⁻¹1110¹10²10³10⁴pivot threshold τgrowth factor · entries in L + U ÷ entries in Agrowth, rows onlygrowth, rows and columnsfill ratio, rows onlyfill ratio, rows and columns8×8 grid, 64 unknowns, 288 entriesτ = 0.001: rows only — 741 entries; growth2209τ = 0.001: rows and columns — 640 entries; growth2.5τ = 0.1: rows only — 875 entries; growth38τ = 0.1: rows and columns — 640 entries; growth2.5τ = 1: rows only — 986 entries; growth1.2τ = 1: rows and columns — 659 entries; growth1.2dashed: a fixed column · solid: the column chosen toothe same threshold rule in both

The column that was never fixed

Every threshold-pivoting measurement so far chose the pivot row in a fixed column, and the routine's own description said that choosing the column as well would change the constants and not the argument. Measured, it changes the argument. On the 8×8 conflict grid the factor shrinks from 875 entries to 640 at the library default, and the growth factor that climbed to 2,209 as the threshold loosened stays at 2.54 at every threshold from 0.3 down to 0.001. What does most of the work is not the column but which of several equally cheap entries is taken — and on random sparse matrices, choosing the column without that makes the growth worse.

sparsity · Sparse pivoting 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.

elimination · Growth 80 matrices, τ = 0.1, ties to the largestno column search, median fill244one column136every column109share of the gain at one column0.804896144192240columns the search may look atmedian fill123468121624allno column search at allmedian fillmedian growthworst growth — see the captionthe fill is the quantity the width buysthe worst case is the draw

How few columns the search needs

A full row-and-column pivot search is quadratic in the active submatrix at every step, and no library performs one. Looking at a single sparsest column takes the median fill from 244 to 136 where the full search reaches 109 — four fifths of the benefit for a linear scan — and that share is 79, 83, 80, 89 and 92 per cent across five thresholds. The worst growth appears to favour the narrow search by a factor of six, and on the next draw it favours the wide one by two.

sparsity · Sparse pivoting a 28 × 28 saddle-point matrixnatural order, entries113sparsest order, entries72two-by-two pivots, loosest test0and at the strictest3023466992115the constant in the pivot testentries in the factor0.10.30.640.80.95133natural ordersparsest ordergrowth, naturalgrowth, sparsestnumbers above the points: two-by-two pivots takenboth curves favour the sparser order

The freedom a symmetric factorisation does not have

Permuting rows and columns together leaves no column to choose, so the conflict between the sparsest pivot and the sound one should be worse rather than better. On a saddle-point matrix whose constraint rows have no diagonal entry at all, it is not there: taking the sparsest available pivot holds 70 entries against the natural order's 113 and a growth of 1.28 against 1.83 — better on both currencies at once, at every setting of the pivot test. The two-by-two blocks that make it legal cost 1.33 entries apiece.

sparsity · Sparse pivoting 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.

elimination · Growth

Named alongside it

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

Gaussian eliminationPartial pivotingBackward errorFill-inPermutationThreshold pivotingComplete pivotingLDLᵀ factorisationWorst-case analysisSaddle-point systemsSymbolic factorisationBackward stability

All concepts