NUM Numbers and algebra Algorithm
Gaussian elimination
Also known as: Gaussian elimination · LU decomposition · elimination method · row reduction · partial pivoting · pivoting · condition number · sparse matrices · fill-in · minimum degree · Kirchhoff's laws
If every street of Almaty is a wire and a battery is connected between two squares, the current spreads over hundreds of streets, and finding out which ones takes solving a system of 7908 equations. Gaussian elimination done head-on would spend a hundred and sixty billion operations on it, the same method in a good order a quarter of a million. Here is why, plus the pivot, conditioning and a factorisation done once.
Let every street of Almaty be a wire, with the wire’s resistance proportional to the length of the street. Connect a battery: plus at Republic Square, minus by the Green Bazaar. The current will not take one road but spread over the city, and how much of it each street gets is told by Kirchhoff’s laws: at every intersection as much flows out as flowed in. That is one linear equation per intersection, and the OpenStreetMap extract the catalogue uses has 7909 intersections. Below, this system is solved in your browser by Gaussian elimination, but in a special order. Move the battery’s poles by tapping the map.
Between Republic Square and the Green Bazaar the city behaves like a wire 222 metres long, although the shortest route along the streets is 3.8 km: hundreds of parallel roads share the current. Put a pole on Kok-Tobe and the resistance jumps to over three kilometres: there are few roads up the hill, and all the current goes along them. The bottom line is the price. Gaussian elimination that ignores zeros would spend 1.6·10¹¹ multiplications on this system, about a minute even in Go, not counting half a gigabyte of memory for the matrix itself. The same method, eliminating intersections in the “minimum degree” order, does a quarter of a million and fits into milliseconds, and each new pair of poles costs a fraction of a millisecond: the heavy part of the work is done once.
The idea
Express the first unknown from the first equation and substitute it into all the others: there is one unknown and one equation fewer. Repeat until one equation with one unknown is left; this is the forward pass, and it turns the system into a staircase. Then the backward pass: the last unknown is already known, substitute it into the second-to-last equation, and so on from the bottom up. Each elimination is subtracting a row multiplied by a number from the other rows; about n³/3 multiplications in all.
In detail and with proofs, why the moves do not change the solutions, when there are no solutions or infinitely many, and what the staircase looks like, it is all taken apart in the chapter on Gaussian elimination of “Mathematics, the Queen of the Sciences”. Here is the engineering side: where real systems come from, why zeros end up on the diagonal, how many digits the answer survives with, and how to solve systems of thousands and millions of equations where almost all the coefficients are zeros.
The equations of a circuit
Let φi be the potential of the i-th intersection and gij = 1/Rij the conductance of the street between i and j. By Ohm’s law the current along the street is gij(φi − φj), and by Kirchhoff’s law the sum of the currents flowing out of an intersection equals the current fed into it from outside: +1 at the battery’s plus, −1 at its minus, 0 everywhere else. Row i of the matrix has the sum of the conductances of all the intersection’s streets on the diagonal and −gij opposite each neighbour. Such a matrix is called the Laplacian of the graph, and it has an awkward property: it is singular.
If the street graph is connected, the Laplacian L has rank n − 1: the solutions of Lφ = 0 are exactly the constant vectors. And if the row and column of one node are deleted (it is “grounded”, φ = 0), the remaining matrix is non-singular.
For any vector φ, the sum of φi(Lφ)i over all nodes, with the brackets opened and the terms collected by streets, equals Σ gij(φi − φj)² over all streets. If Lφ = 0, this sum is zero, and all its terms are non-negative, so on every street φi = φj, and by connectivity the potential is the same over the whole city. Constant vectors really are solutions: in each row the coefficients add up to zero. So the solution space is one-dimensional, and the rank is n − 1.
Now let node z be grounded, and let the vector ψ on the other nodes solve the system without the row and column of z. Extend it by a zero at node z: the same sum Σ φi(Lφ)i, taken over the nodes other than z, still equals the sum of squares over the streets (the term of node z is zero because φz = 0), and so equals zero. The potential is again the same everywhere, and equal to zero, as at z. The only solution of the homogeneous system is zero, so the matrix is non-singular.
Physically this makes sense: voltage is a difference of potentials, and the potentials themselves are defined only up to a common constant until some node is declared “ground”. And if part of the circuit is not connected to the rest at all, its potential is set by nothing, and the system is singular even with a ground. Circuit simulators report a singular matrix in such a case and name the “floating” node.
Code: a Wheatstone bridge
The program solves a Wheatstone bridge, five resistors and a 9 V battery, by Gaussian elimination with partial pivoting, and then adds a resistor not connected to the rest of the circuit. The unknowns are the battery current and the potentials of three nodes; the fourth node is ground.
without row swaps: a zero pivot at step 1 with partial pivoting: potentials 9.0000 6.5138 4.1056 V battery current 47.1095 mA, the whole bridge is 191.044 ohm with a floating resistor 4-5: a zero pivot at step 6 of 6, the matrix is singular
The pivot
The first line of the output is not a typo. The unknowns are written like this: first the battery current, then the potentials, and the first equation, “φ₁ = 9 V”, does not contain the battery current. There is a zero on the diagonal, and you cannot divide by it, although the system solves perfectly well. This happens in any circuit with a voltage source if the equations are written in this order. The cure is partial pivoting: at step k, among the remaining rows take the one with the largest number in absolute value in column k, and swap it with the current one.
Zero is the extreme case. A small number is more dangerous: formally one can divide by it, but the multipliers come out huge, and subtracting a row multiplied by 10¹⁴ wipes out all but the leading digits in neighbouring rows. Pivoting keeps the multipliers no greater than one in absolute value, and for almost all matrices met in practice this is enough. Almost, because in 1961 James Wilkinson built matrices on which the numbers double at every step even with pivoting; in real life such matrices hardly ever turn up. Separately: for symmetric positive definite matrices, and the grounded Laplacian is one, swaps are not needed at all; that is why the rows in the first scene are not swapped, and the order is chosen only for the sake of the zeros.
How many digits the answer survives with
One can go wrong for two different reasons, and it is useful not to confuse them. Below, systems are solved whose answer is known in advance, all the unknowns equal to one, and you can see how many correct digits are left.
A matrix with a tiny number in the corner is a good matrix: with partial pivoting the answer is right to 13–15 digits at any size. Without swaps, zero to three of the same sixteen digits remain. This is the algorithm’s fault. But the Hilbert matrix, with elements 1/(i + j − 1), is bad in itself: its rows are almost proportional, and at n = 10 even with pivoting three digits are right, at n = 12 none. This is the problem’s fault: the slightest change in the right-hand side moves the answer by a huge amount, and no order of operations will fix that. The measure of this sensitivity was introduced by Alan Turing in 1948: the condition number κ(A) = ‖A‖·‖A⁻¹‖.
Let x be the exact solution of Ax = b, x̂ an approximate one, and r = b − Ax̂ its residual. Then
$$\frac{\|x - \hat x\|}{\|x\|} \le \kappa(A)\,\frac{\|r\|}{\|b\|}.$$A(x − x̂) = b − Ax̂ = r, whence x − x̂ = A⁻¹r and ‖x − x̂‖ ≤ ‖A⁻¹‖·‖r‖. On the other hand, ‖b‖ = ‖Ax‖ ≤ ‖A‖·‖x‖, that is, 1/‖x‖ ≤ ‖A‖/‖b‖. Multiplying the two inequalities gives the claim.
Gaussian elimination with partial pivoting almost always gives a residual at the level of rounding errors, about 10⁻¹⁶ of ‖b‖. The claim turns it into an error of the answer: up to κ·10⁻¹⁶. Hence the rule of numerical analysts: with κ ≈ 10k about k digits out of sixteen are lost. A small residual does not yet mean an accurate answer; it only means that the algorithm has done its part of the work.
n Hilbert: kappa digits no swaps tiny corner: digits no swaps 6 2.9e+07 9 9 15 0 8 3.4e+10 6 6 14 1 10 3.5e+13 3 3 15 1 12 4.1e+16 0 0 15 0
The condition number of the Hilbert matrix grows about 30 times with each new row, and the lost digits follow it almost exactly: κ = 3.4·10¹⁰, ten digits lost, six left. Here κ is computed in the row-sum norm, through the exact inverse of the Hilbert matrix, whose elements are integers.
A sparse matrix and fill-in
In the city’s matrix, 7909 × 7909, 62 million cells, only 28,687 are non-zero, 0.046%: an intersection has three or four neighbours. There is no point in storing and processing the zeros. But Gaussian elimination creates them: subtracting row k from row i, we write new numbers into all the columns where row k has non-zeros. This is called fill-in, and it has a simple meaning on the map.
Eliminating node k links every two not yet eliminated neighbours of k to each other, and only them: new non-zeros appear exactly in the places of pairs of neighbours of k that were not linked before.
When k is eliminated, the element aij is replaced by aij − aikakj/akk. The correction is non-zero only if aik ≠ 0 and akj ≠ 0, that is, i and j are both neighbours of k. If they were already neighbours, the non-zero just changes its value; if not, a new one appears. (A new non-zero could also accidentally become zero, but nobody counts on that.)
So the price of elimination depends on the order. Eliminate an out-of-the-way intersection with two neighbours and one new link appears. Eliminate a node in the middle that has forty neighbours by then, and they have to be linked in pairs: hundreds of new non-zeros and thousands of multiplications. See what this looks like on a 3 × 3 km piece of the centre.
In the original order, the order in which the intersections were met when parsing the OpenStreetMap extract, a red web of new links covers the whole district: 41,404 new non-zeros and 1.8 million multiplications. The reverse Cuthill–McKee order goes round the graph in layers, breadth first, and keeps the non-zeros near the diagonal in a narrow band: 7458 new ones and 62 thousand multiplications. The “minimum degree” order each time takes the node that has the fewest neighbours right now, and manages with 1944 new non-zeros and 8404 multiplications. For comparison, dense Gauss on the same 555 equations is 5.7·10⁷. Finding the order with the least fill-in is NP-hard (Mihalis Yannakakis proved this in 1981), so in practice heuristics like these two are used; minimum degree was described in 1967 by Tinney and Walker, power engineers: power systems have equally sparse systems, only the wires in them are real.
Factorise once
The moves of the forward pass are chosen by the matrix and do not depend on the right-hand side. If you record the multipliers by which rows were subtracted, you get a lower triangular matrix L, while the staircase is an upper triangular U, and A = LU (with a row permutation if rows were swapped). Solving Ax = b with a ready factorisation is two passes over triangular matrices, forward over L and back over U: O(n²) operations for a dense matrix and O(number of non-zeros in L and U) for a sparse one. That is exactly why moving a battery pole in the first scene costs almost nothing: the factorisation does not depend on where the current is fed.
The program below factorises the matrix of the whole city once and solves five systems, one for each of five places, and from these five solutions puts together the resistances between all pairs of places.
7909 intersections, 28687 nonzeros in the matrix (0.046%)
dense Gauss: about n^3/3 = 1.65e+11 multiply-adds
min degree: 36500 new nonzeros, 249163 multiply-adds
reverse Cuthill-McKee: 308386 new nonzeros, 7039994 multiply-adds
resistance between places, as a wire of this many km:
Republic Squ Green Bazaar Almaty-2 sta Botanical Ga Kok-Tobe
Republic Square 0.000 0.222 0.177 0.244 3.191
The resistance between two points of a graph is not a toy. It is also the “resistance distance”, which measures how well two places of a network are connected: the more independent paths, the smaller it is. It can also be obtained from random walks: the mean time for a random walker to get from one point to another and back is proportional to the resistance between them.
What it costs
At step k of the forward pass a row of length about n − k is subtracted from each of the n − k rows below: Σ(n − k)² ≈ n³/3 multiplications in all. The backward pass is about n²/2. Below is a benchmark: the same program with three nested loops in Python and in Go and, for comparison, numpy.linalg.solve, which calls LAPACK, the same factorisation with partial pivoting but written in blocks to suit the processor cache. The axes are logarithmic, and the slope of the line is the exponent.
The slopes are close to three: twice the equations, eight times as long. On hundreds of equations Python lags behind LAPACK almost a thousand times, Go 7–15 times: the difference is no longer the language but how the algorithm treats memory. And for the city’s system it is not loops that save the day but zeros.
| Matrix | Time | Memory | Example |
|---|---|---|---|
| dense | ~ n³/3 | n² | a thousand equations: a fraction of a second |
| banded, width w | ~ nw² | nw | tridiagonal (w = 1): O(n), the Thomas algorithm |
| sparse, a good order | by the fill-in | non-zeros of L and U | 7909 intersections: 2.5·10⁵ multiplications |
| with a ready LU | ~ n² or the non-zeros of L and U | — | each new right-hand side |
Where Gauss no longer does is where fill-in is inevitably large: three-dimensional meshes of the finite element method, where a node has dozens of neighbours and after any reordering the factorisation does not fit in memory. There systems are solved by iteration, with the conjugate gradient method and its relatives, which only multiply the matrix by a vector and create not a single new non-zero. Gaussian elimination stays inside them, as a way to build a “preconditioner”, an approximate incomplete factorisation that speeds up the iterations.
Where it works
- Scientific libraries.
numpy.linalg.solvecalls LAPACK’s gesv routine, LU with partial pivoting; on a singular matrix it answers with the error “Singular matrix”.scipy.sparse.linalg.spsolveis sparse LU, by default from the SuperLU library, with column reordering. - Electronics. Circuit simulators of the SPICE family solve a sparse system of Kirchhoff’s laws at every time step, of the same form as in the first scene, only with transistors.
- Power systems and engineering calculations. Network states, structural strength, heat conduction: sparse systems factorised in an order that reduces fill-in.
- Inside other methods. Implicit methods for stiff equations solve a system at every step (see the wall in the frost); least squares in linear regression comes down to a system of normal equations.
History
The oldest known description of the method is in the eighth chapter of “The Nine Chapters on the Mathematical Art”: parts of this book were written around 150 BC, and in the 3rd century Liu Hui provided it with commentaries. The first problem of the chapter is about sheaves of grain of three kinds, and the elimination is carried out on a counting board, in columns of rods. Later the method appears in Isaac Newton’s notes of 1669–1670, which complain that algebra textbooks do not teach how to solve simultaneous equations. In 1810 Gauss devised a notation for elimination in the symmetric systems of least squares, and professional computers used it for over a century; the method began to be called by Gauss’s name only in the 1950s, through a confusion in the history of the subject.
The machine age brought questions that computers with paper did not have. In 1948 Alan Turing wrote elimination as the factorisation A = LU and introduced the condition number. In 1961 James Wilkinson showed that the method with partial pivoting gives the exact answer to a very slightly changed problem, and so backward error analysis appeared. And in the 1960s engineers computing power systems learned not to touch the zeros, and sparse matrices became a science of their own.
Problems
In the first problem a circuit of resistors has to be computed; in one of the tests it is a piece of Almaty’s streets, and in another part of the circuit hangs in the air. In the second, Gaussian elimination works not with fractions but in the field of two elements: a puzzle where each press toggles a light and its neighbours.
The resistance of a circuit Lights outData: Almaty roads © OpenStreetMap contributors, ODbL