Numerical Linear Algebra
where the algebra stops being the arithmetic
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.
20 min read 10 figures Whose fault is itIdentical algebra, different arithmeticMeasured, not assertedAssertions that reject

Worth reading first: The swap that is not optional · The units the matrix is measured in.

The swap that is not optional is the essay this site wrote about why elimination needs a row interchange. Its example is the one everybody’s is:

[ ε  1 ] [x]   [1 + ε]
[ 1  1 ] [y] = [  2  ]

with ε far below the unit roundoff and the exact answer x = y = 1. Without pivoting the multiplier is 1/ε, the (2,2) entry becomes 1 − 1/ε, and at ε = 10⁻¹⁷ that rounds to −1/ε: the 1 is gone, and with it the whole of the second equation’s information. The computed answer is (0, 1), wrong by 100% in the first component.

Partial pivoting compares |ε| against |1|, swaps the rows, and returns the answer exactly.

That is where the essay stopped, and it is where the story usually stops. Here is what happens next.

Multiply the first row by 1/ε

Not a perturbation. Not an approximation. Multiply the first equation through by a constant, which is something anybody may do to any equation, and which changes nothing about the problem:

[ 1  1/ε ] [x]   [1/ε + 1]
[ 1   1  ] [y] = [   2   ]

The same x = y = 1 solves it, exactly.

Partial pivoting now compares |1| against |1|. There is nothing to choose, so it keeps the first row. The multiplier is 1, the (2,2) entry becomes 1 − 1/ε, and the 1 is gone.

It performs the identical elimination it was introduced to prevent, with pivoting fully switched on.

Four pivot rules on [[ε, 1], [1, 1]] with ε = 10⁻¹⁷, before and after its first row is scaled by 1/εEight bars of forward error. On the matrix as given, no pivoting loses the answer entirely and all three pivot rules return it to rounding. Multiply the first row by 1/ε — which does not change the solution — and partial pivoting makes no interchange at all and returns the same wrong answer that no pivoting does. Scaled partial pivoting and complete pivoting are unmoved.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
Fig. 1 Four pivot rules on the same 2×2, twice: as given, and with the first row multiplied by 1/ε. Each bar is the forward error against a true solution of exactly (1, 1), with the number of interchanges beside it. Partial pivoting makes one interchange on the left and none on the right. Drag the exponent of ε.

What the figure says

Eight bars, and the four on the left are the textbook result: no pivoting loses the answer; all three pivot rules recover it exactly.

The four on the right are the same four rules on the same problem in different units.

  • No pivoting: loses it, as before.
  • Partial pivoting: zero interchanges, and the same loss, to four significant figures. The assertion checks that: the scaled system’s partial pivoting and the unscaled system’s no pivoting agree on their forward error to within a factor of two, because they are the same elimination.
  • Scaled partial pivoting: one interchange, answer exact.
  • Complete pivoting: one interchange — a column interchange, not a row one — answer exact.

The last of those is a detail worth keeping rather than folding away. Complete pivoting searches the whole remaining submatrix, finds 1/ε sitting in the (1,2) position, and moves it to the pivot by swapping columns. It is a different repair reaching the same answer, and counting it as “a swap” alongside the scaled rule’s row interchange would hide that.

Four pivot rules on [[ε, 1], [1, 1]] with ε = 10⁻¹², before and after its first row is scaled by 1/εEight bars of forward error. On the matrix as given, no pivoting loses the answer entirely and all three pivot rules return it to rounding. Multiply the first row by 1/ε — which does not change the solution — and partial pivoting makes no interchange at all and returns the same wrong answer that no pivoting does. Scaled partial pivoting and complete pivoting are unmoved.forward error, relative to a solution of exactly (1, 1)no pivoting · as given1.3·10⁻⁴ 0 interchangesno pivoting · rows scaled1.2·10⁻⁴ 0 interchangespartial · as given2.2·10⁻¹⁶ 1 interchangepartial · rows scaled1.2·10⁻⁴ 0 interchangesscaled partial · as given2.2·10⁻¹⁶ 1 interchangescaled partial · rows scaled2.2·10⁻¹⁶ 1 interchangecomplete · as given1.1·10⁻¹⁶ 1 interchangecomplete · rows scaled2.2·10⁻¹⁶ 1 interchangethe same problem twicepartial, as given2.2·10⁻¹⁶partial, rows scaled1.2·10⁻⁴its relative residual1.2·10⁻¹⁶complete, rows scaled2.2·10⁻¹⁶the two systems have the same solutionand one pivot rule cannot see it
Fig. 2 The same eight bars at ε = 10⁻¹², where four digits survive rather than none. The pivot counts are identical to the figure above: one interchange on the left, none on the right, at every ε.

The rule reads numbers, and numbers carry units

The reason is one sentence and it survives being said plainly.

Partial pivoting’s rule is compare |a_ik| for i ≥ k and take the largest. Those entries are coefficients of different equations. If one equation is in newtons and another in kilonewtons, the comparison is between a number and a thousand times a number that means the same thing.

A comparison of magnitudes across rows is a comparison that depends on the units of the rows. There is no way to write a rule of that shape which is scale-invariant, and partial pivoting is exactly that shape.

Scaled partial pivoting fixes it by dividing each candidate by the largest entry in its own row, so what is compared is a dimensionless ratio. Complete pivoting fixes it by searching over columns too, which happens to find the large entry wherever the scaling put it. Both are invariant under row scaling, and both are asserted to be.

And the residual cannot see it

Here is the part that makes this a trap rather than a curiosity.

The wrong answer partial pivoting returns on the scaled system has a relative residual of 10⁻¹⁷. It is at rounding. The right answer’s residual is 0.

Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u40×40, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 3 The site’s standing measurement of the gap between what a solver can see and what a reader wants. The scaled 2×2 in this essay is that gap at its most extreme: a residual at rounding, and a forward error of exactly 1.

So the failure is invisible from the only quantity a solver without the true answer can compute. The computed x really is the exact solution of a nearby problem — nearby in the norm the residual is measured in — and that is precisely what backward stability promises. The promise is kept and the answer is wrong, which is a small residual is not a small error arriving through a door that essay did not mention.

The site’s identity accounts for it exactly: forward error ⪅ κ × backward error, and κ of the scaled matrix is 1/ε. The backward error is 10⁻¹⁷, the condition number is 10¹⁷, and the product is 1. Every term is doing what it says. What is not obvious from the identity is that the condition number in it was created by the scaling, which is the subject of the units the matrix is measured in.

Gaussian elimination on a 4×4, one step at a timeFour copies of the same matrix: as given, and after each of the three elimination steps. The pivot in use is outlined and the entries reduced to zero are greyed.2114-60-272as givenrows in the order 1 2 34-60041042after step 1pivot 44-60041001after step 2pivot 4‖PA − LU‖/‖A‖0largest multiplier1row order 2 1 3the pivot is chosen
Fig. 4 An elimination drawn as an elimination, which is how this site’s first essay in the field introduced the pivot. Every choice that picture makes is a comparison of magnitudes down a column.

Two questions this raises about the growth factor

The growth factor is the quantity the stability of elimination is usually stated in — max|u_ij| / max|a_ij|, bounded by 2ⁿ⁻¹ under partial pivoting, and this site has measured Wilkinson’s matrix attaining that bound exactly.

On the scaled 2×2 the growth factor under partial pivoting is about 1/ε, so the theory is not being violated: the bound of 2ⁿ⁻¹ = 2 does not hold, because a growth factor is a scale-invariant quantity only when the pivoting is. Under scaled partial pivoting the growth is 1.

Wilkinson's 12×12 matrix and its upper triangular factorThe matrix on the left has ones on the diagonal, minus ones below it and a column of ones at the right. On the right, its U factor, whose last column doubles down the rows to sixty-four.1··········1-11·········1-1-11········1-1-1-11·······1-1-1-1-11······1-1-1-1-1-11·····1-1-1-1-1-1-11····1-1-1-1-1-1-1-11···1-1-1-1-1-1-1-1-11··1-1-1-1-1-1-1-1-1-11·1-1-1-1-1-1-1-1-1-1-111-1-1-1-1-1-1-1-1-1-1-11A1··········1·1·········2··1········4···1·······8····1······16·····1·····32······1····64·······1···128········1··256·········1·512··········11024···········2048U‖PA − LU‖/‖A‖0growth factor2048the 2ⁿ⁻¹ bound2048Row interchanges performed: 0.Partial pivoting had nothing to choose.every entry is 0, 1 or −1the bound is attained here
Fig. 5 The matrix that attains the growth bound exactly, under the pivot rule that is supposed to prevent growth. It gets there by making zero interchanges — which is the same mechanism as this essay’s example, arrived at by construction rather than by a scaling.

The two examples are worth holding together. Wilkinson’s matrix makes partial pivoting choose the diagonal at every step because the matrix was built that way. The scaled 2×2 makes it choose the diagonal because a scaling put a large number there. In both cases the rule behaves exactly as specified and the elimination is unstable, and in both cases the number of interchanges is zero.

Zero interchanges is the signal, and no library reports it.

How much survives, and the one stop where nothing is lost

The drag on the hero figure runs ε from 10⁻¹² down to 10⁻¹⁸, and it is worth being clear about what it changes and what it does not.

What it does not change is the pivot counts. Partial pivoting makes one interchange on the unscaled system at every ε and none on the scaled one at every ε, because the comparison it performs — |1| against |1| — has no ε in it at all. The scaled and complete rules make one interchange each at every ε.

What it changes is how much of the answer is left. At ε = 10⁻¹² the elimination 1 − 10¹² keeps about four significant digits of the 1 and the forward error is 1.2·10⁻⁴; at 10⁻¹⁷ the 1 falls off the end and the error is exactly 1.

And there is one stop where nothing is lost at all, which the site’s own stop sweep found rather than its author.

At ε = 10⁻¹⁴ the scaled system’s partial pivoting is exact. 1 − 10¹⁴ is exactly representable — 99999999999999 fits in a double with room to spare — so the elimination introduces no error, and the back-substitution then divides by 1 rather than by ε. Meanwhile the unscaled system without pivoting still loses four digits, because its back-substitution computes ((1 + ε) − 1)/ε and the subtraction in the numerator cancels.

So at that one setting the scaling makes the answer better, by a factor of 10¹², and the figure’s assertion is written in two halves saying so. Asserting the loss unconditionally would have been asserting something false at a fifth of the slider’s positions — which is the shape of mistake this site has recorded three times now under a different name each time.

A detail about how the scaled matrix is formed

Worth a paragraph because it is the same class of error one level down.

The scaled row is built by dividing each entry by ε rather than by multiplying by 1/ε. Over the reals those are the same operation. In binary floating point they are not: 1/ε rounds, and ε · fl(1/ε) comes out as 1.0000000000000002 for ε = 10⁻¹⁸.

That is enough. Partial pivoting compares the (1,1) entry against the (2,1) entry, finds 1.0000000000000002 against 1, and makes the interchange the whole demonstration depends on it not making. One stop of the slider silently stopped demonstrating anything.

x/x is exactly 1 in IEEE-754 for every finite nonzero x, so dividing gives the intended matrix at every ε. The comment in lauchliPair says so, and the sweep that caught it is the same gate that found the other two defects in this phase.

The same failure, one dimension up

The 2×2 is a demonstration. The mechanism is not confined to it, and the general statement is short.

Partial pivoting selects, at step k, the row maximising |a_ik|. A row scaling by D replaces that with |d_i a_ik|. For any target sequence of pivots there exists a positive diagonal D producing it — take d_i large for the rows to be chosen early — so every pivot order is reachable by a row scaling, including the worst one.

The 2×2 is the smallest case in which the worst order is catastrophic. Larger cases are not harder to construct, only harder to read; the site draws the small one because the whole elimination fits in two lines of arithmetic and a reader can check it by hand, which is the point of the example everybody uses.

Elimination with and without pivoting, ε = 10⁻¹⁷The same two-by-two system solved twice. With a row swap the answer is exact; without one the upper triangular factor contains an entry of order one over epsilon and the second component of the answer is wrong.[ ε 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¹⁷
Fig. 6 The original: what the same 2×2 does with pivoting and without it, at ε = 10⁻¹⁷. Everything in this essay is the observation that a row scaling moves a problem from the left bar to the right one without changing the problem.
Growth factor under partial pivoting to n = 14: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.03691215110¹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 = 14: bound 8192, worst 4.8
Fig. 7 The growth factor against the size, with the bound partial pivoting is supposed to guarantee. A row scaling moves a matrix along this axis without changing the problem it represents.

What real codes do

Every serious dense solver equilibrates, and the reason is this essay.

LAPACK’s expert driver xGESVX computes row and column scalings, applies them if the row norms are extreme enough, factorises the scaled matrix, and reports the condition number of what it factorised. The simple driver xGESV does none of that. Both are one call, and the difference between them is a letter.

The equilibration is by powers of two, so the scaled matrix is the unscaled one with exponents adjusted and no rounding is introduced. equilibrateRowsByPowersOfTwo does the same thing here for the same reason: a repair that costs accuracy to apply has to be argued for, and this one does not.

Sparse solvers go further, because they have to. UMFPACK and MUMPS both scale as a matter of course, and both offer several scalings, because on a sparse matrix the pivot choice is also a fill choice and structure and stability stop being separable — the threshold in a sparse LU is a negotiation between two objectives, and a row scaling moves the boundary between them.

Two condition numbers of one 8×8 system, as its rows are put into different unitsFour curves against the spread of the row units, in decades. κ_∞ of the scaled matrix rises from 9.83 to 1.9·10⁸ while the componentwise condition number stays at 6.98 throughout — the same system, the same solution, and one of the two numbers is a fact about the units. Hilbert's two numbers are drawn flat beside them at 3.4·10¹⁰ and 1.2·10¹⁰: a matrix whose sensitivity no scaling repairs.02468110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 8 decades1.9·10⁸cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it
Fig. 8 The condition numbers of the badly scaled system this essay’s argument runs on. Equilibrating is the repair for both the κ in this figure and the pivot choice in the one above, which is why it is one routine rather than two.

Scaled partial pivoting, and why it is not the default

If dividing by the row’s largest entry fixes this, why does the dense default not do it?

Two answers, and both are honest.

It is equilibration by another name, done later. Computing the row scales once and dividing every candidate by them is exactly what equilibrating first and then pivoting normally would do. A library that equilibrates has already applied the fix, and applying it twice is not applying it more.

And the scales must be read once, from the original matrix. Recomputing them from the shrinking submatrix as the elimination proceeds is a different rule, and a worse one: the trailing submatrix’s row norms change as the elimination introduces fill, so the comparison drifts. luPivot reads them once, from A, before anything is eliminated, and the comment saying so is there because the other version is the natural way to write the loop.

Complete pivoting, and the reason nobody uses it

Complete pivoting is invariant under row scaling, invariant under column scaling, and has a growth bound that is polynomial rather than exponential. It is, on every stability criterion, the better rule. It is used essentially nowhere.

The reason is the search. Partial pivoting looks at n − k entries at step k, which is n²/2 comparisons over the whole elimination and is lost in the noise beside 2n³/3 arithmetic operations. Complete pivoting looks at (n − k)², which is n³/3 comparisons — the same order as the elimination itself, and comparisons do not vectorise, do not block, and do not run at the rate a multiply-add does.

So the cost is not “a bit more”; it is a solve whose runtime is dominated by a search. The same arithmetic at a different price argument applies with full force: the operation count barely moves and the time doubles or worse, because the search touches the whole trailing submatrix at every step and does nothing useful with it.

Rook pivoting is the compromise that gets used when partial pivoting is not trusted — search along the column, then along that entry’s row, alternating until an entry is largest in both. It is scale-invariant in the same way complete pivoting is, and its expected cost is a small multiple of partial pivoting’s rather than a factor of n.

What is worth noticing is that all three of these rules cost search, and the alternative — equilibrate first, then pivot partially — costs one pass over the matrix and no search at all. That is why the repair that shipped is the cheap one, and why it is applied before the elimination rather than inside it.

Four pivot rules on [[ε, 1], [1, 1]] with ε = 10⁻¹⁶, before and after its first row is scaled by 1/εEight bars of forward error. On the matrix as given, no pivoting loses the answer entirely and all three pivot rules return it to rounding. Multiply the first row by 1/ε — which does not change the solution — and partial pivoting makes no interchange at all and returns the same wrong answer that no pivoting does. Scaled partial pivoting and complete pivoting are unmoved.forward error, relative to a solution of exactly (1, 1)no pivoting · as given1.2 0 interchangesno pivoting · rows scaled1 0 interchangespartial · as given1.1·10⁻¹⁶ 1 interchangepartial · rows scaled1 0 interchangesscaled partial · as given1.1·10⁻¹⁶ 1 interchangescaled partial · rows scaled2.2·10⁻¹⁶ 1 interchangecomplete · as given2.2·10⁻¹⁶ 1 interchangecomplete · rows scaled2.2·10⁻¹⁶ 1 interchangethe same problem twicepartial, as given1.1·10⁻¹⁶partial, rows scaled1its relative residual5·10⁻¹⁷complete, rows scaled2.2·10⁻¹⁶the two systems have the same solutionand one pivot rule cannot see it
Fig. 9 At ε = 10⁻¹⁶ the scaled system’s partial pivoting is wrong by 100% and its residual is at rounding. Nothing a solver can compute distinguishes this bar from the one beside it.
Two perturbation bounds and the error that was measured, on a 8×8 matrix spread over 8 decades of unitsThree curves against the size of an entrywise relative perturbation. The normwise bound κ_∞·ε is a valid bound and sits 4·10⁷ times above the componentwise one cond(A, x)·ε, which is also a bound and is nearly attained by the worst of forty random perturbations at each size.10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ_∞ · εcond(A,x) · εmeasuredboth bounds holdκ_∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement
Fig. 10 The perturbation model the scaled system’s small residual is small in. A backward error measured normwise is satisfied by an answer that is wrong componentwise, which is what the bar chart above is a single instance of.

What is worth carrying

A pivot rule that compares across rows is reading the units. That is not a defect of partial pivoting in particular; it is a property of the shape of the rule, and no rule of that shape can be scale-invariant.

“Use partial pivoting and elimination is stable” is a statement about a matrix whose rows are comparable in size. That is nearly always true and it is not a theorem.

The failure is invisible from the residual, because the algorithm really is backward stable — what the scaling did was manufacture a condition number for the backward error to be multiplied by.

And zero interchanges is the tell. On a general matrix, partial pivoting that never swaps has either been handed a matrix that needs no swaps or has been handed a matrix whose scaling has already made the choice.

Where this goes next

[ ε 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. 0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it The units the matrix is measured in One linear system, written twice. The rows of the second are the rows of the first in different units, the solution is identical to the last bit, and the condition number has moved by eight orders of magnitude. One of those two numbers is a fact about the problem and the other is a fact about the notation. 110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method A small residual is not a small error Substituting the answer back and finding that it fits is the most natural check there is, and it verifies the wrong thing. A residual of 10⁻¹⁷ is entirely compatible with an answer whose second digit is wrong.

What links here

Computed from the collection, not written here: the essays that point at this one.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Backward errorComplete pivotingEquilibrationGaussian eliminationGrowth factorPartial pivotingRow scalingSilent failureWilkinson's matrix