Solving Systems of Equations — Direct Methods and Iterative Methods

How much a bridge sags, how heat spreads through a room, what a Gaussian process predicts — once a computer gets hold of them they all turn into the same shape, Ax = b. This article builds up elimination (LU) and approximation (conjugate gradients) from zero, through why a million-unknown system can't be solved by elimination, all the way to condition numbers, preconditioning, and matrix-free solvers.

Working backwards from a web of relationships

How much does a bridge sag under load? How does temperature distribute across a room? What voltage sits at each node of a circuit? Different fields entirely — but the moment you hand any of them to a computer, they collapse into the same shape.

$$Ax = b$$

$A$ is a table (a matrix) of how the quantities relate to one another, $b$ is what the outside world imposes, and $x$ is the unknown you want. Put plainly: you know the mechanism $A$ and the outcome $b$, and you are working backwards to the cause $x$.

That is exactly what you did with simultaneous equations in school. The only difference is scale. A textbook problem has two or three unknowns and yields to substitution. Chop a bridge into a fine mesh and you have tens of thousands; a three-dimensional fluid gives you millions. What happens when you run the same elimination on a few million unknowns is where numerical computing begins.

And none of this is far from machine learning. Least squares, Gaussian process regression, Newton-style second-order optimization — every one of them has $Ax = b$ sitting at its core. Later in the article we'll look at exactly where the simulation toolbox and the ML toolbox part ways.

Two philosophies — eliminate, or approach

A direct method is one that produces the answer after a fixed, finite number of steps. It is school elimination scaled up: Gaussian elimination, and LU factorization, which stores that elimination for reuse. Ignore rounding error and what comes out is the exact solution.

An iterative method starts from a guess $x_0$ and nudges it toward the truth. It's like focusing a camera: every turn of the ring makes the image sharper, and you stop when it looks good enough.

The two behave in fundamentally different ways. A direct method gives you nothing if you stop halfway — a half-eliminated matrix is neither the solution nor an approximation of it, just an intermediate state. An iterative method always holds an answer, at every moment. You get to dial in accuracy against runtime after the fact.

If your simulation's discretization error already shows up in the third digit, solving the linear system to fifteen digits buys you nothing. That refusal to pay for accuracy you can't use is one big reason iterative methods took over practice.

Intuition — "solving" is the same as "finding the bottom of a valley"

Before the iterative machinery, it helps to have one picture in your head.

Suppose $A$ is symmetric ($a_{ij} = a_{ji}$) and positive definite ($v^\top A v > 0$ for every nonzero $v$). Consider this function:

Here $x^\top A x$ sums the products of $x$'s components weighted by $A$, and $b^\top x$ is the ordinary dot product of $b$ and $x$. With two variables, $f$ is a bowl.

The important part: the gradient (slope) of $f$ is $\nabla f(x) = Ax - b$, so the point where the slope vanishes — the bottom of the bowl — is precisely the solution of $Ax = b$. Solving a symmetric positive-definite system and finding the bottom of a bowl are the same job, which means the whole optimization toolbox transfers over intact.

But the bowl isn't necessarily round. When $A$'s eigenvalues are spread out, the bowl squashes into a long narrow valley, and steepest descent — step in whatever direction is locally steepest — becomes painfully inefficient. The slope points almost sideways, so "steepest downhill" means across the valley, and you make almost no progress along it.

That zigzag is the enemy the conjugate gradient method was built to kill. For optimization in general, see Convexity and Optimization.

Direct methods — Gaussian elimination and LU

Gaussian elimination repeats one move: use equation 1 to remove $x_1$ from every later equation, then use equation 2 to remove $x_2$, and so on until the system is triangular; then back-substitute from the bottom. LU factorization stores that elimination as $A = LU$ (lower triangular times upper triangular). Storing it pays off when you solve many problems with the same $A$ and different right-hand sides $b$: the factorization costs work proportional to $n$ cubed, but each subsequent solve costs only $n$ squared. If $A$ is symmetric positive definite you can use $A = LL^\top$ (Cholesky) instead, at roughly half the work. A Tour of Matrix Decompositions covers when to reach for which.

One rule of practice while we're here: never form the inverse. On paper we write $x = A^{-1}b$, but explicitly building $A^{-1}$ on a computer is both more expensive and less accurate. You always call a solve routine instead (np.linalg.solve in NumPy).

Sparsity, and the fill-in problem

Matrices that come out of simulation share a striking property: almost every entry is zero. The temperature at a grid point only relates directly to its neighbors, so a row typically holds five or seven nonzeros, not a million.

Store a million-by-million matrix densely and you need a trillion numbers — eight terabytes in double precision, which simply doesn't fit. Store only the seven nonzeros per row and you have seven million numbers, a few hundred megabytes including the index bookkeeping. Hopeless dense, laptop-sized sparse. Hence the specialized sparse storage formats.

The trouble starts when you eliminate. Running elimination on a sparse matrix creates nonzeros where zeros used to be: subtract one row from another and an empty cell fills in. This is called fill-in, and as it spreads, the sparse structure dissolves and both memory and arithmetic balloon.

The countermeasure is to reorder the eliminations (AMD, nested dissection, and friends), and for two-dimensional problems that works well enough to keep direct methods competitive. In three dimensions, no ordering saves you — fill-in dominates and direct methods hit a memory wall. Iterative methods never modify $A$ at all. All they ask for is "multiply $A$ by a vector," so the sparsity survives from beginning to end. That is the fork in the road.