Gaussian elimination

In the last chapter a matrix became a transformation: $A$ takes a vector $\mathbf x$ and carries it to $A\mathbf x$. Solving the system $A\mathbf x = \mathbf b$ means travelling that road backwards and finding out which vector the transformation sent to $\mathbf b$. For two unknowns the formula for the inverse matrix is enough. But an astronomer looking for a planet’s orbit has six unknowns, and an engineer designing the frame of a bridge has hundreds of thousands. We need a method that doesn’t care how many equations there are and never fails.

Such a method was known long before matrices, and what made it famous was a disappearance. On the first night of the nineteenth century, 1 January 1801, Giuseppe Piazzi at the Palermo observatory noticed a faint star that wasn’t in the catalogue. The next night it had moved. Piazzi took it for a comet, although for a comet it crept too slowly and too evenly. He observed it 24 times until he fell ill on 11 February, and then the new body, later named Ceres, vanished into the glare of the Sun. To find it again in the autumn an orbit was needed, and six weeks of observations gave only a tiny piece of its path.

An orbit is given by six numbers, the orbital elements: the size and elongation of the ellipse, the tilt of its plane and the line along which that plane cuts the plane of the Earth’s orbit, the turn of the ellipse within its plane, and the planet’s position at a chosen moment. The job of recovering them from Piazzi’s observations was taken on by the twenty-four-year-old Carl Friedrich Gauss. How the hunt ended we’ll tell at the end of the chapter. First, the tool Gauss couldn’t have managed without.

When refining an orbit from many observations, a computer (in those days a person who did calculations) takes approximate elements and looks for small corrections to them. For small corrections the equations become linear, just as the tangent from the chapter on the derivative replaces a curve near a point. The result is a system of linear equations in six unknowns, which has to be solved by hand, again and again. Gauss solved such systems by eliminating the unknowns one after another. Let’s see how this method works, why it always works, and what it costs.

Three permitted moves

Let’s start with a toy version of the astronomer’s problem. A point was observed three times: at the moments $t = 1$, $2$ and $3$ it was at heights $6$, $11$ and $18$. Suppose the height changes by the law $y = a + bt + ct^2$, like a thrown stone. Each observation gives an equation in the unknown coefficients:

$$\left\{\begin{aligned} a + b + c &= 6, \\ a + 2b + 4c &= 11, \\ a + 3b + 9c &= 18. \end{aligned}\right.$$

Subtract the first equation from the second and from the third. The unknown $a$ disappears from both: we’re left with $b + 3c = 5$ and $2b + 8c = 12$. From the last, subtract twice the one before and we get $2c = 2$. Now we climb back up: $c = 1$, then $b = 5 - 3c = 2$, then $a = 6 - b - c = 3$. The law of motion is $y = 3 + 2t + t^2$, and at the moment $t = 4$ the point will be at height $3 + 8 + 16 = 27$. We’ve predicted the next observation, a miniature of what Gauss had to do.

The letters in these calculations only take up space: all the work is done by the coefficients. Let’s write them out as a table, with a bar separating the right-hand sides. Each operation on the equations becomes an operation on the rows of the table:

$$\left(\begin{array}{ccc|c} 1 & 1 & 1 & 6 \\ 1 & 2 & 4 & 11 \\ 1 & 3 & 9 & 18 \end{array}\right) \xrightarrow{\substack{R_2 - R_1 \\ R_3 - R_1}} \left(\begin{array}{ccc|c} 1 & 1 & 1 & 6 \\ 0 & 1 & 3 & 5 \\ 0 & 2 & 8 & 12 \end{array}\right) \xrightarrow{R_3 - 2R_2} \left(\begin{array}{ccc|c} 1 & 1 & 1 & 6 \\ 0 & 1 & 3 & 5 \\ 0 & 0 & 2 & 2 \end{array}\right).$$

On equations the moves do what we did by hand: swap equations, multiply an equation by a number, add one equation to another. But why, after such moves, are no solutions lost and no extra ones gained? This needs proving, and the proof is the most visual part of the whole method.

The trick is that each move can be undone by another move. So it’s enough to check one thing: after a move, the old solutions remain solutions.

Let’s start with two equations in two unknowns. In the picture they’re two lines, $R_1$ and $R_2$ (at first $x + y = 3$ and $x - y = 1$), and their common point $P$ is the solution of the system. Replace the second equation with $R_2 + kR_1$. The point $P$ turns both equations into true equalities. Add these equalities, multiplying the first by $k$, and we get a true equality, so $P$ satisfies the new equation too. For every $k$ the new line passes through $P$: by changing $k$ we rotate the second line about the solution. Such a family of lines through one point is called a pencil. Choose $k$ so that $x$ disappears from the equation (at first that’s $k = -1$: $(x - y) - (x + y) = 1 - 3$, that is, $-2y = -2$). The line has become horizontal and still passes through $P$. That’s what eliminating an unknown looks like in a picture: rotating a line about the solution until it lies parallel to an axis. In a system of any number of equations the argument is the same. The move $R_i \to R_i + kR_j$ changes only equation $i$, and if $\mathbf x$ satisfied equations $i$ and $j$, it satisfies the new one too: the left side of the new equation is the sum of the left sides with factors $1$ and $k$, and the right side the sum of the right sides with the same factors. Multiplying an equation by a number and swapping equations preserve solutions all the more. So every solution of the old system is a solution of the new one. Conversely: the move $R_i \to R_i + kR_j$ is undone by the move $R_i \to R_i - kR_j$ (the row $R_j$ hasn’t changed), multiplication by $k \ne 0$ by multiplication by $\frac1k$, a swap by the same swap. Apply the undoing move to the new system: by what’s already proved, its solutions remain solutions of the old one. Both systems have the same solutions. Carry the elimination through to the end, and the lines become vertical and horizontal, $x = 2$ and $y = 1$: the answer can be read off at once. It’s also clear why a row mustn’t be multiplied by zero. Such a move can’t be undone, and every point satisfies the resulting equation $0 = 0$.

A system has two equations: $x + y = 3$ and $x - y = 1$. The second equation was multiplied by $0$. What happened to the set of solutions?

It didn’t change: multiplication is a permitted move It became infinite: one line, $x + y = 3$, is left It became empty: $0 = 1$ Empty, because both unknowns disappeared

The equation $0x + 0y = 0$ holds for all $x$ and $y$, so what’s left of the system is the single equation $x + y = 3$, a whole line of solutions. You may multiply only by a non-zero number: such a move is reversible, and the solutions don’t change.

The staircase

We have the moves; we need a strategy. The goal is a matrix with nothing but zeros below a “staircase”, where each next step begins further right than the one before:

$$\left(\begin{array}{cccc|c} \boxed{2} & 1 & 0 & 3 & 5 \\ 0 & 0 & \boxed{-1} & 4 & 1 \\ 0 & 0 & 0 & \boxed{7} & 14 \\ 0 & 0 & 0 & 0 & 0 \end{array}\right)$$

A system in echelon form is solved from the bottom up, as we’ve already done: the last non-zero equation contains the fewest unknowns, so we solve it first and substitute what we find into the ones above. This is back substitution. The staircase itself is obtained by forward elimination. We go through the columns from left to right. In the current column we look for a non-zero entry among the rows that haven’t yet become steps; if there’s none, we skip the column. We bring the entry found to the top of the remaining part by a swap: it’s the new pivot. From each row below we subtract the pivot row with a factor chosen to make a zero under the pivot: if the pivot is $p$ and the entry under it is $q$, the factor is $\frac qp$. Then we go down a row and move on to the next column.

We argue by induction on the number of columns: forward elimination deals with the first column and then works with a smaller matrix.

If there’s one column, the procedure does this: if the column has a non-zero entry, it brings it to the top and zeros everything below it; if all the entries are zero, it does nothing. In both cases we get a column in echelon form. Suppose the statement is proved for matrices with fewer columns. If the first column is entirely zero, set it aside: the procedure applies to the rest of the matrix in the same way, and by the inductive hypothesis it becomes echelon. A zero column on the left doesn’t spoil the staircase: all the steps are simply shifted one column to the right. If the first column has a non-zero entry $p$, a swap brings it to the first row, and subtractions turn all the entries under it into zeros: from the row with entry $q$ we subtract the first row multiplied by $\frac qp$, and in place of $q$ we’re left with $q - \frac qp \cdot p = 0$. After that the procedure touches neither the first row nor the first column (which has zeros under $p$), and works with the matrix made of the remaining rows without the first column. By the inductive hypothesis it becomes echelon. At the top is a row with the pivot $p$ in the first column, and below it an echelon matrix whose first column is zero. The steps of the lower part begin to the right of the first column, so the whole matrix is in echelon form.

The staircase can be polished. Divide each pivot row by its pivot so that the steps become ones, and zero everything above them by working back through the columns. For our parabola system we get

$$\left(\begin{array}{ccc|c} 1 & 1 & 1 & 6 \\ 0 & 1 & 3 & 5 \\ 0 & 0 & 2 & 2 \end{array}\right) \longrightarrow \left(\begin{array}{ccc|c} 1 & 0 & 0 & 3 \\ 0 & 1 & 0 & 2 \\ 0 & 0 & 1 & 1 \end{array}\right),$$

and the answer $a = 3$, $b = 2$, $c = 1$ stands right there in the last column.

It’s named after Wilhelm Jordan, a German geodesist who described this variant in a surveying textbook at the end of the nineteenth century. The method is often mistakenly linked with Camille Jordan, whose normal form of matrices we’ll meet in chapter 38.

Sometimes a zero turns up in the place of the next pivot. Take the system $x + y + 2z = 5$, $2x + 2y + 5z = 11$, $3x + 4y + 7z = 17$. After the subtractions $R_2 - 2R_1$ and $R_3 - 3R_1$ both $x$ and $y$ have disappeared from the second equation: it’s become $z = 1$. The pivot for the second column has to be looked for lower down: it’s in the third row, $y + z = 2$. Swap the rows, and the staircase is ready. Try doing it yourself.

Solve the system $x + y + z = 4$, $2x + 3y + z = 4$, $x - y + 2z = 9$. Write the answer as a triple $(x, y, z)$.

$R_2 - 2R_1$ gives $y - z = -4$, and $R_3 - R_1$ gives $-2y + z = 5$. Add twice the second to the third: $-z = -3$, that is, $z = 3$. Then $y = -4 + z = -1$ and $x = 4 - y - z = 2$. A check with the third equation: $2 + 1 + 6 = 9$.

One, infinitely many or none

So far the staircase has come out complete: as many steps as unknowns. It isn’t always so. Take the system

$$\left\{\begin{aligned} x + 2y + z &= 4, \\ 2x + 4y + 3z &= 9, \\ 3x + 6y + 4z &= 13. \end{aligned}\right.$$

The subtractions $R_2 - 2R_1$ and $R_3 - 3R_1$ give identical rows $z = 1$, and after $R_3 - R_2$ the third equation turns into $0 = 0$. The backward pass leads to reduced row echelon form:

$$\left(\begin{array}{ccc|c} 1 & 2 & 0 & 3 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 \end{array}\right), \qquad\text{that is,}\qquad x + 2y = 3, \quad z = 1.$$

A step in the second column never appeared, and the unknown $y$ isn’t constrained by anything. Set it however you like, and the others adjust: $x = 3 - 2y$, $z = 1$. There are infinitely many solutions, and they all lie on one line in space.

Any solution of the system: a column $(x, y, z)$. A particular solution: the free unknown was given the value $0$, and the leading ones were found by back substitution. The parameter: the value of the free unknown $y$. It runs through all real numbers. The direction of the line of solutions. This column solves the homogeneous system with zero right-hand sides: moving along it doesn’t change the left sides of the equations. Example: at $t = 2$ we get $(-1, 2, 1)$. Check: $-1 + 4 + 1 = 4$, $-2 + 8 + 3 = 9$, $-3 + 12 + 4 = 13$.

If there are two free unknowns, the general solution has two parameters, and the solutions fill a plane. And if the right side of the third equation is $14$ instead of $13$, the same chain of moves ends with the row $0 = 1$. This equation is false for every $x$, $y$, $z$, and the system is inconsistent.

Everything is decided by the number of steps. It can only be counted once the staircase has been built, and a fair question arises: doesn’t it depend on which moves we made? Different people choose pivots differently and get different echelon matrices. It turns out that after polishing to the reduced form, all roads lead to the same place.

The idea: the reduced form can be read from the relations between the columns of the matrix, and the moves don’t change these relations.

Write the columns of $A$ as $\mathbf a_1, \dots, \mathbf a_n$. The equation $c_1\mathbf a_1 + \dots + c_n\mathbf a_n = \mathbf 0$ means exactly that the column $(c_1, \dots, c_n)$ solves the homogeneous system $A\mathbf x = \mathbf 0$: the $i$th row of this vector equation is the $i$th equation of the system. If the matrix $B$ is obtained from $A$ by elementary row operations, then by the previous theorem the systems $A\mathbf x = \mathbf 0$ and $B\mathbf x = \mathbf 0$ have the same solutions. So any linear combination of the columns of $A$ is zero if and only if the same combination of the columns of $B$ is zero. In particular, a column of $A$ is expressed through the preceding columns with some coefficients exactly when the same column of $B$ is expressed through the preceding columns of $B$ with the same coefficients. Look at the reduced row echelon form $R$. The column with the $k$th pivot is the $k$th unit column $\mathbf e_k$: a one in row $k$, zeros elsewhere. It isn’t expressed through the preceding columns, because they all have zeros in row $k$. A column without a pivot has non-zero entries only in the rows of pivots already met, so it equals $r_1\mathbf e_1 + \dots + r_k\mathbf e_k$, where $r_1, \dots, r_k$ are its entries, and it’s expressed through the preceding pivot columns. So the pivot columns of $R$ are exactly those not expressed through the preceding ones. By the second step a column of $R$ has this property exactly when the column of $A$ with the same number does, so the positions of the steps are determined by the matrix $A$ itself. A column of $R$ without a pivot equals $r_1\mathbf e_1 + \dots + r_k\mathbf e_k$, and this expression is unique: the coefficient of $\mathbf e_i$ is simply the $i$th entry of the column. By the second step the corresponding column of $A$ is expressed through the pivot columns of $A$ with the same coefficients $r_1, \dots, r_k$. If it could be expressed through the pivot columns of $A$ in two ways, it could also be expressed through $\mathbf e_1, \dots, \mathbf e_k$ in two ways, which is impossible. So the numbers $r_i$ are determined by the matrix $A$. The positions of the leading ones and all the other entries of the reduced form are determined by the matrix $A$, not by the choice of moves. Two reduced forms of the same matrix coincide.

Any echelon form can be polished into the reduced one without changing the number of steps, so the rank can be counted from any staircase. The parabola system has $\operatorname{rank} A = 3$; the system with the free $y$ has rank $2$. Now we can say when a system has solutions by comparing two numbers.

The trick is that forward elimination doesn’t distinguish the last column from the others: everything depends on whether a step lands there.

Bring $(A \mid \mathbf b)$ to echelon form. Without the last column we get an echelon form of $A$ itself: all the moves were chosen column by column from left to right. So the ranks differ exactly when a pivot lands in the last column, that is, when a row $(0\ \dots\ 0 \mid c)$ with $c \ne 0$ appears. In the picture there are three equations in two unknowns: three lines. The first two meet at the point $P$ and give steps in the columns of $x$ and $y$. Forward elimination subtracts from the third equation the combination $\alpha R_1 + \beta R_2$ with the same left side. All the combinations of $R_1$ and $R_2$ are lines of the pencil through $P$, and the one we need is the line through $P$ parallel to the third. If the third line misses $P$, what’s left of the third equation is $0 \cdot x + 0 \cdot y = c$ with $c \ne 0$. No point satisfies this equation: there are no solutions, and $\operatorname{rank}\,(A \mid \mathbf b) = 3 > 2 = \operatorname{rank} A$. The further the third line is from $P$, the larger $|c|$. If there’s no pivot in the last column, all the steps stand in the columns of unknowns. Give the free unknowns any values and find the leading ones by back substitution: each non-zero equation expresses its leading unknown through those further right. A solution exists. Drag the third line through $P$: the row becomes entirely zero, $0 = 0$, and $P$ becomes a solution of the whole system. The system is consistent exactly when the ranks are equal. If then $r = n$, there are no free unknowns, back substitution determines each unknown uniquely, and there’s one solution. If $r < n$, there are $n - r$ free unknowns; each takes any value, and different values give different solutions, infinitely many of them.

This theorem goes by different names in different books. In English-language books it’s the Rouché–Capelli theorem; others name it after Leopold Kronecker and Alfredo Capelli, or call it the Rouché–Fontené or the Rouché–Frobenius theorem.

In space each equation in three unknowns gives a plane, and the theorem turns into a catalogue of the ways three planes can be arranged:

The “book” row also includes the case where two planes coincide and the third cuts them. A “prism” is three planes that meet in pairs along three parallel lines; the other inconsistent case of the same rank is two parallel planes and a third that cuts both. In the last row all three are parallel, and at least two of them are distinct.

An especially simple case is a homogeneous system, where all the right-hand sides are zero. It’s always consistent: the zero column solves it. The interesting question is whether it has other solutions.

This modest corollary will do a job in the next chapter: the proof that a space has a dimension rests on it.

A system has three equations and four unknowns. Can it have exactly one solution?

Yes, if the equations are independent No: it has either no solutions or infinitely many No: it always has infinitely many solutions

The rank is at most the number of rows, $3 < 4$, so there’s always a free unknown. If the system is consistent, it takes any value. But it may not be consistent: for example, if two equations differ only in their right-hand sides.

Find the rank of the matrix $\begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix}$.

$R_2 - 4R_1$ and $R_3 - 7R_1$ give the rows $(0, -3, -6)$ and $(0, -6, -12)$. The second of these is twice the first, and after $R_3 - 2R_2$ a zero row is left. There are two steps: the rank is $2$.

For which $a$ does the system $x + y + z = 1$, $x + 2y + 3z = 2$, $2x + 3y + (a + 6)z = 4$ have no solutions?

$R_2 - R_1$ gives $y + 2z = 1$, and $R_3 - 2R_1$ gives $y + (a + 4)z = 2$. Subtract one from the other: $(a + 2)z = 1$. For $a \ne -2$ this gives $z$, and the solution is unique. For $a = -2$ we get $0 = 1$: the rank of the matrix is $2$, the rank of the augmented matrix $3$, and there are no solutions.

What a solution costs

Gauss computed by hand, and every multiplication took time. How many operations does Gaussian elimination need to solve a system of $n$ equations in $n$ unknowns?

The trick is to see where the work is concentrated: in squares that shrink with each step.

The first step of forward elimination. There are $n - 1$ rows under the pivot. Each needs a factor (one division) and an update of its other $n$ numbers: $n - 1$ entries of the matrix and the right-hand side. Each update is one multiplication and one subtraction. The main work goes on in the $(n - 1) \times (n - 1)$ square in the bottom right corner of the matrix. The second step works with an $(n - 2) \times (n - 2)$ square, the third with an even smaller one. The squares are nested inside each other and shrink towards the corner. Altogether they contain $(n - 1)^2 + (n - 2)^2 + \dots + 1^2$ cells, and each cell is a multiplication plus a subtraction. Stand these numbers up as bars of width $1$: the bar of height $m^2$ over the interval $[m, m + 1]$ lies entirely under the parabola $y = x^2$, while the same bar over the interval $[m - 1, m]$ covers it completely. So the sum of the bars is squeezed between areas under the parabola: $\frac{(n - 1)^3}{3} \le 1^2 + 2^2 + \dots + (n - 1)^2 \le \frac{n^3}{3}$. We computed the area under the parabola in the chapter on the integral. The exact count. At the step with $m$ rows under the pivot, $m(m + 2)$ multiplications and divisions and $m(m + 1)$ subtractions are done. Back substitution: unknown number $i$ costs $n - i$ multiplications, as many subtractions and one division. Adding up with the formulas $1 + 2 + \dots + k = \frac{k(k + 1)}{2}$ and $1^2 + \dots + k^2 = \frac{k(k + 1)(2k + 1)}{6}$ from the chapter on sequences, we get $\sum_{m=1}^{n-1} m(m + 2) + \frac{n(n + 1)}{2} = \frac{n^3}{3} + n^2 - \frac n3$ and $\sum_{m=1}^{n-1} m(m + 1) + \frac{n(n - 1)}{2} = \frac{n^3}{3} + \frac{n^2}{2} - \frac{5n}{6}$. Everything apart from the squares, the factors, the right-hand sides, back substitution, adds only terms of order $n^2$. For large $n$ the main part is $\frac{2n^3}{3}$ operations: twice as many unknowns, eight times as much work. The main part: updating the shrinking squares during forward elimination, one multiplication and one subtraction per cell. The factors, the right-hand sides and back substitution: a few operations for each entry of the matrix. A correction left over from the exact sums; for large $n$ it’s negligible. Example: at $n = 3$ we get $18 + 13.5 - 3.5 = 28$ operations ($17$ multiplications and divisions and $11$ additions and subtractions). For Gauss’s six unknowns it’s $191$ operations; for a thousand, $668{,}165{,}500$, of which $666{,}666{,}667$ come from the first term.

Is that a lot or a little? Compare it with Cramer’s rule from the chapter on systems: each unknown is a ratio of two determinants. If a determinant is computed by expanding along a row, it has $n!$ terms. For $n = 20$ that’s about $10^{21}$ multiplications: at a billion operations a second, thirty thousand years. Gaussian elimination solves the same system in $5910$ operations, that is, in a few microseconds.

Cubic growth is the main law of this section. Double the size of the mesh in an engineering model, and the work grows eightfold. That’s why, for huge systems where almost all the coefficients are zero, many workarounds have been devised that exploit those zeros. But a dense system with no special structure is still solved by Gaussian elimination. It’s also how the power of supercomputers is measured: in the TOP500 ranking, kept since 1993, the machines compete in a single exercise, solving a huge dense system of linear equations.

How many multiplications and divisions will Gaussian elimination need for a system of $10$ equations in $10$ unknowns (with a unique solution and no swaps)?

By the formula $\frac{n^3}{3} + n^2 - \frac n3$ at $n = 10$: $\frac{1000}{3} + 100 - \frac{10}{3} = \frac{990}{3} + 100 = 330 + 100 = 430$.

Thin ice

A computer stores numbers with a finite number of significant digits, usually about sixteen. After each operation the result is rounded, and the errors accumulate. For Gaussian elimination this is no trifle: an unlucky pivot can ruin the whole answer. Let’s show this on a calculator that remembers only three significant digits. Solve the system

$$\left\{\begin{aligned} 0.0001x + y &= 1, \\ x + y &= 2. \end{aligned}\right.$$

Its exact solution is $x = \frac{10{,}000}{9999} \approx 1.0001$, $y = \frac{9998}{9999} \approx 0.9999$. The pivot is $0.0001$, the factor $10{,}000$. Subtract: $(1 - 10{,}000)y = 2 - 10{,}000$, that is, $-9999y = -9998$. The calculator’s three digits turn both sides into $-10{,}000$, and we get $y = 1$. Substitute into the first equation: $0.0001x = 1 - 1 = 0$, that is, $x = 0$. The answer is completely wrong. The two from the second equation drowned in the huge $10{,}000$, and then the tiny pivot blew up the remaining error.

Swap the equations. The pivot is now $1$, the factor $0.0001$: $(1 - 0.0001)y = 1 - 0.0002$, after rounding $y = 1$, and $x = 2 - y = 1$. The answer is right to three digits.

Pivoting fixes the errors of the algorithm, but it can’t fix the problem itself. Recall the system from the chapter on systems: $x + y = 2$, $x + 1.01y = 2.01$, with solution $(1, 1)$. Shift the right side by a hundredth, to $2.02$, and the solution moves to $(0, 2)$. The right side changed by half a per cent, and the solution by a hundred per cent. The lines are almost parallel, and no order of operations can fix that: that’s how the data themselves are.

For the matrix in our example $\kappa \approx 402$, and half a per cent turned into a hundred per cent, two hundred times more, within the promised bound. A rough rule of numerical analysts: when $\kappa \approx 10^k$, about $k$ correct digits are lost in the answer. For the $10 \times 10$ Hilbert matrix, with entries $\frac{1}{i + j - 1}$, the condition number is about $1.6 \cdot 10^{13}$, and of the computer’s sixteen digits only two or three will be correct.

Recording the moves: LU decomposition

An astronomer gets new observations, but the moments at which they’re made are the same. In our toy this means the matrix $A$ is the same, and only the right-hand sides change. The moves of forward elimination are chosen by the matrix and don’t depend on the right-hand sides. They’re worth recording once.

$L$ (for lower): the log of moves. Rows $2$ and $3$ had $R_1$ subtracted with factor $1$, and then $R_2$ was subtracted from row $3$ with factor $2$. $U$ (for upper): the staircase with which forward elimination ended in the parabola problem. Example: new observations $2$, $4$, $8$ at the same moments. Solve $L\mathbf z = \mathbf b$ from the top down: $z_1 = 2$, $z_2 = 4 - 2 = 2$, $z_3 = 8 - 2 - 2 \cdot 2 = 2$. Then $U\mathbf x = \mathbf z$ from the bottom up: $c = 1$, $b = 2 - 3 = -1$, $a = 2 + 1 - 1 = 2$. The law $y = 2 - t + t^2$ passes through all three points.

Two triangular systems are solved in about $2n^2$ operations, against $\frac23 n^3$ for the full method. For a thousand unknowns that’s two million operations instead of almost seven hundred million. If swaps are needed, the decomposition is written as $PA = LU$: the matrix $P$ rearranges the rows of $A$ as pivoting would, after which forward elimination goes through without swaps. That’s the form in which Gaussian elimination lives in numerical libraries. LU decomposition was described in 1938 by the Polish astronomer Tadeusz Banachiewicz. In 1948 Alan Turing analysed how rounding errors accumulate in it and introduced the term “condition number”.

Ceres found

Back to 1801. Gauss worked out how to recover an orbit from three observations, computed where to look for Ceres, and sent the result to the astronomer Franz Xaver von Zach. On 31 December 1801 von Zach, and the next night Heinrich Olbers too, found it less than half a degree from the place Gauss had indicated. The twenty-four-year-old mathematician became a famous astronomer overnight. A few years later he took charge of the observatory in Göttingen.

He set out his method in 1809 in the book Theory of the Motion of the Heavenly Bodies Moving about the Sun in Conic Sections. The same book has the method of least squares: how to reconcile dozens of observations, each with its own error. It, and Gauss’s priority dispute with Legendre, will be the subject of chapter 39. The method of least squares leads to systems of linear equations, and Gauss explained how to solve them by successive elimination. In 1810, calculating the orbit of Pallas, the second asteroid, discovered by Olbers in the spring of 1802, he devised a compact notation for elimination. Professional computers used it throughout the nineteenth century.

The method itself is almost two thousand years older than Gauss: elimination of unknowns is described in the eighth chapter of the Chinese Nine Chapters on the Mathematical Art, and we saw the counting board with its columns of rods in the chapter on systems. In Europe it was taught from Isaac Newton’s Universal Arithmetick, published in 1707 from his lectures. By the end of the eighteenth century elimination of unknowns had become a routine algebra lesson. Gauss’s name became attached to the school algorithm only in the 1950s, in the era of the first computers, largely through confusion about its history.

Today Ceres is a dwarf planet, the largest body in the asteroid belt between Mars and Jupiter. From 2015 to 2018 the Dawn probe orbited it and photographed its bright spots, deposits of salts. The probe’s trajectory was refined from measurements by the method of least squares, a direct descendant of Gauss’s calculations.

Where next

Very different things were being added up in this chapter. We added the rows of a matrix to one another like vectors. Solutions of a homogeneous system can be added or multiplied by a number, and the result is again a solution. The parabolas $y = a + bt + ct^2$ add coefficient by coefficient: the sum of the parabolas from our two examples, $3 + 2t + t^2$ and $2 - t + t^2$, is again a polynomial of the same form, and it passes through points whose heights are the sums of the heights. Sound signals add when two instruments play, and are multiplied by a number when the volume knob is turned. Arrows, rows, polynomials, signals: why do they all behave the same way, and what can be understood about all such worlds at once? The answer is a single definition eight lines long. It’s where chapter 37 begins.