The pivot that reads the units

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.

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.

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ᵢₖ| 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.

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.

Two questions this raises about the growth factor #

The growth factor is the quantity the stability of elimination is usually stated in — max|uᵢⱼ| / max|aᵢⱼ|, 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.

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 setting the scaling makes the answer better, and the figure’s assertion is written in two halves saying so.

Which stops, and why the stated reason is not the reason #

That is one stop, and the sweep says how many there are and whether the explanation predicts them:

    ε         scaled, no pivoting   scaled, partial   1 − 1/ε representable?
  10⁻¹¹            0                    0                   yes
  10⁻¹²            1.22·10⁻⁴            1.22·10⁻⁴           yes
  10⁻¹³            2.22·10⁻¹⁶           2.22·10⁻¹⁶          yes
  10⁻¹⁴            0                    0                   yes
  10⁻¹⁵            1.25·10⁻¹            1.25·10⁻¹           yes
  10⁻¹⁶            1.00                 1.00                no
  10⁻¹⁷            1.00                 1.00                no
  10⁻¹⁸            1.00                 1.00                no

Three of the eight lose nothing, not one — 10⁻¹¹, 10⁻¹³ and 10⁻¹⁴ — so it is nearer a third of the slider than a fifth.

And the criterion above does not sort them. 1 − 1/ε is exactly representable at five of the eight stops, and at two of those five the answer is lost anyway: 1.22·10⁻⁴ at 10⁻¹² and 1.25·10⁻¹ at 10⁻¹⁵. Exact representability of the pivot’s neighbour is necessary and it is not sufficient.

What decides the other two is the section immediately below this one, applied to the right-hand side rather than to the matrix. The scaled b is (1 + ε)/ε, and at ε = 10⁻¹² that evaluates to 1000000000001.0001221 against an exact 1000000000001. The absolute error is 1.22·10⁻⁴ — which is the measured forward error, to three digits. The elimination was exact and the data it was handed was not.

That is worth stating as a rule rather than as a curiosity: when a construction is exact, the next place to look is the right-hand side, and this essay’s own closing paragraph about dividing rather than multiplying is the same observation one column to the left.

The sweep also strengthens the essay’s central claim rather than qualifying it. The comparison above says the scaled system’s partial pivoting and the unscaled system’s no pivoting agree to within a factor of two. Inside the scaled system the two rules agree exactly, at all eight stops, because they perform the identical elimination — which is what “partial pivoting has been disarmed” means when it is said precisely.

What the exact stops are evidence for, which is less than it looks #

It would be easy to read three exact stops out of eight as a mitigation — the failure is intermittent, so perhaps it is rare. It is not, and the table says why.

The stops where nothing is lost are the ones where (1 + ε)/ε and 1 − 1/ε both happen to round to integers a double holds exactly. That is a property of the decimal ε somebody chose for a slider, not of the problem: at ε = 10⁻¹² the loss is 1.22·10⁻⁴ and at 10⁻¹³ it is at the rounding, and the difference between those two problems is nothing a modeller would notice or control.

So the honest reading is the opposite of a mitigation. The rule’s failure is not intermittent in any useful sense; the demonstration’s arithmetic is. Partial pivoting makes zero interchanges at every one of the eight stops, which is the defect, and whether that defect costs anything at a particular ε is decided by whether some intermediate landed on a representable integer. A code that got the right answer at 10⁻¹³ has not avoided the problem; it has been lucky about a rounding it did not choose.

That is the same distinction the growth factor section above draws between a rule behaving as specified and an elimination being stable, and it is why the pivot count is the signal to watch rather than the error. Zero interchanges is the defect, and it is there at every stop. Whether the answer survives is a separate question with a separate answer, and only one of the two is a property of the code.

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ᵢₖ|. A row scaling by D replaces that with |dᵢ aᵢₖ|. For any target sequence of pivots there exists a positive diagonal D producing it — take dᵢ 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.

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.

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.

The scaling changes the pivot choice at every stop: one interchange as given and none after, at all five ε. That is the mechanism, and it is invariant — the row scaling is not occasionally misleading partial pivoting, it is misleading it every time.

Complete pivoting is unmoved at every stop, which is the other half. It reads the same error before and after scaling at all five, because it searches the whole submatrix and a row scaling cannot change which entry of the whole submatrix is largest in the way it changes which entry of a column is.

And there is one ε on this slider at which none of it happens. At 10⁻¹⁴ the scaled answer’s relative error is exactly 0 — every route on the page is exact, and a reader who measured there would conclude that row scaling costs nothing. It is not at the easy end of the slider either: it sits between a stop where scaling costs four digits and a stop where it costs all of them.

That is worth putting on the page rather than smoothing over, because it is the failure mode this collection keeps finding in its own earlier work — a claim measured at the one parameter value somebody happened to be looking at. The general statement survives it: the error under scaling reads 1.2·10⁻⁴, 0, 1, 1 and 1, so four of the five stops are catastrophic and the exception is a single value at which the arithmetic happens to come out even. A defect that disappears at one sampled point has not disappeared, and the way to know that is to have sampled more than one.

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.