Appearance
4.3 — Solving Systems: Gaussian Elimination
Chapter 2.4 solved two equations in two unknowns by substitution and elimination. Both worked. Neither scales — try it with twenty equations and you will lose track of what you have already used.
Real problems are large. A structural engineer solving for the forces in a truss has one equation per joint. A circuit simulator has one per node. A weather model has millions. Every one of them is solved by the same procedure, invented in China by 200 BCE, described in Europe by Gauss, and running right now in every numerical library on Earth.
1. Writing a system as a matrix equation
\begin{cases} 2x + y - z = 8 \\ -3x - y + 2z = -11 \\ -2x + y + 2z = -3\end{cases}
Separate the coefficients, the unknowns and the results:
\underbrace{\begin{bmatrix} 2 & 1 & -1 \\ -3 & -1 & 2 \\ -2 & 1 & 2\end{bmatrix}}_{A} \underbrace{\begin{bmatrix} x\\y\\z\end{bmatrix}}_{\mathbf{x}} = \underbrace{\begin{bmatrix} 8\\-11\\-3\end{bmatrix}}_{\mathbf{b}}
A\mathbf{x} = \mathbf{b}
Three characters, and they carry a real meaning. Chapter 4.2 said a matrix is a transformation. So A\mathbf{x} = \mathbf{b} asks: which input vector \mathbf{x} does the transformation A send to the known output \mathbf{b}? Solving a system is running a transformation backwards.
That reading immediately explains the three outcomes of Chapter 2.4:
- If A is invertible, exactly one input produces \mathbf{b}. One solution.
- If A squashes space (singular) and \mathbf{b} is not in the squashed region, nothing maps there. No solution.
- If A squashes space and \mathbf{b} is in the squashed region, a whole family maps there. Infinitely many solutions.
The geometry of Chapter 2.4 — lines crossing, parallel, or identical — is this in two dimensions.
2. Row operations and why they are safe
Write the system as an augmented matrix, coefficients and results together with a bar:
\left[\begin{array}{ccc|c} 2 & 1 & -1 & 8 \\ -3 & -1 & 2 & -11 \\ -2 & 1 & 2 & -3\end{array}\right]
Three operations are permitted, and each one is something you were already doing with equations:
- Swap two rows. Writing the equations in a different order changes nothing.
- Multiply a row by a nonzero number. Multiplying both sides of an equation by 3 changes nothing.
- Add a multiple of one row to another. This is elimination — the same move as Chapter 2.4.
Each is reversible, which is exactly why the solution set is unchanged. Nothing is lost and nothing is invented. Note the "nonzero" in rule 2: multiplying a row by zero would destroy an equation, which is the same trap as Chapter 2.1's warning about multiplying an equation by zero.
3. Forward elimination
The goal is row echelon form: zeros below the diagonal, so the last equation has one unknown, the second-last has two, and so on. Then you work backwards.
Start:
\left[\begin{array}{ccc|c} 2 & 1 & -1 & 8 \\ -3 & -1 & 2 & -11 \\ -2 & 1 & 2 & -3\end{array}\right]
Clear the first column below the top. The top-left entry, 2, is the pivot.
Row 2 becomes Row 2 + \frac{3}{2} Row 1: (-3 + 3, -1 + 1.5, 2 - 1.5, -11 + 12) = (0, 0.5, 0.5, 1).
Row 3 becomes Row 3 + Row 1: (-2+2, 1+1, 2-1, -3+8) = (0, 2, 1, 5).
\left[\begin{array}{ccc|c} 2 & 1 & -1 & 8 \\ 0 & 0.5 & 0.5 & 1 \\ 0 & 2 & 1 & 5\end{array}\right]
Clear the second column below the diagonal. The new pivot is 0.5.
Row 3 becomes Row 3 $- 4,Row 2: (0, 2-2, 1-2, 5-4) = (0,0,-1,1)$.
\left[\begin{array}{ccc|c} 2 & 1 & -1 & 8 \\ 0 & 0.5 & 0.5 & 1 \\ 0 & 0 & -1 & 1\end{array}\right]
Done. The matrix is triangular.
Back-substitution. Read the last row: -z = 1, so z = -1.
Second row: 0.5y + 0.5(-1) = 1, so 0.5y = 1.5, so y = 3.
First row: 2x + 3 - (-1) = 8, so 2x = 4, so x = 2.
Solution: (2, 3, -1). Substitute back into all three original equations to check: 4+3+1 = 8 ✓, -6-3-2 = -11 ✓, -4+3-2 = -3 ✓.
The procedure never required cleverness. That is the point — it is an algorithm, and a machine can run it on a matrix of any size.
4. When there is no solution, or too many
The algorithm tells you which case you are in, and you should be able to read the signal.
No solution shows up as a row like
\left[\begin{array}{ccc|c} 0 & 0 & 0 & 5\end{array}\right]
which reads 0x + 0y + 0z = 5, that is, 0 = 5. False, so the system is inconsistent and no assignment of values can satisfy it.
Infinitely many solutions show up as a row of all zeros, including the right-hand side:
\left[\begin{array}{ccc|c} 0 & 0 & 0 & 0\end{array}\right]
which reads 0 = 0 — true but empty. One equation was a combination of the others and told you nothing new.
When that happens, some unknowns become free variables: you may set them to anything, and the rest are then determined. Suppose elimination leaves
\left[\begin{array}{ccc|c} 1 & 2 & 1 & 5 \\ 0 & 1 & 1 & 3 \\ 0&0&0&0\end{array}\right]
Let z = t, any value. Then y = 3 - t, and x = 5 - 2(3-t) - t = -1 + t. The solution set is a line in three-dimensional space:
(x,y,z) = (-1, 3, 0) + t(1,-1,1)
a point plus a direction. That direction vector is precisely the null space of Chapter 4.2 — the set of things A sends to zero, which you can add freely without changing the output.
The rank test. The rank of A is the number of nonzero rows after elimination.
- Rank of A equals rank of the augmented matrix, and equals the number of unknowns → one solution.
- Rank of A equals rank of the augmented matrix, but is less than the number of unknowns → infinitely many, with (unknowns − rank) free variables.
- Rank of A is less than the rank of the augmented matrix → no solution.
5. What can go wrong on a computer
The algorithm is exact in arithmetic and troublesome in floating point, and the reasons are worth knowing because they explain a whole class of real bugs.
A zero pivot. If the pivot position holds a zero, you cannot divide by it. The fix is to swap in a row below that has a nonzero entry there.
A tiny pivot, which is worse. Dividing by 10^{-15} multiplies every subsequent number by 10^{15}, and the rounding errors in those numbers grow just as much. The answer comes out looking like a number and being nonsense.
The standard defence is partial pivoting: at each step, swap in the row whose entry in that column has the largest absolute value. It costs almost nothing and it is what every serious library does by default.
Ill-conditioning, which no algorithm can fix. Consider:
\begin{cases} x + y = 2 \\ x + 1.0001y = 2.0001\end{cases}
The solution is x = y = 1. Now change the last number to 2.0002 — a change of one part in twenty thousand. The solution becomes x = 0, y = 2. A tiny change in the input produced a total change in the answer.
Geometrically the two lines are nearly parallel, so their crossing point slides enormously when either line twitches. The condition number measures this, and a large one means your answer is only as trustworthy as your inputs, no matter how good your arithmetic. Chapter 10.1 develops it properly.
This is a genuine engineering hazard, not a textbook curiosity. It is why structural analyses are checked for conditioning, why a regression with nearly-dependent columns produces wildly unstable coefficients, and why "the code ran and produced numbers" is not the same as "the answer is right".
6. The cost, and why direct solving stops being an option
Gaussian elimination on an n \times n system takes roughly \frac{2}{3}n^3 arithmetic operations. Cubic growth, and Chapter 1.4 warned about what that means.
| Size n | Operations, roughly |
|---|---|
| 100 | 700 thousand |
| 1,000 | 700 million |
| 10,000 | 700 billion |
| 1,000,000 | 7 \times 10^{17} |
Ten times bigger is a thousand times slower. A million-unknown system — perfectly ordinary in weather modelling or structural simulation — is out of reach by direct elimination on any machine.
Two things save it.
Sparsity. Real large systems are almost all zeros. In a structural model each joint connects to a handful of neighbours, not to all million. Specialised sparse solvers skip the zeros and reduce the cost enormously.
Iterative methods. Rather than solving exactly, start from a guess and improve it repeatedly until it is good enough. Each step is one cheap matrix-vector multiply. Methods such as conjugate gradient converge fast for well-behaved systems, and Chapter 10.2 covers the idea. Approximately right in a minute beats exactly right in a year, and for most physical problems the input data is not exact anyway.
7. Least squares: solving systems that have no solution
Here is a situation that sounds contradictory and is enormously important.
You have 500 data points and want to fit a straight line through them. That is 500 equations in 2 unknowns (slope and intercept). There is no solution — no line passes through 500 scattered points.
So change the question. Instead of "which \mathbf{x} makes A\mathbf{x} = \mathbf{b}", ask "which \mathbf{x} makes A\mathbf{x} as close as possible to \mathbf{b}?" Measure closeness by the length of the error vector \|A\mathbf{x} - \mathbf{b}\|.
The answer is one formula, the normal equations:
A^\mathsf{T}A\,\mathbf{x} = A^\mathsf{T}\mathbf{b}
which is a small square system you solve by elimination.
Why this is the right answer, geometrically. The reachable outputs A\mathbf{x} form a plane (the column space of Chapter 4.2), and \mathbf{b} is a point off that plane. The closest point on a plane to an external point is the perpendicular foot — drop a perpendicular. "Perpendicular" means the error is orthogonal to every column of A, and by Chapter 4.1 that means every dot product is zero, which is exactly A^\mathsf{T}(A\mathbf{x} - \mathbf{b}) = 0. Rearrange and you have the normal equations.
This is linear regression — Chapter 7.8 approaches it from the statistics side and arrives at the same formula. It is behind every trend line, every calibration curve, every GPS position fix (which has more satellites than unknowns and averages out the inconsistencies), and the final layer of many machine learning models.
Every formula above, built from scratch
None of the results in this chapter are worth memorising, because each one can be rebuilt in under a minute from something simpler. What follows is that rebuilding, one result at a time, so the formula and the reason for it sit on the same page as the explanation that needed them.
Systems of equations
A\mathbf{x} = \mathbf{b}
| \det A | Rank condition | Solutions |
|---|---|---|
| \ne 0 | full rank | exactly one |
| =0 | equations agree | infinitely many |
| =0 | equations contradict | none |
The middle and last rows are told apart by rank: write the coefficients and the right-hand side side by side as one wide matrix, and compare its rank with the rank of the coefficients alone. Equal ranks mean the extra column said nothing new, so the system is consistent and has a whole family of solutions. A larger rank for the wide matrix means the right-hand side demanded something the columns cannot produce, and there is no solution at all.
The rank is the number of genuinely independent rows — how many of the equations say something new. Two equations describing the same line have rank 1, however they are written.
Gaussian elimination, in three moves: swap two rows, multiply a row by a non-zero number, add a multiple of one row to another. None of these changes the solution set, because each is reversible and each is a legal operation on equations. Repeat until the matrix is triangular, then substitute backwards.
The cost. Elimination on an n\times n system takes about \frac{2}{3}n^3 arithmetic operations. Cramer's rule, done naively, takes about (n+1)! — for n = 20 that is the difference between an instant and longer than the age of the universe.
8. Where this shows up in your life
Every structural engineering calculation. Forces at joints, one equation per joint, solved by elimination.
Every circuit simulation. Kirchhoff's current law gives one equation per node; SPICE and every EDA tool solve exactly this. Volume III, 1.3 covers nodal analysis.
Every weather forecast and crash simulation. Millions of unknowns, sparse and iterative.
Every trend line and calibration. Least squares.
Every GPS fix. More satellites than unknowns, solved in the least-squares sense so the noisy measurements average out.
Every economic input-output model. Leontief's model of an economy is literally \mathbf{x} = A\mathbf{x} + \mathbf{d}, rearranged to (I - A)\mathbf{x} = \mathbf{d} and solved by elimination. It won the 1973 Nobel Prize in Economics.
Elimination tells you whether a system has a unique solution, but only by running it. There is a single number that answers the question in advance, and it also measures exactly how much a transformation stretches space. That is the determinant.