Lecture 05: The Algebra of Matrices

The idea in one sentence

There are four ways to read the product $AB$ , they always agree on the numbers, and the three you were not taught in school are the ones you will think with.

One product, four readings

Everything in this chapter runs on one pair:

A = \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{bmatrix}, \qquad B = \begin{bmatrix} 1 & 0 \\ 2 & 1 \\ 0 & 3 \end{bmatrix}, \qquad AB = \begin{bmatrix} 5 & 11 \\ 14 & 23 \end{bmatrix}.

Four times now you will get $(5, 11, 14, 23)$ back. What changes each time is which piece of $A$ meets which piece of $B$ , and therefore what the product is telling you.

Fig. 1

Four cameras on one product. Orange is what you take from $A$ and $B$ , green is what comes out, and that pairing holds for the rest of the chapter. Each reading hands back a different shape: one entry, one whole layer, one column, one row. Three of them are pieces of $[[5, 11], [14, 23]]$ you can read straight off; the second is a layer, and three layers add to the same answer.

Row times column. The school rule. Entry $(i,j)$ of $AB$ is row $i$ of $A$ dotted with column $j$ of $B$ :

(AB)_{11} = 1 \cdot 1 + 2 \cdot 2 + 3 \cdot 0 = 5, \qquad (AB)_{22} = 4 \cdot 0 + 5 \cdot 1 + 6 \cdot 3 = 23.

Four dot products, four numbers. This is how you compute by hand and how the innermost loop of a matrix library is written. It says nothing about what the product means.

Draw it once and the mechanics stop being a rule to remember. Put $B$ above the answer and $A$ beside it, and the row and the column you are pairing run straight into the cell they build.

Fig. 2

Where one entry comes from. The orange row of $A$ runs right, the orange column of $B$ runs down, and the green cell is where they meet. The three products are stacked on the right, one for each row of $B$ , so the index $k$ that vanishes in the sum is the one running down the page.

Column times row. Take column $k$ of $A$ and row $k$ of $B$ , multiply them the wide way, and add up over $k$ :

AB = \begin{bmatrix} 1 \\ 4 \end{bmatrix}\begin{bmatrix} 1 & 0 \end{bmatrix} + \begin{bmatrix} 2 \\ 5 \end{bmatrix}\begin{bmatrix} 2 & 1 \end{bmatrix} + \begin{bmatrix} 3 \\ 6 \end{bmatrix}\begin{bmatrix} 0 & 3 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 4 & 0 \end{bmatrix} + \begin{bmatrix} 4 & 2 \\ 10 & 5 \end{bmatrix} + \begin{bmatrix} 0 & 9 \\ 0 & 18 \end{bmatrix}.

Add the three and you get $AB$ again. Each piece is a whole matrix built from one column and one row, so both of its columns are multiples of $(1,4)$ , or of $(2,5)$ , or of $(3,6)$ . A matrix like that is called rank one: it has one direction in it, repeated. This reading says $AB$ is a sum of three rank-one layers. It is the least obvious of the four and by far the most valuable later; the whole of Part VI is about which layers of a matrix carry the information and which can be thrown away.

Fig. 3

The same product as a sum of three flat pieces. Each orange layer is one column of $A$ meeting one row of $B$ , so it carries a single direction and nothing more. Part VI lives here.

Matrix times column. Column $j$ of $AB$ is $A$ acting on column $j$ of $B$ , and Lecture 2 already told you what that is: a combination of the columns of $A$ , with the entries of that column as the weights. So column 1 of $AB$ is

1\begin{bmatrix} 1 \\ 4 \end{bmatrix} + 2\begin{bmatrix} 2 \\ 5 \end{bmatrix} + 0\begin{bmatrix} 3 \\ 6 \end{bmatrix} = \begin{bmatrix} 5 \\ 14 \end{bmatrix}.

That has a consequence you get for free: every column of $AB$ is a mix of columns of $A$ , so $AB$ can never reach anywhere $A$ could not reach on its own. Multiplying on the right can only narrow the set of reachable outputs.

Row times matrix. The mirror image. Row $i$ of $AB$ is row $i$ of $A$ used as weights on the rows of $B$ , so row 1 of $AB$ is

1\begin{bmatrix} 1 & 0 \end{bmatrix} + 2\begin{bmatrix} 2 & 1 \end{bmatrix} + 3\begin{bmatrix} 0 & 3 \end{bmatrix} = \begin{bmatrix} 5 & 11 \end{bmatrix}.

This is the reading Lecture 4 leaned on the whole way through. An elimination matrix $E$ is the identity with one entry changed, and row 2 of $EA$ is row 2 of $E$ used as weights on the rows of $A$ . Change one weight and you change one row of the answer. That is why a matrix on the left moves rows.

Move between the readings until they stop feeling like four rules and start feeling like one object seen from four sides:

Try it. One product, four readings

Orange is what you take in, green is what comes out. The numbers never change.

Order matters

Take two matrices that both do something simple to the rows of whatever they touch. $S$ adds row 2 into row 1. $T$ adds row 1 into row 2:

S = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix}, \qquad T = \begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix}, \qquad ST = \begin{bmatrix} 2 & 1 \\ 1 & 1 \end{bmatrix}, \qquad TS = \begin{bmatrix} 1 & 1 \\ 1 & 2 \end{bmatrix}.

Different matrices, from the same two factors. The reason is easy to say out loud. In $ST$ the matrix $T$ is closer to whatever comes next, so $T$ acts first: row 1 gets added into row 2, and then the already-changed row 2 gets added into row 1. In $TS$ the first step is the other one, so the second step works on a row nobody has touched yet. Operations that interfere cannot be reordered.

Fig. 4

The yellow square is where both products start, and the blue arrows show where its two sides land. The green regions have the same area and different shapes. Two shears in two orders are two different machines.

Sometimes the failure is louder than that. With the $A$ and $B$ from the top of the chapter, $AB$ is 2 by 2 and $BA$ is 3 by 3, so the question of whether they are equal never gets asked. The one pair that always commutes is a matrix and the identity: $AI = IA = A$ , because doing nothing before and doing nothing after are the same nothing.

Grouping does not matter

Order is fixed. Bracketing is free:

A(BC) = (AB)C.

The reason is that both sides compute the same double sum. Entry $(i,j)$ is $\sum_k \sum_l a_{ik} b_{kl} c_{lj}$ either way; the brackets only decide which sum you do first. This law looks like bookkeeping and pays like an investment.

It is what let Lecture 4 collapse the chain $E_{32}E_{31}E_{21}A$ into one matrix $E$ acting on $A$ alone. It is also the entire reason an inverse is unique, in a single line: if $BA = I$ and $AC = I$ , then $B = BI = B(AC) = (BA)C = IC = C$ . Two candidate inverses have no room to differ.

And it decides how fast your code runs.

Fig. 5

Both groupings give the same answer for $1000 \times 1000$ matrices. Pink is the pair to group first if you enjoy waiting; green is the one to group first otherwise. Never build a matrix you are only going to hit with a vector.

Cutting a matrix into blocks

Nothing in the four readings cared that the pieces were single rows and columns. Cut $A$ into rectangular blocks instead, cut $B$ to match, and the same rule holds with numbers replaced by blocks:

\begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix} \begin{bmatrix} B_{11} \\ B_{21} \end{bmatrix} = \begin{bmatrix} A_{11}B_{11} + A_{12}B_{21} \\ A_{21}B_{11} + A_{22}B_{21} \end{bmatrix}.

The one condition is that the vertical cuts in $A$ line up with the horizontal cuts in $B$ . If $A$ is split after column 2, then $B$ has to be split after row 2, so that $A_{11}$ has as many columns as $B_{11}$ has rows.

Fig. 6

Block multiplication is the ordinary rule with numbers replaced by blocks. Orange meets orange and green meets green, which is what the blue cuts have to guarantee.

You have already used this. The augmented matrix $[\,A \;|\; b\,]$ from Lecture 4 is a block matrix with one skinny block glued on the right, and $E[\,A \;|\; b\,] = [\,EA \;|\; Eb\,]$ is block multiplication written out. Deep learning uses the same idea at scale: a batch of inputs is stacked into one wide block so that one call to the hardware does the work of a thousand.

The shape the hardware actually wants

Blocks are not only a proof device. A library never multiplies two matrices entry by entry: it cuts the answer into tiles and builds one tile at a time, because a tile is small enough to sit in fast memory while the panels stream past it.

Fig. 7

How the machine multiplies. Orange is what you take in, green is what comes out, and the yellow edge is the shared dimension $k$ that meets and vanishes. One tile of $C$ needs a panel of rows from $A$ and a panel of columns from $B$ ; the loop then takes the next tile, and the next, until $C$ is full. This loop has a name, GEMM, and it is where the time goes.

Now put a real layer through it. A batch of 32 images, each flattened to 784 numbers, meets a weight matrix that turns 784 numbers into 128.

Fig. 8

One layer, one call. The 784 that both operands share is the yellow edge, and it is gone from the answer, exactly as $k$ was. Every forward pass in every network you will meet is this picture, repeated once per layer.

The four readings from the top of this chapter are four ways to describe what happens inside one tile. The machine picks whichever one its memory prefers, and the answer is the same either way.

Undoing a matrix

A matrix takes $x$ to $Ax$ . The natural question is whether anything takes $Ax$ back to $x$ . If some matrix $A^{-1}$ does, for every $x$ , then

A^{-1}A = AA^{-1} = I,

and $A$ is called invertible. For square matrices either equation implies the other, so you only ever have to check one.

When is there no such matrix? Suppose $Ax = 0$ for some $x$ that is not the zero vector. $A$ also sends $0$ to $0$ . So two different inputs produce the same output, and no machine on earth can send that one output back to both places. Lecture 2 met this already: the ring road machine sent $(1, 1, 1)$ to zero because raising all three cities by a meter changed no climb, and that is exactly why it could not be undone. Singular and “kills a nonzero vector” are the same condition wearing different clothes.

Fig. 9

What singular costs you. Blue is the map. On the left two different inputs land on the same pink point, so an arrow back would have to choose between them. On the right nothing was lost, and every blue arrow has a green one returning beside it.

At size 2 you can write the answer down. For

A = \begin{bmatrix} a & b \\ c & d \end{bmatrix}, \qquad A^{-1} = \frac{1}{ad - bc}\begin{bmatrix} d & -b \\ -c & a \end{bmatrix}.

Try it on $A = [[3, 1], [4, 2]]$ , where $ad - bc = 3 \cdot 2 - 1 \cdot 4 = 2$ :

A^{-1} = \frac{1}{2}\begin{bmatrix} 2 & -1 \\ -4 & 3 \end{bmatrix} = \begin{bmatrix} 1 & -0.5 \\ -2 & 1.5 \end{bmatrix}, \qquad \begin{bmatrix} 3 & 1 \\ 4 & 2 \end{bmatrix} \begin{bmatrix} 1 & -0.5 \\ -2 & 1.5 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}.

The number $ad - bc$ is the determinant. The formula divides by it, so a zero determinant is a matrix with no inverse. At size 2 that happens exactly when the two columns point along the same line and the machine flattens the plane. Lecture 14 builds the determinant properly, for every size, out of that one idea.

Inverses come in reverse order

Put on socks, then shoes. To undo that you take off the shoes first. Matrices are the same:

(AB)^{-1} = B^{-1}A^{-1}.

Fig. 10

Orange is doing, green is undoing. Walking a chain backwards means the last thing done is the first thing undone, so the bottom row is the top row reversed.

The proof is one line of associativity used twice:

(AB)(B^{-1}A^{-1}) = A(BB^{-1})A^{-1} = AIA^{-1} = AA^{-1} = I.

Check it on the two shears. $(ST)^{-1}$ and $T^{-1}S^{-1}$ both come out as $[[1, -1], [-1, 2]]$ , while $S^{-1}T^{-1}$ comes out as $[[2, -1], [-1, 1]]$ , a different matrix. The order in the answer really is the reverse of the order in the question.

Computing an inverse, once

Here is the honest way to find $A^{-1}$ by hand. Elimination does all the work; what changes is where you point it. $AA^{-1} = I$ says that $A$ times column $j$ of $A^{-1}$ equals column $j$ of $I$ . So the columns of the inverse are the solutions of $n$ systems that all share the same matrix. Solve them together: glue the whole identity on as extra columns and run elimination on the block $[\,A \;|\; I\,]$ . When the left half becomes $I$ , the right half has become $A^{-1}$ . Elimination downward makes the staircase, elimination upward clears above the pivots, and one division per row finishes it. That last upward sweep is the only new move, and it is why the method carries Jordan’s name alongside Gauss’s.

Run it on the matrix this course will keep coming back to, the second-difference matrix:

K = \begin{bmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{bmatrix}.

Downward, the pivots are $2$ , then $\tfrac{3}{2}$ , then $\tfrac{4}{3}$ . Upward and scaled, the right half settles at

K^{-1} = \frac{1}{4}\begin{bmatrix} 3 & 2 & 1 \\ 2 & 4 & 2 \\ 1 & 2 & 3 \end{bmatrix}.

That answer has more in it than the nine numbers.

The pivots multiply to $2 \cdot \tfrac{3}{2} \cdot \tfrac{4}{3} = 4$ , and 4 is the determinant of $K$ . Every entry of the inverse is then divided by that 4. A matrix with determinant zero has nothing to divide by, which is the same sentence as “no inverse”, arrived at from a third direction.

$K$ is symmetric and so is $K^{-1}$ . That always happens, and it will matter a great deal from Lecture 15 onward.

And $K$ is as sparse as a matrix can be while still being connected: seven nonzeros out of nine, and every nonzero sits on or beside the diagonal. Its inverse has nine out of nine. Grow the matrix and the contrast grows with it.

Fig. 11

Sparse in, dense out. Every entry of the green inverse depends on every entry of the orange original, so the band spreads to fill the whole square. Depth of colour is size.

Why nobody calls inv

Start with the picture above. A sparse matrix of size a million has a few million nonzeros and fits in memory. Its inverse has a trillion entries and does not exist on any machine you own.

Then there is speed. The inverse costs a full elimination plus $n$ back substitutions, and then you still have to multiply by $b$ . Solving directly costs one elimination and one back substitution. On a $3000 \times 3000$ system, measured on the machine that wrote this chapter, inverting and multiplying took 0.20 seconds while solving took 0.067.

The reason that bites silently is accuracy. The code below solves a $12 \times 12$ system whose true answer is all ones, using a matrix chosen to be nasty. Going through the inverse gives an answer wrong by 12.6 in its worst coordinate. Solving directly is wrong by 0.22. Both are bad, because the matrix is genuinely awful, and one of them is nearly sixty times worse for no benefit.

$A^{-1}$ is how you think. solve is how you compute. The code below shows the four readings agreeing, then the inverse story end to end.

python

import numpy as np

A = np.array([[1., 2., 3.], [4., 5., 6.]])
B = np.array([[1., 0.], [2., 1.], [0., 3.]])

print(A @ B)                                                  # [[ 5. 11.] [14. 23.]]
print(np.array([[A[i] @ B[:, j] for j in range(2)]            # rows dot columns
                for i in range(2)]))                          # [[ 5. 11.] [14. 23.]]
print(sum(np.outer(A[:, k], B[k]) for k in range(3)))         # rank-one layers
                                                              # [[ 5. 11.] [14. 23.]]
print(np.column_stack([A @ B[:, j] for j in range(2)]))       # columns of A, mixed
                                                              # [[ 5. 11.] [14. 23.]]
print(np.vstack([A[i] @ B for i in range(2)]))                # rows of B, mixed
                                                              # [[ 5. 11.] [14. 23.]]

S = np.array([[1., 1.], [0., 1.]])
T = np.array([[1., 0.], [1., 1.]])
print(S @ T)                    # [[2. 1.] [1. 1.]]
print(T @ S)                    # [[1. 1.] [1. 2.]]   order is not negotiable

python

import torch

A = torch.tensor([[1., 2., 3.], [4., 5., 6.]])
B = torch.tensor([[1., 0.], [2., 1.], [0., 3.]])

print(A @ B)                                                  # [[ 5. 11.] [14. 23.]]
print(torch.stack([torch.stack([A[i] @ B[:, j] for j in range(2)])
                   for i in range(2)]))                       # rows dot columns
                                                              # [[ 5. 11.] [14. 23.]]
print(sum(torch.outer(A[:, k], B[k]) for k in range(3)))      # rank-one layers
                                                              # [[ 5. 11.] [14. 23.]]
print(torch.stack([A @ B[:, j] for j in range(2)], dim=1))    # columns of A, mixed
                                                              # [[ 5. 11.] [14. 23.]]
print(torch.stack([A[i] @ B for i in range(2)]))              # rows of B, mixed
                                                              # [[ 5. 11.] [14. 23.]]

S = torch.tensor([[1., 1.], [0., 1.]])
T = torch.tensor([[1., 0.], [1., 1.]])
print(S @ T)                    # [[2. 1.] [1. 1.]]
print(T @ S)                    # [[1. 1.] [1. 2.]]   order is not negotiable

python

import jax
jax.config.update("jax_enable_x64", True)   # jax is single precision unless you ask
import jax.numpy as jnp

A = jnp.array([[1., 2., 3.], [4., 5., 6.]])
B = jnp.array([[1., 0.], [2., 1.], [0., 3.]])

print(A @ B)                                                  # [[ 5. 11.] [14. 23.]]
print(jnp.einsum('ik,kj->ij', A, B))                          # rows dot columns
                                                              # [[ 5. 11.] [14. 23.]]
print(sum(jnp.outer(A[:, k], B[k]) for k in range(3)))        # rank-one layers
                                                              # [[ 5. 11.] [14. 23.]]
print(jnp.stack([A @ B[:, j] for j in range(2)], axis=1))     # columns of A, mixed
                                                              # [[ 5. 11.] [14. 23.]]
print(jnp.stack([A[i] @ B for i in range(2)]))                # rows of B, mixed
                                                              # [[ 5. 11.] [14. 23.]]

S = jnp.array([[1., 1.], [0., 1.]])
T = jnp.array([[1., 0.], [1., 1.]])
print(S @ T)                    # [[2. 1.] [1. 1.]]
print(T @ S)                    # [[1. 1.] [1. 2.]]   order is not negotiable

python

import tensorflow as tf

A = tf.constant([[1., 2., 3.], [4., 5., 6.]])
B = tf.constant([[1., 0.], [2., 1.], [0., 3.]])

print(A @ B)                                                  # [[ 5. 11.] [14. 23.]]
print(tf.einsum('ik,kj->ij', A, B))                           # rows dot columns
                                                              # [[ 5. 11.] [14. 23.]]
print(sum(tf.einsum('i,j->ij', A[:, k], B[k]) for k in range(3)))
                                                              # rank-one layers
                                                              # [[ 5. 11.] [14. 23.]]
print(tf.stack([A @ B[:, k:k+1] for k in range(2)], axis=1)[..., 0])
                                                              # columns of A, mixed
                                                              # [[ 5. 11.] [14. 23.]]
print(tf.concat([A[i:i+1] @ B for i in range(2)], axis=0))    # rows of B, mixed
                                                              # [[ 5. 11.] [14. 23.]]

S = tf.constant([[1., 1.], [0., 1.]])
T = tf.constant([[1., 0.], [1., 1.]])
print(S @ T)                    # [[2. 1.] [1. 1.]]
print(T @ S)                    # [[1. 1.] [1. 2.]]   order is not negotiable

The second block is Gauss-Jordan written out, then the measurements behind the paragraphs above. Only NumPy is timed, because a stopwatch reading is a fact about one machine.

python

import numpy as np

K = np.array([[2., -1., 0.], [-1., 2., -1.], [0., -1., 2.]])
M = np.hstack([K, np.eye(3)])          # the block [K | I], three rows of six

for k in range(3):                     # down: clear below every pivot
    for i in range(k + 1, 3):
        M[i] -= (M[i, k] / M[k, k]) * M[k]
pivots = np.diag(M[:, :3]).copy()
print(pivots)                          # [2.  1.5  1.33333333]
print(np.prod(pivots), np.linalg.det(K))   # 4.0 4.0   pivots multiply to det

for k in (2, 1):                       # up: clear above every pivot
    for i in range(k):
        M[i] -= (M[i, k] / M[k, k]) * M[k]
M /= np.diag(M[:, :3])[:, None]        # one division per row

print(M[:, :3])                        # the identity, so the left half is done
print(4 * M[:, 3:])                    # [[3. 2. 1.] [2. 4. 2.] [1. 2. 3.]]
print(np.allclose(M[:, 3:], np.linalg.inv(K)))   # True

K8 = (np.diag(2 * np.ones(8)) - np.diag(np.ones(7), 1) - np.diag(np.ones(7), -1))
print(np.count_nonzero(K8))                             # 22   sparse going in
print(np.sum(np.abs(np.linalg.inv(K8)) > 1e-9))          # 64   dense coming out

python

import torch

K = torch.tensor([[2., -1., 0.], [-1., 2., -1.], [0., -1., 2.]])
M = torch.cat([K, torch.eye(3)], dim=1)   # the block [K | I], three rows of six

for k in range(3):                     # down: clear below every pivot
    for i in range(k + 1, 3):
        M[i] -= (M[i, k] / M[k, k]) * M[k]
pivots = torch.diag(M[:, :3]).clone()
print(pivots)                          # tensor([2.0000, 1.5000, 1.3333])
print(pivots.prod(), torch.linalg.det(K))  # tensor(4.) tensor(4.)

for k in (2, 1):                       # up: clear above every pivot
    for i in range(k):
        M[i] -= (M[i, k] / M[k, k]) * M[k]
M /= torch.diag(M[:, :3])[:, None]     # one division per row

print(M[:, :3])                        # the identity, so the left half is done
print(4 * M[:, 3:])                    # [[3. 2. 1.] [2. 4. 2.] [1. 2. 3.]]
print(torch.allclose(M[:, 3:], torch.linalg.inv(K), atol=1e-6))   # True

K8 = (torch.diag(2 * torch.ones(8)) - torch.diag(torch.ones(7), 1)
      - torch.diag(torch.ones(7), -1))
print(torch.count_nonzero(K8))                       # tensor(22)  sparse going in
print((torch.linalg.inv(K8).abs() > 1e-9).sum())     # tensor(64)  dense coming out

python

import jax
jax.config.update("jax_enable_x64", True)   # jax is single precision unless you ask
import jax.numpy as jnp

K = jnp.array([[2., -1., 0.], [-1., 2., -1.], [0., -1., 2.]])
M = jnp.hstack([K, jnp.eye(3)])        # the block [K | I], three rows of six

for k in range(3):                     # down: clear below every pivot
    for i in range(k + 1, 3):
        M = M.at[i].add(-(M[i, k] / M[k, k]) * M[k])
pivots = jnp.diag(M[:, :3])
print(pivots)                          # [2.  1.5  1.33333333]
print(jnp.prod(pivots), jnp.linalg.det(K))   # 4.0 4.0

for k in (2, 1):                       # up: clear above every pivot
    for i in range(k):
        M = M.at[i].add(-(M[i, k] / M[k, k]) * M[k])
M = M / jnp.diag(M[:, :3])[:, None]    # one division per row

print(M[:, :3])                        # the identity, so the left half is done
print(4 * M[:, 3:])                    # [[3. 2. 1.] [2. 4. 2.] [1. 2. 3.]]
print(jnp.allclose(M[:, 3:], jnp.linalg.inv(K), atol=1e-6))   # True

K8 = (jnp.diag(2 * jnp.ones(8)) - jnp.diag(jnp.ones(7), 1) - jnp.diag(jnp.ones(7), -1))
print(jnp.count_nonzero(K8))                            # 22   sparse going in
print(jnp.sum(jnp.abs(jnp.linalg.inv(K8)) > 1e-9))      # 64   dense coming out

python

import numpy as np
import tensorflow as tf

K = tf.constant([[2., -1., 0.], [-1., 2., -1.], [0., -1., 2.]])
M = np.hstack([K.numpy(), np.eye(3)])  # loop in numpy: tf tensors do not mutate

for k in range(3):                     # down: clear below every pivot
    for i in range(k + 1, 3):
        M[i] -= (M[i, k] / M[k, k]) * M[k]
pivots = tf.constant(np.diag(M[:, :3]).copy())
print(pivots)                          # [2.  1.5  1.33333333]
print(tf.reduce_prod(pivots), tf.linalg.det(K))   # 4.0 4.0

for k in (2, 1):                       # up: clear above every pivot
    for i in range(k):
        M[i] -= (M[i, k] / M[k, k]) * M[k]
M /= np.diag(M[:, :3])[:, None]        # one division per row

print(M[:, :3])                        # the identity, so the left half is done
print(4 * M[:, 3:])                    # [[3. 2. 1.] [2. 4. 2.] [1. 2. 3.]]
print(np.allclose(M[:, 3:], tf.linalg.inv(K).numpy()))   # True

K8 = (np.diag(2 * np.ones(8)) - np.diag(np.ones(7), 1) - np.diag(np.ones(7), -1))
print(tf.math.count_nonzero(K8))                     # 22   sparse going in
print(tf.math.count_nonzero(tf.abs(tf.linalg.inv(K8)) > 1e-9))
                                                     # 64   dense coming out

python

import numpy as np
from time import perf_counter

def best(f, reps=3):
    out = float("inf")
    for _ in range(reps):
        t0 = perf_counter(); f(); out = min(out, perf_counter() - t0)
    return out

n, rng = 3000, np.random.default_rng(0)
A = rng.standard_normal((n, n)) + np.sqrt(n) * np.eye(n)
b = rng.standard_normal(n)

print(best(lambda: np.linalg.inv(A) @ b))   # 0.20  seconds, on the machine that wrote this
print(best(lambda: np.linalg.solve(A, b)))  # 0.067 seconds

n = 12                                      # a matrix chosen to be numerically awful
H = 1.0 / (np.arange(1, n + 1)[:, None] + np.arange(1, n + 1) - 1)
x_true = np.ones(n)
c = H @ x_true
print(np.abs(np.linalg.inv(H) @ c - x_true).max())    # 12.59375
print(np.abs(np.linalg.solve(H, c) - x_true).max())   # 0.22008972067760946

Where this is going

Elimination gave you a way to solve one system. This chapter gave you the algebra to talk about what elimination did: a product of simple matrices sitting to the left of $A$ , regrouped into one matrix, and inverted by walking the chain backwards. Put those together and a question appears. If the expensive part of solving $Ax = b$ is making the staircase, and the staircase is a product, why not compute the product once, write it down, and reuse it for every $b$ that ever arrives?

That is a factorization, and the next chapter delivers two of them. $A = LU$ is elimination’s receipt: $U$ is the staircase from Lecture 4 and $L$ holds the multipliers, unmixed, for the reason this chapter half revealed. $A = CR$ goes the other way, building $A$ out of its independent columns, and it answers a question you have not been able to ask yet: how much genuinely different information is in a matrix at all?

Check yourself

With $A$ and $B$ from the top of the chapter, compute row 2 of $AB$ twice: once as two dot products, once as a combination of the rows of $B$ . The second way should feel faster.
The three rank-one layers of $AB$ are $[[1,0],[4,0]]$ , $[[4,2],[10,5]]$ and $[[0,9],[0,18]]$ . Which single layer would you keep if you had to throw two away, and what did you use to decide?
Find a pair of $2 \times 2$ matrices, neither one the identity, that do commute. What do yours have in common?
Invert $[[2, 1], [1, 1]]$ with the formula, then check your answer by multiplying. What is the determinant, and what does it say about the columns?
$(ABC)^{-1}$ is which product of $A^{-1}$ , $B^{-1}$ and $C^{-1}$ ? Prove it by applying the two-matrix rule twice.
Redraw Figure 1 from memory, all four panels and the green piece each one builds. The drawing is the test.