Numerical Zoo · Matrix

Hilbert matrix

Every entry looks harmless. The linear system does not.

Why numerical analysts keep this one around

The Hilbert matrix has entries H[i,j] = 1 / (i + j - 1). That is about as innocent a matrix formula as we could ask for. It is also famously ill conditioned as the dimension grows.

This makes it a useful antidote to the idea that difficult numerical problems have to look difficult. We can construct a right-hand side from a known solution, hand the resulting system to a solver, and then watch small perturbations in the data produce surprisingly large changes in the recovered coefficients.

The point is not that R cannot solve a Hilbert system. The point is that the numerical problem itself becomes sensitive. Once that happens, extra digits in the input may matter more than another clever line in the solver.

R setup

Start with the actual object.

n <- 8
H <- outer(1:n, 1:n, function(i, j) 1 / (i + j - 1))

x.true <- rep(1, n)
b <- as.vector(H %*% x.true)

x.hat <- solve(H, b)
cbind(x.true, x.hat, error = x.hat - x.true)

By hand first

Do enough arithmetic to see the trap.

For n = 2, the matrix is [[1, 1/2], [1/2, 1/3]]. Nothing in those entries announces a catastrophe.

Construct b = H * 1. In exact arithmetic, the solution is the vector of ones by definition.

Now perturb one component of b by a tiny amount and solve again. The difference in x is the behavior we care about.

What to look for

Now make it earn its reputation.

Increase n gradually instead of jumping straight to a large matrix. Conditioning is more instructive when we watch it deteriorate.
Compare the residual ||Hx - b|| with the actual error ||x - x.true||. A tiny residual does not guarantee that the recovered x is close to the known solution.
Try the iterative linear algebra lab if you want to see why 'the solver converged' and 'the problem was well conditioned' are separate statements.

Keep this distinction:

This is a good place to insist on two different questions: did the algorithm solve the equations it was given, and were those equations numerically capable of telling us the answer we wanted?

← Back to the Numerical Zoo