Solving Systems: Elimination and Inverses

Discover: a puzzle from two thousand years ago

About two thousand years ago, Chinese scholars compiled The Nine Chapters on the Mathematical Art, a handbook of practical problems. Its eighth chapter opens with one about grain. In modern words:

Three bundles of top-grade grain, two of medium grade and one of low grade yield 39 dou of grain. Two top, three medium and one low yield 34 dou. One top, two medium and three low yield 26 dou. How much does one bundle of each grade yield?

Call the yields of one bundle $x_1$ (top), $x_2$ (medium) and $x_3$ (low). The problem becomes three equations in three unknowns:

The scribes laid the numbers out on a counting board, one equation per column of rods, and called the grid a fangcheng, a “rectangular array”. Then they combined columns to wipe out one unknown after another. It is, step for step, the method of this chamber, written down more than fifteen centuries before Gauss was born.

You can try the key move right now.

One subtraction and $x_3$ has vanished. Do that systematically and you can solve any system. This chamber turns the trick into an algorithm, then asks the deeper questions behind it: when does a system have exactly one answer, when can a matrix be undone, and how much does it stretch space? That last question turns out to be the engine of a family of generative models, the normalising flows, which you'll meet at the end.

Learn: two pictures of Ax = b

Stack the coefficients into a matrix and the grain problem becomes one equation, $A\mathbf{x} = \mathbf{b}$, read “A x equals b”:

There are two ways to see what this equation asks, and you need both.

The row picture. Each equation on its own is a flat shape. In two unknowns, $2x_1 - x_2 = 1$ is a line in the plane; in three unknowns, each equation is a plane in space. A solution must lie on all of them at once, so solving means finding where they meet. The lines $2x_1 - x_2 = 1$ and $x_1 + x_2 = 5$ cross at the single point $(2, 3)$.

The column picture. In Chamber 4 you read $A\mathbf{x}$ as a combination of the columns of $A$, weighted by the entries of $\mathbf{x}$. So the same system asks: which combination of the columns makes $\mathbf{b}$?

The weights $2$ and $3$ are the solution again, seen from the other side. The column picture also tells you when a solution exists at all: exactly when $\mathbf{b}$ lies in the span of the columns (Chamber 3). Chamber 6 gives that span a name.

In the row picture you can see every possible outcome with two lines. They cross once (one solution), they are parallel (no solution), or they are the same line (infinitely many). Three planes in space behave the same way: there are never exactly two solutions, or exactly five.

Learn: Gaussian elimination

Elimination works on the augmented matrix $[A \mid \mathbf{b}]$ (“A augmented with b”): the coefficients with the right-hand side attached as an extra column. It never touches the unknowns directly. It only changes rows, and it uses just three row operations:

Why is this allowed? Because every row operation can be undone by another one: a swap by the same swap, $R_i \leftarrow c\,R_i$ by $R_i \leftarrow \tfrac{1}{c}R_i$ (which is why $c = 0$ is banned), and $R_i \leftarrow R_i + c\,R_j$ by $R_i \leftarrow R_i - c\,R_j$. Any solution of the old system also solves the new one, and since the step can be reversed, any solution of the new one solves the old. The set of solutions never changes; only the equations get simpler.

The plan: go column by column, left to right. In each column choose a non-zero pivot, then subtract multiples of the pivot's row to create zeros below it. Here is the grain problem.

1. 1

   Put a convenient pivot on top

   The first pivot could be the $3$, but a $1$ keeps the arithmetic clean. Swap rows 1 and 3:

   $$\left[\begin{array}{ccc|c} 3 & 2 & 1 & 39 \\ 2 & 3 & 1 & 34 \\ 1 & 2 & 3 & 26 \end{array}\right] \xrightarrow{R_1 \leftrightarrow R_3} \left[\begin{array}{ccc|c} 1 & 2 & 3 & 26 \\ 2 & 3 & 1 & 34 \\ 3 & 2 & 1 & 39 \end{array}\right]$$
2. 2

   Clear the first column below the pivot

   The multipliers are the entry to be cleared divided by the pivot: $2/1$ and $3/1$.

   $$\xrightarrow[R_3 \leftarrow R_3 - 3R_1]{R_2 \leftarrow R_2 - 2R_1} \left[\begin{array}{ccc|c} 1 & 2 & 3 & 26 \\ 0 & -1 & -5 & -18 \\ 0 & -4 & -8 & -39 \end{array}\right]$$
3. 3

   Move to the second column

   The second pivot is the $-1$. The multiplier for row 3 is $(-4)/(-1) = 4$:

   $$\xrightarrow{R_3 \leftarrow R_3 - 4R_2} \left[\begin{array}{ccc|c} 1 & 2 & 3 & 26 \\ 0 & -1 & -5 & -18 \\ 0 & 0 & 12 & 33 \end{array}\right]$$

   The matrix is now upper triangular: zeros below the diagonal. The pivots are $1$, $-1$ and $12$.

4. 4

   Back substitution, from the bottom up

   The last row says $12x_3 = 33$, so $x_3 = \tfrac{11}{4}$. The middle row says $-x_2 - 5x_3 = -18$, so $x_2 = 18 - \tfrac{55}{4} = \tfrac{17}{4}$. The top row says $x_1 + 2x_2 + 3x_3 = 26$, so $x_1 = 26 - \tfrac{34}{4} - \tfrac{33}{4} = \tfrac{37}{4}$.

One bundle yields $9\tfrac14$, $4\tfrac14$ and $2\tfrac34$ dou: the answer printed in the Nine Chapters. Check it in the first equation: $3 \cdot \tfrac{37}{4} + 2 \cdot \tfrac{17}{4} + \tfrac{11}{4} = \tfrac{156}{4} = 39$.

A pivot must be non-zero, since you divide by it. When the pivot position holds a zero, swap in a row from below that has a non-zero entry there. When there is no such row, that column simply gets no pivot, and as you'll see next, that is where infinitely many solutions (or none) come from.

Learn: echelon forms and the three outcomes

Elimination aims for a staircase shape called row echelon form. A matrix is in row echelon form when

1. any all-zero rows sit at the bottom, and
2. each row's first non-zero entry, its pivot, is strictly to the right of the pivot in the row above.

It follows that every entry below a pivot is zero. Stopping here and back-substituting is Gaussian elimination. Keep going and you reach the reduced row echelon form, $\operatorname{rref}$: in addition, every pivot equals $1$ and is the only non-zero entry in its column. Clearing upwards as well as downwards like this is called Gauss–Jordan elimination. Every matrix has exactly one reduced row echelon form, however you get there, so it's a fingerprint of the system.

Once $[A \mid \mathbf{b}]$ is in reduced form, you can read the outcome straight off it:

• Left: a pivot in every column of $A$. Exactly one solution, and the last column is it: $\mathbf{x} = (2, -1, 4)$.
• Middle: a column without a pivot. Column 2 has no pivot, so $x_2$ is a free variable: choose it freely, say $x_2 = t$, and the pivot rows fix the rest, $x_1 = 2 - 3t$ and $x_3 = 5$. Infinitely many solutions, $\mathbf{x} = (2, 0, 5) + t\,(-3, 1, 0)$: a line. The zero row, $0 = 0$, is harmless.
• Right: a pivot in the last column. The bottom row says $0x_1 + 0x_2 + 0x_3 = 1$, that is, $0 = 1$. No solution: the system is inconsistent.

That's the whole story, for any number of equations and unknowns. A pivot in the $\mathbf{b}$ column means none; otherwise, every free column means infinitely many; no free columns means exactly one.

Now do it yourself. The workbench records each operation in the notation above, and it recognises both echelon forms. One of the three systems has a single solution, one has infinitely many and one has none. Reduce all three to find out which is which. If you're stuck, the Hint button loads the next Gauss–Jordan step.

Learn: undoing a matrix

In Chamber 4, a matrix was a transformation of space. Some transformations can be reversed: rotate by $30°$ and you can rotate back. Others can't: a projection squashes the plane onto a line, and once two different points land in the same place, no transformation can tell them apart again.

A square matrix $A$ is invertible if there is a matrix $A^{-1}$, read “A inverse”, with

Then every system $A\mathbf{x} = \mathbf{b}$ has exactly one solution, $\mathbf{x} = A^{-1}\mathbf{b}$. (You'll put the proof of that together in the last part of this chamber.) A square matrix without an inverse is called singular.

Computing an inverse. The $j$-th column of $A^{-1}$ is the solution of $A\mathbf{x} = \mathbf{e}_j$, where $\mathbf{e}_j$ is the $j$-th standard basis vector. You could solve $n$ systems one at a time, but they share the same $A$, so do them all at once: reduce the wide augmented matrix $[A \mid I]$. If $A$ is invertible, the left half becomes $I$ and the right half becomes $A^{-1}$. If elimination runs out of pivots first, $A$ is singular.

Let's run that on a general $2 \times 2$ matrix and watch the famous formula appear. (Here $b$ is just the top-right entry of the matrix, not a right-hand side.)

1. 1

   Set up, assuming a ≠ 0 for now

   $$\left[\begin{array}{cc|cc} a & b & 1 & 0 \\ c & d & 0 & 1 \end{array}\right]$$
2. 2

   Clear below the first pivot

   Use $R_2 \leftarrow R_2 - \tfrac{c}{a}R_1$. The new $(2,2)$ entry is $d - \tfrac{bc}{a} = \tfrac{ad - bc}{a}$:

   $$\left[\begin{array}{cc|cc} a & b & 1 & 0 \\ 0 & \frac{ad - bc}{a} & -\frac{c}{a} & 1 \end{array}\right]$$

   The second pivot is $\tfrac{ad - bc}{a}$. If $ad - bc = 0$ there is no second pivot, and $A$ has no inverse. From here on, assume $ad - bc \neq 0$.

3. 3

   Make the second pivot 1

   Scale by its reciprocal, $R_2 \leftarrow \tfrac{a}{ad - bc}R_2$:

   $$\left[\begin{array}{cc|cc} a & b & 1 & 0 \\ 0 & 1 & \frac{-c}{ad - bc} & \frac{a}{ad - bc} \end{array}\right]$$
4. 4

   Clear above the second pivot

   $R_1 \leftarrow R_1 - bR_2$. The new top-left entry on the right is $1 + \tfrac{bc}{ad - bc} = \tfrac{ad}{ad - bc}$:

   $$\left[\begin{array}{cc|cc} a & 0 & \frac{ad}{ad - bc} & \frac{-ab}{ad - bc} \\ 0 & 1 & \frac{-c}{ad - bc} & \frac{a}{ad - bc} \end{array}\right]$$
5. 5

   Make the first pivot 1 and read off the inverse

   $R_1 \leftarrow \tfrac{1}{a}R_1$ turns the left half into $I$, and the right half is the inverse:

   $$\begin{bmatrix} a & b \\ c & d \end{bmatrix}^{-1} = \frac{1}{ad - bc}\begin{bmatrix} d & -b \\ -c & a \end{bmatrix}$$

Swap the diagonal, negate the off-diagonal, divide by $ad - bc$. (If $a = 0$, swap the rows first; the formula comes out the same.) Always check by multiplying:

That number $ad - bc$ decides everything: non-zero, the inverse exists; zero, it doesn't. It has a name, and a meaning, which is the subject of the next section.

Inverting a product. Here is a complete proof of a rule you'll use constantly.

Claim. If $A$ and $B$ are invertible $n \times n$ matrices, then $AB$ is invertible and $(AB)^{-1} = B^{-1}A^{-1}$.

Proof. To show that a matrix $C$ is the inverse of $AB$, it is enough to check that $(AB)C = I$ and $C(AB) = I$: that's the definition. Take $C = B^{-1}A^{-1}$. Using associativity (Chamber 4) to regroup,

Both products are $I$, so $B^{-1}A^{-1}$ is the inverse of $AB$. $\blacksquare$

The order reverses for the reason your shoes come off before your socks: $AB$ means “first $B$, then $A$”, so undoing it means “first undo $A$, then undo $B$”. Notice the proof's shape. It didn't derive the answer. It guessed a candidate and verified the definition, a pattern you'll see again and again.

Learn: the determinant as a volume scale

A $2 \times 2$ matrix sends the unit square to the parallelogram spanned by its columns (Chamber 4's picture: the columns are where $\mathbf{e}_1$ and $\mathbf{e}_2$ land). The determinant, $\det A$ (read “det A”), is the area of that parallelogram, with a sign:

Since the transformation is linear, every region is scaled by the same factor, so $|\det A|$ is the factor by which $A$ scales all areas. For example, $\begin{bmatrix} 3 & 1 \\ 1 & 2 \end{bmatrix}$ sends the unit square to a parallelogram of area $3 \cdot 2 - 1 \cdot 1 = 5$, and a circle of area $\pi$ to an ellipse of area $5\pi$. A negative determinant means the transformation also flips orientation, like a mirror. A zero determinant means the parallelogram is flat: the plane has been squashed onto a line, which is exactly the singular case you met in the $2 \times 2$ inverse.

For a $3 \times 3$ matrix the determinant is the signed volume of the box (a parallelepiped) spanned by its three columns. One formula expands along the first row, each entry times the $2 \times 2$ determinant left when you delete its row and column, with alternating signs:

For the grain matrix: $3(9 - 2) - 2(6 - 1) + 1(4 - 3) = 21 - 10 + 1 = 12$.

Expanding like this takes roughly $n!$ multiplications for an $n \times n$ matrix, which is hopeless beyond tiny sizes. Four rules give a far better way:

• Products multiply: $\det(AB) = \det A \cdot \det B$. Doing one transformation after another multiplies their scale factors.
• Triangular matrices are easy: the determinant is the product of the diagonal entries.
• Row operations have simple effects: a swap flips the sign, scaling a row by $c$ scales the determinant by $c$, and adding a multiple of one row to another changes nothing.
• The identity has determinant $1$.

Put together, elimination computes determinants. It turns $A$ into a triangular matrix using only swaps and replacements, so

For the grain problem: one swap and pivots $1$, $-1$, $12$ give $\det A = -(1 \cdot (-1) \cdot 12) = 12$, just as the expansion said. It also explains the central fact of this chamber. $\det A = 0$ exactly when a pivot is missing, which is exactly when elimination can't reach $I$, which is exactly when there is no inverse:

Why machine learning cares. A normalising flow is a generative model that turns simple noise into data, such as images, through an invertible function. Training it needs the exact probability of each training image, and when an invertible function $f$ maps data $x$ to noise $z = f(x)$, probabilities get rescaled by exactly how much $f$ stretches volume near $x$:

The matrix inside the determinant is the Jacobian of $f$, its matrix of partial derivatives (Mathematics for Machine Learning, Chamber 7). Near $x$, $f$ behaves like that matrix, so its determinant is $f$'s local volume scale. The problem is cost. A small colour image has $32 \times 32 \times 3 = 3072$ numbers, so the Jacobian is $3072 \times 3072$, and a general determinant of that size costs around twenty billion arithmetic operations, even by elimination. For every image. At every training step.

Real NVP's answer is to design each layer so that its Jacobian is triangular. A coupling layer copies the first $d$ inputs unchanged and scales and shifts the rest, using amounts computed from the first $d$. The paper writes out the Jacobian:

Zeros above the diagonal, the identity and a diagonal matrix on it. The complicated bottom-left block, which depends on neural networks, doesn't matter at all: the determinant is the product of the diagonal, $\exp\big(\sum_j s(x_{1:d})_j\big)$, and its logarithm is a plain sum. And because $\det(AB) = \det A \det B$ (the paper's Equation 11), a whole stack of layers just adds up their log-determinants. A few facts from this section turn an impossible computation into a cheap one.

Read beyond the course

Each of these covers this chamber's ideas from a different angle. Pick at least one, and mark it read when you're done.

Papers and lectures

Real NVP is a good first paper about generative models because its central idea is pure linear algebra. You can read the core today. Start with Section 3.1, the change of variable formula (Equations 2 and 3, which you've now seen). Then read Section 3.2 on coupling layers, which defines the layer in Equations 4 and 5, and Section 3.3, where Equation 6 is the triangular Jacobian and Equations 7–8 invert the layer. Finish with the short Section 3.5, whose Equations 11 and 12 are $\det(AB) = \det A \det B$ and the socks-and-shoes rule. Skip the masking scheme, the multi-scale architecture, batch normalisation and the experiments for now, but do look at the generated faces and bedrooms in the figures: every one came from inverting a stack of these layers.

Watch

The best visual explanation of the determinant there is. Watch for the moment a region's area is scaled by the same factor everywhere, and for why a flipped grid gives a negative sign.

Gilbert Strang's second lecture from MIT's 18.06. He runs elimination on a $3 \times 3$ system, then shows that each step is itself a matrix multiplication, an idea that leads to the LU factorisation used inside np.linalg.solve.

Your turn

First, translate between notation, code and the rules of determinants.

Proofs

The first proof comes in pieces: put the lines in order and leave out the ones that don't belong. The second you write yourself, and it's short.

Code it up

Three problems: count the singular weight matrices, beat floating point with exact elimination, and evaluate a Real NVP likelihood the way the paper does.