🧮Linear Systems: Direct Methods

Pivoting, Norms and Condition Numbers

Partial pivoting, vector and matrix norms, κ(A), residual versus error, backward stability and refinement

A pessimistic prediction, 1943

In 1943 Harold Hotelling, one of the founders of mathematical statistics, analysed how rounding errors might accumulate in Gaussian elimination. His bound suggested that each extra equation could multiply the error by up to 4 — so a system of

n

equations would need about

0.6n

extra decimal digits just to break even. For 100 equations, 60 digits. He recommended iterative methods instead. If he had been right, large-scale scientific computing would never have happened. At Princeton, John von Neumann and Herman Goldstine produced a far more careful analysis in 1947. At the NPL, Turing pointed out in 1948 that in practice — including Fox's computations by hand — no such exponential build-up was ever observed. The full explanation came in 1961, when James Wilkinson showed that elimination with pivoting is backward stable: the computed answer is the exact answer to a problem changed by a few units of roundoff. This lesson is that explanation — and the one case where Hotelling's fear is real.

Why a small pivot is fatal. Eliminate in

\begin{bmatrix} \varepsilon & 1 \\ 1 & 1 \end{bmatrix}\mathbf x = \begin{bmatrix} 1 \\ 2 \end{bmatrix}

. The multiplier is

m = 1/\varepsilon

, and the new (2,2) entry is

1 - 1/\varepsilon

. For

\varepsilon = 10^{-17}

, the 1 is lost entirely in rounding:

\mathrm{fl}(1 - 10^{17}) = -10^{17}

. Back substitution gives

x_2 = 1

correctly, but then

x_1 = (1 - x_2)/\varepsilon = 0

— completely wrong.

The matrix is perfectly well behaved. The algorithm made a large number (

10^{17}

) out of a small one, and the rounding of that large number destroyed the information.

💡Partial pivoting

At stage

k

, before eliminating, swap rows so that the entry of largest magnitude in column

k

(on or below the diagonal) becomes the pivot. Every multiplier then satisfies

|m_{rk}| \le 1

, and no entry is amplified by more than a factor of 2 per stage. The row swaps are recorded in a permutation matrix:

PA = LU

. This is what np.linalg.solve, scipy.linalg.lu_factor and MATLAB's backslash do by default, and why a library solver returns

x_1 = x_2 = 1

for the example above.

Derivation

Why pivoting keeps the numbers in check

Stage

k

replaces each entry below and to the right of the pivot by

a_{ij} - m_{ik}a_{kj}

. With partial pivoting

|m_{ik}| \le 1

, so

|a_{ij} - m_{ik}a_{kj}| \le |a_{ij}| + |a_{kj}| \le 2\max_{p,q}|a_{pq}|.

No stage can more than double the largest entry. After

n - 1

stages the entries are at most

2^{n-1}

times those of

A

— the bound on the growth factor quoted below. Without pivoting

m_{ik} = a_{ik}/a_{kk}

is unbounded, and in the

\varepsilon

example it was

10^{17}

: the new entry

1 - 10^{17}

wiped out the 1, which was the information.

To say how accurate a solution is we need to measure the size of vectors and matrices. A norm does that:

$\|\mathbf x\|_2 = \sqrt{\sum x_i^2}$ , $\|\mathbf x\|_\infty = \max |x_i|$ , $\|\mathbf x\|_1 = \sum |x_i|$ .
The matrix norm induced by a vector norm is the largest stretch: $\|A\| = \max_{\mathbf x \ne 0} \|A\mathbf x\|/\|\mathbf x\|$ . Concretely, $\|A\|_1$ is the largest absolute column sum, $\|A\|_\infty$ the largest absolute row sum, and $\|A\|_2$ the largest singular value (Unit 5).

Derivation

Why

\|A\|_\infty

is the largest row sum

Take any

\mathbf x

with

\|\mathbf x\|_\infty = 1

, so every

|x_j| \le 1

. Each component of

A\mathbf x

is bounded by its row's absolute sum:

|(A\mathbf x)_i| = \Big|\sum_j a_{ij}x_j\Big| \le \sum_j |a_{ij}|\,|x_j| \le \sum_j |a_{ij}|.

So no vector is stretched by more than the largest row sum. And that much is achieved: in the row with the largest sum, choose

x_j = \pm1

with the sign of

a_{ij}

; then every term is

|a_{ij}|

and the bound is reached. The same argument with columns, for

\|\cdot\|_1

, gives the largest column sum.

💡The condition number of a matrix

Perturb

\mathbf b

\delta\mathbf b

. The solution changes by

\delta\mathbf x = A^{-1}\delta\mathbf b

, and a few lines of norm inequalities give

\dfrac{\|\delta\mathbf x\|}{\|\mathbf x\|} \;\le\; \kappa(A)\,\dfrac{\|\delta\mathbf b\|}{\|\mathbf b\|}, \qquad \kappa(A) = \|A\|\,\|A^{-1}\|.

The same

\kappa

governs perturbations of

A

itself. It is lesson 1-4's condition number for the problem 'solve

A\mathbf x = \mathbf b

': a property of the matrix, not of the algorithm.

\kappa \ge 1

always;

\kappa = \infty

for a singular matrix. Geometrically, in two dimensions, the equations are two lines and

\kappa

is large when they are nearly parallel.

Derivation

The condition-number bound, line by line

The exact solution satisfies

A\mathbf x = \mathbf b

; the perturbed one satisfies

A(\mathbf x + \delta\mathbf x) = \mathbf b + \delta\mathbf b

. Subtract:

A\,\delta\mathbf x = \delta\mathbf b \quad\Longrightarrow\quad \delta\mathbf x = A^{-1}\delta\mathbf b \quad\Longrightarrow\quad \|\delta\mathbf x\| \le \|A^{-1}\|\,\|\delta\mathbf b\|.

The last step is the definition of the matrix norm:

A^{-1}

stretches no vector by more than

\|A^{-1}\|

. The same definition applied to

\mathbf b = A\mathbf x

gives

\|\mathbf b\| \le \|A\|\,\|\mathbf x\|

, that is

1/\|\mathbf x\| \le \|A\|/\|\mathbf b\|

. Multiply the two inequalities:

\frac{\|\delta\mathbf x\|}{\|\mathbf x\|} \le \|A\|\,\|A^{-1}\|\,\frac{\|\delta\mathbf b\|}{\|\mathbf b\|} = \kappa(A)\,\frac{\|\delta\mathbf b\|}{\|\mathbf b\|}.

The residual bound is the same calculation. For a computed

\hat{\mathbf x}

\mathbf r = \mathbf b - A\hat{\mathbf x} = A(\mathbf x - \hat{\mathbf x})

, so

\hat{\mathbf x} - \mathbf x = -A^{-1}\mathbf r

: the residual plays the role of

\delta\mathbf b

, which gives the inequality in the equation below.

And $\kappa \ge 1$ :

1 = \|I\| = \|AA^{-1}\| \le \|A\|\,\|A^{-1}\|

\frac{\|\hat{\mathbf x} - \mathbf x\|}{\|\mathbf x\|} \;\le\; \kappa(A)\,\frac{\|\mathbf r\|}{\|\mathbf b\|}, \qquad \mathbf r = \mathbf b - A\hat{\mathbf x}

The residual r is cheap to compute and is always small for a stable solver. The error can still be κ(A) times larger: a small residual is not a small error.

📐Hilbert matrices

The Hilbert matrix

H_{ij} = 1/(i + j - 1)

arises when fitting polynomials by least squares in the monomial basis. Its condition number grows roughly like

e^{3.5n}

. Solve

H\mathbf x = H\mathbf 1

, whose exact answer is all ones (table below). For

n = 8

\kappa \approx 1.5 \times 10^{10}

and the error is about

10^{-7}

; for

n = 12

\kappa \approx 1.7 \times 10^{16}

and the answer has no correct digits at all — yet the residual is

10^{-16}

every time. The solver did its job perfectly. The problem itself has no more digits to give.

n	κ₂(H)	relative residual	relative error
4	$1.6 \times 10^{4}$	$\sim 10^{-16}$	$4 \times 10^{-14}$
6	$1.5 \times 10^{7}$	$\sim 10^{-16}$	$1 \times 10^{-10}$
8	$1.5 \times 10^{10}$	$\sim 10^{-16}$	$6 \times 10^{-8}$
10	$1.6 \times 10^{13}$	$\sim 10^{-16}$	$9 \times 10^{-5}$
12	$1.7 \times 10^{16}$	$\sim 10^{-16}$	$3 \times 10^{-1}$

Solving H x = H·1 in double precision with partial pivoting (np.linalg.solve)

One of the main virtues of an electronic computer from the point of view of the numerical analyst is its ability to "do arithmetic fast." Need the arithmetic be so bad!

– James H. Wilkinson, Some Comments from a Numerical Analyst, Turing Award lecture (1970)

💡Backward stability of elimination

Wilkinson's theorem (stated here without its proof, which tracks the standard model of lesson 1-2 through every stage of elimination — see Higham, chapter 9): Gaussian elimination with partial pivoting produces

\hat{\mathbf x}

that exactly solves

(A + \delta A)\hat{\mathbf x} = \mathbf b

with

\|\delta A\| \le c\,n\,\rho\,u\,\|A\|,

where

\rho

is the growth factor — the largest entry appearing during elimination divided by the largest entry of

A

. With pivoting,

\rho \le 2^{n-1}

. In theory that is Hotelling's exponential. In practice

\rho

is almost always below 10, so the backward error is a few units of roundoff and the forward error is about

\kappa(A)\,u

— the best any algorithm could do.

⚠The worst case is real, and rare

Wilkinson constructed matrices — ones on the diagonal,

-1

below it, ones in the last column — for which partial pivoting does no swaps and the last column doubles at every stage:

\rho = 2^{n-1}

. At

n = 60

that is

5.8 \times 10^{17}

, and the computed solution is wrong in the first digit even though

\kappa

is modest. Over decades of practice, such growth has almost never occurred in real problems; random

60 \times 60

matrices give

\rho

below 10. When certainty matters, complete pivoting or a QR factorisation (Unit 5) removes the risk at extra cost.

Iterative refinement. Given a computed

\hat{\mathbf x}

and the LU factors already in hand, compute the residual

\mathbf r = \mathbf b - A\hat{\mathbf x}

(ideally in higher precision), solve

A\mathbf d = \mathbf r

with the same factors, and update

\hat{\mathbf x} \leftarrow \hat{\mathbf x} + \mathbf d

. Each pass costs only

O(n^2)

. It recovers accuracy lost to instability, and with an extended-precision residual it can even beat

\kappa u

. Modern GPUs use it in reverse: factorise in fast low precision, then refine to double-precision accuracy.

Derivation

Why one correction can fix the answer

In exact arithmetic the correction is perfect. The error of

\hat{\mathbf x}

\mathbf e = \mathbf x - \hat{\mathbf x}

, and

A\mathbf e = A\mathbf x - A\hat{\mathbf x} = \mathbf b - A\hat{\mathbf x} = \mathbf r.

So the error itself solves a linear system with the same matrix and the residual as right-hand side, and

\hat{\mathbf x} + \mathbf e = \mathbf x

. In floating point

\mathbf d

is only an approximation to

\mathbf e

— computed with the same slightly wrong factors — but it is accurate to a few digits relative to

\mathbf e

, so each pass multiplies the error by a factor well below 1, as long as

\kappa(A)u

is small. The one thing that must be accurate is

\mathbf r

: it is a difference of nearly equal vectors (lesson 1-2), which is why it is computed in higher precision when possible.

💡Where you will use this

Structures: a structure with a nearly unrestrained mechanism, or stiffness values spanning ten orders of magnitude (steel beside rubber bearings), produces a huge

\kappa

; FE codes warn about 'small pivots' for exactly this reason. GNSS Navigation: poor satellite geometry makes the positioning equations ill-conditioned, and the dilution of precision (DOP) is the factor by which that geometry amplifies range errors into position errors. Attitude Control: Kalman filter covariance matrices drift towards indefinite in single precision unless the update is written in a stable form.

📚Further reading

Lloyd N. Trefethen and David Bau, Numerical Linear Algebra (SIAM, 1997), Lectures 21–22, 'Pivoting' and 'Stability of Gaussian Elimination' — including why growth factors are small in practice.
Nicholas Higham, Accuracy and Stability of Numerical Algorithms (2nd ed., SIAM, 2002), chapter 9, 'LU Factorization and Linear Equations', and chapter 12, 'Iterative Refinement'.
James H. Wilkinson, Rounding Errors in Algebraic Processes (Prentice-Hall, 1963; Dover reprint 1994) — backward error analysis, from the man who invented it.

↓But there is a problem

Elimination with pivoting is stable, its cost is

\tfrac23 n^3

, and conditioning tells us what accuracy to expect. For a dense matrix, that is the end of the story. But the stiffness matrix of a beam has only five nonzero diagonals. The matrix of a heated rod has three. Treating them as dense means doing

\tfrac23 n^3

operations to process about

5n

meaningful numbers. Before the First World War, a French artillery officer surveying Crete and North Africa found a way to halve the work for the symmetric matrices his survey produced. He was killed in 1918. His method was published by a colleague six years later.