Determinants

The idea in one sentence

The determinant of a square matrix is the factor by which that matrix scales volume, and it is zero exactly when the matrix flattens space.

The number you have already dragged

In Lecture 2 you watched a machine with columns $(2, 0)$ and $(1, 1)$ pick up the unit square and put it down as a parallelogram. Underneath the picture a small readout said $\det A = 2$. Here is what it was telling you.

The unit square has area 1. Its image has area 2. That is the whole content of the number: feed the machine any region, and whatever area went in comes out multiplied by $\det A$. In three dimensions the same number multiplies volume. Nothing about the shape of the region matters, which is the surprising part. A square or a circle or the outline of a country all scale by the same factor, because a linear map treats every small patch of the plane the same way.

So the determinant carries one number and one bit. The number $|\det A|$ is the volume factor. The sign says whether the machine kept the handedness of space or turned it over, the way a mirror does. And when the columns line up, the parallelogram has no area at all, the machine has crushed the plane onto a line, and the determinant is zero. Zero determinant means singular, which is the test this whole course has been circling.

Drag the two columns together and watch the number die:

Three rules

You could define the determinant by a formula. It is better to define it by what it does, because then the formulas can be derived instead of memorized. Three rules are enough.

Rule 1. $\det I = 1$. The identity does nothing, so the unit box keeps its volume.

Rule 2. Exchanging two rows reverses the sign. Swapping two edges of the box turns it inside out and the volume factor changes sign.

Rule 3. The determinant is linear in each row separately, with the other rows held still. Two halves:

The first half is stretching one edge of the box while the others stay put, which multiplies the volume by $t$. The second half is hard to see in words and immediate in a picture.

Now use the rules. Rule 3 says scaling one row scales the determinant, so scaling every row of a $2 \times 2$ matrix scales it twice:

Pull the 2 out of row one, then out of row two, then apply Rule 1. In $n$ dimensions the same argument gives $\det(tI) = t^n$, which is the statement that doubling every edge of a box in three dimensions multiplies its volume by 8. This is the one place where $\det(2I) = 2\det I$ looks right and is not. Linear in each row separately does not mean linear in the matrix.

What falls out

Everything else is a consequence. Each one is worth checking on a $2 \times 2$ before you trust it in general.

Two equal rows give zero. Swap the two equal rows. The matrix is unchanged, so its determinant is unchanged. Rule 2 says the determinant flipped sign. The only number equal to its own negative is 0:

Geometrically the two edges point the same way, so the box is flat.

A row of zeros gives zero. Scale that row by $t = 0$ and use Rule 3.

Elimination changes nothing. Subtracting a multiple of one row from another leaves the determinant alone:

Rule 3 splits the changed row into the original plus a multiple of the row you subtracted, and the second piece is a matrix with two proportional rows, which is zero. In pictures it is a shear.

A triangular matrix is the product of its diagonal. Pull each diagonal entry out of its row with Rule 3, clear the off-diagonal entries with elimination, and land on $\det I = 1$:

Singular means zero. If $A$ is singular, elimination produces a row of zeros, so $\det A = 0$. If $A$ is invertible, elimination produces nonzero pivots all the way down, so $\det A \neq 0$. One number decides invertibility, which is what makes the determinant worth having at all.

Determinants multiply. $\det(AB) = (\det A)(\det B)$. Volume factors compose: do one machine then the other and the scalings multiply. Take $A$ with rows $(3, 1)$ and $(1, 2)$, and $B$ with rows $(0, 2)$ and $(-1, 4)$. Their determinants are 5 and 2, the product $AB$ has rows $(-1, 10)$ and $(-2, 10)$, and its determinant is 10. One corollary drops out at once: $AA^{-1} = I$ forces $\det(A^{-1}) = 1/\det A$.

Transposing changes nothing. $\det(A\T) = \det A$. That is why every rule stated for rows is equally true for columns, and why the column pictures in this lecture were legitimate all along.

import numpy as np

A = np.array([[2., 4., -2.],
              [4., 9., -3.],
              [-2., -3., 7.]])       # the elimination example from Lecture 4
r = lambda x: round(float(x), 9)     # det comes out of elimination: it rounds

print(r(np.linalg.det(A)))                     # 8.0
print(r(np.linalg.det(A[[1, 0, 2]])))     # -8.0   a row swap flips the sign

E = A.copy()
E[1] = E[1] - 2 * E[0]               # exactly what elimination does
E[2] = E[2] + E[0]
print(r(np.linalg.det(E)))                     # 8.0    unchanged

U = np.array([[2., 4., -2.], [0., 1., 1.], [0., 0., 4.]])
print(r(np.prod(np.diag(U))))            # 8.0    the pivots multiply to it

D = A.copy(); D[2] = D[0]            # make two rows equal
print(r(np.linalg.det(D)))                     # 0.0    flat

B = np.array([[0., 2., 1.], [1., 1., 0.], [3., 0., 1.]])
print(r(np.linalg.det(A @ B)))                 # -40.0
print(r(np.linalg.det(A) * np.linalg.det(B)))  # -40.0   the same
print(r(np.linalg.det(A.T)))                   # 8.0    same as det A

import torch

A = torch.tensor([[2., 4., -2.],
                  [4., 9., -3.],
                  [-2., -3., 7.]], dtype=torch.float64)   # from Lecture 4
r = lambda x: round(float(x), 9)

print(r(torch.linalg.det(A)))                  # 8.0
print(r(torch.linalg.det(A[[1, 0, 2]])))  # -8.0   a row swap flips the sign

E = A.clone()                        # clone, not copy
E[1] = E[1] - 2 * E[0]
E[2] = E[2] + E[0]
print(r(torch.linalg.det(E)))                  # 8.0    unchanged

urows = [[2., 4., -2.], [0., 1., 1.], [0., 0., 4.]]
U = torch.tensor(urows, dtype=torch.float64)
print(r(torch.prod(torch.diagonal(U))))  # 8.0    the pivots multiply to it

D = A.clone(); D[2] = D[0]
print(r(torch.linalg.det(D)))                  # 0.0    flat

brows = [[0., 2., 1.], [1., 1., 0.], [3., 0., 1.]]
B = torch.tensor(brows, dtype=torch.float64)
print(r(torch.linalg.det(A @ B)))                          # -40.0
print(r(torch.linalg.det(A) * torch.linalg.det(B)))       # -40.0   the same
print(r(torch.linalg.det(A.T)))                # 8.0    same as det A

import jax
jax.config.update("jax_enable_x64", True)
import jax.numpy as jnp

A = jnp.array([[2., 4., -2.],
               [4., 9., -3.],
               [-2., -3., 7.]])      # the elimination example from Lecture 4
r = lambda x: round(float(x), 9)

print(r(jnp.linalg.det(A)))                    # 8.0
print(r(jnp.linalg.det(A[jnp.array([1, 0, 2])])))  # -8.0   a row swap flips

E = A.at[1].add(-2 * A[0])           # jax arrays are immutable: .at[].add()
E = E.at[2].add(A[0])
print(r(jnp.linalg.det(E)))                    # 8.0    unchanged

U = jnp.array([[2., 4., -2.], [0., 1., 1.], [0., 0., 4.]])
print(r(jnp.prod(jnp.diag(U))))          # 8.0    the pivots multiply to it

D = A.at[2].set(A[0])                # make two rows equal
print(r(jnp.linalg.det(D)))                    # 0.0    flat

B = jnp.array([[0., 2., 1.], [1., 1., 0.], [3., 0., 1.]])
print(r(jnp.linalg.det(A @ B)))                   # -40.0
print(r(jnp.linalg.det(A) * jnp.linalg.det(B)))   # -40.0   the same
print(r(jnp.linalg.det(A.T)))                  # 8.0    same as det A

import tensorflow as tf

A = tf.constant([[2., 4., -2.],
                 [4., 9., -3.],
                 [-2., -3., 7.]], dtype=tf.float64)   # from Lecture 4
r = lambda x: round(float(x), 9)

print(r(tf.linalg.det(A)))                     # 8.0
print(r(tf.linalg.det(tf.gather(A, [1, 0, 2]))))  # -8.0   gather to reorder

E = tf.tensor_scatter_nd_update(A, [[1], [2]], [A[1] - 2 * A[0], A[2] + A[0]])
print(r(tf.linalg.det(E)))                     # 8.0    unchanged

U = tf.constant([[2., 4., -2.], [0., 1., 1.], [0., 0., 4.]], dtype=tf.float64)
print(r(tf.reduce_prod(tf.linalg.diag_part(U))))  # 8.0   the pivots multiply

D = tf.tensor_scatter_nd_update(A, [[2]], [A[0]])   # make two rows equal
print(r(tf.linalg.det(D)))                     # 0.0    flat

B = tf.constant([[0., 2., 1.], [1., 1., 0.], [3., 0., 1.]], dtype=tf.float64)
print(r(tf.linalg.det(A @ B)))                        # -40.0
print(r(tf.linalg.det(A) * tf.linalg.det(B)))         # -40.0   the same
print(r(tf.linalg.det(tf.transpose(A))))            # 8.0    same as det A

The determinant your computer actually computes

Put two of those consequences side by side. Elimination does not change the determinant, and a triangular matrix is the product of its diagonal. So run elimination, then multiply the pivots. That is the whole algorithm, and it is what every library does.

Take the elimination example from Lecture 4:

The pivots were 2, 1, 4 when you first met them, and their product is the determinant. If elimination needed row exchanges, count them: $PA = LU$ gives $\det A = \pm d_1 d_2 \cdots d_n$, with the sign fixed by whether the number of exchanges was even or odd.

The tridiagonal matrix from Lecture 5 makes the same point. Its pivots were $2$, $\tfrac{3}{2}$, $\tfrac{4}{3}$, and they multiply to 4:

One more reading of the same fact. The $k$-th pivot is a ratio of determinants:

where $A_k$ is the top-left $k \times k$ corner of $A$. The product telescopes, which is the formula above run backwards.

Cost: elimination is about $n^3/3$ operations, so a $1000 \times 1000$ determinant takes a fraction of a second. Hold that number. It is about to be worth a great deal.

Counting every term honestly

Pivots are excellent for computing and useless for seeing structure. They are ratios of things that came out of elimination, and you cannot look at them and see the original entries $a_{ij}$. So derive a formula that uses the entries directly, from the three rules alone.

Split each row into its coordinate pieces and expand by linearity. For a $3 \times 3$ matrix each row splits three ways, giving 27 terms, but any term that picks two entries from the same column has two proportional rows and dies. What survives is one entry from each row and one from each column. Count those: three choices for row one, two left for row two, one for row three. Six terms, or $n!$ in general.

Each surviving term is a signed product, and the sign comes from Rule 2. Write the column numbers in row order; the sign is $+$ if that ordering can be sorted back to $1, 2, 3$ with an even number of swaps, and $-$ if odd.

In general this is the big formula:

summing over all $n!$ orderings $\sigma$ of the columns. It is a real formula, it is exact, and nobody computes with it. Here is why, measured on one laptop:

import numpy as np, time
from itertools import permutations

def sign(p):
    s = 1
    for i in range(len(p)):
        for j in range(i + 1, len(p)):
            if p[i] > p[j]: s = -s
    return s

def big_formula(A):
    n = len(A)
    return sum(sign(p) * np.prod([A[i, p[i]] for i in range(n)])
               for p in permutations(range(n)))

rng = np.random.default_rng(0)
for n in [3, 8, 10, 11]:
    A = rng.standard_normal((n, n))
    t = time.perf_counter(); d1 = big_formula(A)
    big = time.perf_counter() - t
    t = time.perf_counter(); d2 = np.linalg.det(A)
    lu = time.perf_counter() - t
    print(n, f"{big:9.4f}s", f"{lu:.6f}s", abs(d1 - d2) < 1e-9)
#  3      0.0000s 0.000043s True
#  8      0.1253s 0.000023s True         40,320 terms
# 10     12.8782s 0.000030s True      3,628,800 terms
# 11    153.1656s 0.000115s True     39,916,800 terms, at 11 by 11

Two and a half minutes against a tenth of a millisecond, and the gap grows by a factor of $n$ every time you add a row. Your machine will print different seconds; the shape of the table is what matters. At $n = 20$ the big formula needs $2.4 \times 10^{18}$ terms, which at the rate measured above is about three hundred thousand years. Elimination on the same matrix needs about 2700 operations.

Cofactors, or the same sum folded

The six terms have a pattern. Group them by which entry they took from the first row:

Every bracket is a $2 \times 2$ determinant, and each one is what remains of $A$ after deleting row 1 and the column that $a_{1j}$ came from. That leftover determinant is the minor $M_{1j}$, and with its sign attached it is the cofactor:

The signs form a checkerboard starting with $+$ in the top-left corner. Expanding along any row gives $\det A = a_{i1}C_{i1} + a_{i2}C_{i2} + \cdots + a_{in}C_{in}$, which is worth doing by hand when a row is mostly zeros. For the matrix above, along the first row:

Same 8, third time.

Cofactors also give a formula for the inverse, which almost every course leaves out. Collect all $n^2$ cofactors, transpose the result, and you have the adjugate $\operatorname{adj}(A)$, whose $(i, j)$ entry is $C_{ji}$. Then

Test it on the tridiagonal $K$ from Lecture 5. Its determinant is 4, and its nine cofactors turn out to be

which is the inverse Gauss-Jordan produced back in Lecture 5, arrived at from the opposite direction. The fraction $1/\det A$ out front says the same thing again: no inverse when the determinant is zero.

The cross product is a determinant

Here is a formula from physics and graphics that turns out to be this lecture in disguise. For vectors $u$ and $v$ in three dimensions, the cross product is

with the standard basis vectors written into the top row and cofactor expansion applied along it. The result is a vector, perpendicular to both inputs, and its length is the area of the parallelogram they span. That last fact is not a coincidence. It is this lecture: the area of a parallelogram is a determinant.

The right-hand rule, that awkward business with your fingers, is the sign of the determinant. Swap $u$ and $v$ and you swap two rows, so $v \times u = -(u \times v)$ and the arrow points the other way.

Cramer’s rule, and its price

Determinants also solve $Ax = b$ outright. The idea is one trick. Take the identity matrix, replace its first column with the unknown $x$, and call the result $X_1$. That matrix is triangular with $x_1$ on the diagonal, so $\det X_1 = x_1$. Now multiply by $A$:

The first column became $Ax$, which is $b$. The other columns are copies of $A$. Take determinants of both sides and use the product rule:

Repeat with the unknown placed in column $i$ and you get Cramer’s rule: $x_i = \det B_i / \det A$, where $B_i$ is $A$ with column $i$ replaced by $b$. On the running example with $b = (2, 8, 10)$, the four determinants are $\det A = 8$ and $\det B_1 = -8$, $\det B_2 = 16$, $\det B_3 = 16$, so

which is the answer back substitution gave in Lecture 4.

Now the price. Cramer’s rule needs $n+1$ determinants of size $n$. Each one is an elimination, so the total is about $n+1$ times the cost of just solving the system, and that is the charitable accounting. There is a second problem that is worse and less often mentioned.

import numpy as np

def cramer(A, b):
    d = np.linalg.det(A)
    x = []
    for i in range(len(b)):
        B = A.copy()
        B[:, i] = b            # column i swapped for b
        x.append(np.linalg.det(B) / d)
    return np.array(x)

A = np.array([[2., 4., -2.], [4., 9., -3.], [-2., -3., 7.]])
b = np.array([2., 8., 10.])
print(cramer(A, b))            # [-1.  2.  2.]
print(np.linalg.solve(A, b))   # [-1.  2.  2.]   they agree

u, v = np.array([1., 2., 3.]), np.array([4., 5., 6.])
i, j = np.array([1., 0., 0.]), np.array([0., 1., 0.])
print(np.cross(u, v))          # [-3.  6. -3.]
print(np.cross(i, j))          # [0. 0. 1.]     i cross j is k

rng = np.random.default_rng(1)
big = rng.standard_normal((200, 200)) + 200 * np.eye(200)
print(np.linalg.det(big))      # inf, with an overflow warning: |det| is ~1e460

import torch

def cramer(A, b):
    d = torch.linalg.det(A)
    x = []
    for i in range(len(b)):
        B = A.clone()          # clone, not copy
        B[:, i] = b            # column i swapped for b
        x.append(torch.linalg.det(B) / d)
    return torch.stack(x)

rows = [[2., 4., -2.], [4., 9., -3.], [-2., -3., 7.]]
A = torch.tensor(rows, dtype=torch.float64)
b = torch.tensor([2., 8., 10.], dtype=torch.float64)
print(cramer(A, b))               # tensor([-1., 2., 2.])
print(torch.linalg.solve(A, b))   # tensor([-1., 2., 2.])   they agree

u, v = torch.tensor([1., 2., 3.]), torch.tensor([4., 5., 6.])
i, j = torch.tensor([1., 0., 0.]), torch.tensor([0., 1., 0.])
print(torch.linalg.cross(u, v))   # tensor([-3., 6., -3.])   linalg.cross
print(torch.linalg.cross(i, j))   # tensor([0., 0., 1.])

torch.manual_seed(1)
I200 = torch.eye(200, dtype=torch.float64)
big = torch.randn(200, 200, dtype=torch.float64) + 200 * I200
print(torch.linalg.det(big))      # inf, the determinant overflows float64

import jax
jax.config.update("jax_enable_x64", True)
import jax.numpy as jnp

def cramer(A, b):
    d = jnp.linalg.det(A)
    x = []
    for i in range(len(b)):
        B = A.at[:, i].set(b)  # jax is immutable: .at[].set()
        x.append(jnp.linalg.det(B) / d)
    return jnp.array(x)

A = jnp.array([[2., 4., -2.], [4., 9., -3.], [-2., -3., 7.]])
b = jnp.array([2., 8., 10.])
print(cramer(A, b))             # [-1.  2.  2.]
print(jnp.linalg.solve(A, b))   # [-1.  2.  2.]   they agree

u, v = jnp.array([1., 2., 3.]), jnp.array([4., 5., 6.])
i, j = jnp.array([1., 0., 0.]), jnp.array([0., 1., 0.])
print(jnp.cross(u, v))          # [-3.  6. -3.]
print(jnp.cross(i, j))          # [0. 0. 1.]     i cross j is k

key = jax.random.key(1)
big = jax.random.normal(key, (200, 200)) + 200 * jnp.eye(200)
print(jnp.linalg.det(big))      # inf, the determinant overflows float64

import tensorflow as tf

def cramer(A, b):
    d = tf.linalg.det(A)
    n = int(b.shape[0])
    x = []
    for i in range(n):
        cols = [b if k == i else A[:, k] for k in range(n)]
        x.append(tf.linalg.det(tf.stack(cols, axis=1)) / d)
    return tf.stack(x)

rows = [[2., 4., -2.], [4., 9., -3.], [-2., -3., 7.]]
A = tf.constant(rows, dtype=tf.float64)
b = tf.constant([2., 8., 10.], dtype=tf.float64)
print(cramer(A, b))               # [-1.  2.  2.]
rhs = tf.reshape(b, [3, 1])       # tf wants a 2-D right-hand side
print(tf.linalg.solve(A, rhs))    # [[-1.] [2.] [2.]]   they agree

u, v = [[1., 2., 3.]], [[4., 5., 6.]]
print(tf.linalg.cross(u, v))      # [[-3.  6. -3.]]   cross wants batches
print(tf.linalg.cross([[1., 0., 0.]], [[0., 1., 0.]]))   # [[0. 0. 1.]]

noise = tf.random.stateless_normal([200, 200], seed=[1, 1], dtype=tf.float64)
big = noise + 200 * tf.eye(200, dtype=tf.float64)
print(tf.linalg.det(big))         # inf, the determinant overflows float64

At size 200 the determinant of the matrix above is around $10^{460}$, and float64 stops at about $10^{308}$. So $\det A$ comes back as infinity, every $\det B_i$ comes back as infinity, and Cramer’s rule returns a vector of nan. Timed on the same matrix, it spends about a hundred times longer than solve to arrive at that answer. The elimination path never forms a determinant at all, so it never notices.

The receipt a flow model keeps

A normalizing flow is a stack of invertible maps. It starts from a simple distribution, pushes it through the stack, and lands on the shape of your data. To train it you need the density it assigns at the far end, and a density is an amount of probability per unit volume. So every layer has to declare how much volume it stretched or squeezed on the way through, and that declaration is the determinant of its Jacobian.

Watch two layers do it with real numbers. The first has Jacobian $J_1$ with rows $(1.5, 0.5)$ and $(0, 1)$, so $\det J_1 = 1.5$. The second has rows $(0.8, 0)$ and $(0.4, 0.5)$, so $\det J_2 = 0.4$.

The areas multiply, $1.5 \times 0.4 = 0.6$, because determinants multiply. Taking logs turns that product into a sum, which is why a flow reports $\log|\det J|$ per layer and adds:

Now you can see the shape of the whole field. A general $n \times n$ Jacobian costs $n^3/3$ just to get one determinant, per layer, per training step, which no one can afford. So flow architectures are designed backwards from this chapter: build layers whose Jacobian is triangular, where P7 says the determinant is the product of the diagonal and the whole receipt is a sum of $n$ logarithms. Coupling layers, autoregressive flows, and the triangular tricks behind them all exist to keep that one number cheap.

Where this is going

One number, and it has told you whether a matrix is invertible, how much it scales volume, whether it preserves handedness, and what its inverse is. The next lecture puts it to the use it was really built for.

Ask which directions a matrix leaves pointing the same way. If $Ax = \lambda x$ for some nonzero $x$, then $(A - \lambda I)x = 0$, so the matrix $A - \lambda I$ has a nonzero null space and must be singular. The test for singular is a determinant:

That is a polynomial in $\lambda$, its roots are the eigenvalues, and the whole of Lecture 15 follows from it. There is a second hint waiting there too. The big formula’s diagonal term makes the product of the eigenvalues come out equal to $\det A$, so the volume factor is exactly how much all the special directions stretch, multiplied together.