13. Advanced Topics and Applications

The previous chapters built a complete theory of vectors, linear maps, eigenvalues, and their computation. This chapter shows how far that theory reaches. It covers multilinear algebra (forms, tensors, wedge and Kronecker products), the spectral theory of nonnegative matrices and graphs, six applied domains where linear algebra does most of the work, and the branches of mathematics that start where a first course ends. The goal is to recognize the same few ideas (bases, eigenvectors, orthogonality, factorizations) in settings that look different at first.

<a id="131-multilinear-algebra"></a>

13.1 Multilinear Algebra

Linear algebra studies functions that are linear in one vector argument. Multilinear algebra studies functions that are linear in each of several arguments separately (inner products, determinants, stress tensors), and the vector spaces built to turn them back into ordinary linear maps.

Bilinear and Sesquilinear Forms

Why it matters. Many core objects in this guide take two vectors and return a scalar: the dot product (Section 2.3), general inner products (Section 9.1), and the quadratic forms $\mathbf{x}^TA\mathbf{x}$ of Section 10.2. Each is linear in each argument when the other is held fixed. Studying this property on its own, without assuming symmetry or positivity, gives a single framework for all of them. It also covers forms that are not inner products at all, such as the Minkowski form of special relativity and the symplectic form of classical mechanics.

Analogy. A linear functional is a machine with one input slot that turns a vector into a number. A bilinear form is a machine with two slots. Freeze either input, and what remains is an ordinary linear functional of the other.

Definition. Let $V$ be a vector space over a field $F$. A bilinear form is a map $B : V \times V \to F$ such that, for all vectors and scalars,

$$B(a\mathbf{u} + b\mathbf{v}, \mathbf{w}) = aB(\mathbf{u}, \mathbf{w}) + bB(\mathbf{v}, \mathbf{w}), \qquad B(\mathbf{u}, a\mathbf{v} + b\mathbf{w}) = aB(\mathbf{u}, \mathbf{v}) + bB(\mathbf{u}, \mathbf{w}).$$

Matrix representation. Fix a basis $\mathbf{e}_1, \dots, \mathbf{e}n$ and set $A{ij} = B(\mathbf{e}_i, \mathbf{e}_j)$. Bilinearity then gives, for coordinate vectors $\mathbf{x}$ and $\mathbf{y}$,

$$B(\mathbf{x}, \mathbf{y}) = \sum_{i,j} x_iy_jB(\mathbf{e}_i, \mathbf{e}_j) = \mathbf{x}^TA\mathbf{y}.$$

So bilinear forms on $F^n$ correspond one to one with $n \times n$ matrices. The matrix, however, depends on the basis, and it transforms differently from the matrix of a linear operator. If the new basis vectors are the columns of an invertible $P$, then $\mathbf{x} = P\mathbf{x}'$ and

$$B(\mathbf{x}, \mathbf{y}) = (P\mathbf{x}')^TA(P\mathbf{y}') = \mathbf{x}'^T(P^TAP)\mathbf{y}'.$$

	Linear operator $T : V \to V$	Bilinear form $B : V \times V \to F$
Matrix entry	$A_{ij}$ = $i$-th coordinate of $T\mathbf{e}_j$	$A_{ij} = B(\mathbf{e}_i, \mathbf{e}_j)$
Change of basis	similarity $P^{-1}AP$	congruence $P^TAP$
Invariants	eigenvalues, trace, determinant	rank, signature (for symmetric forms); $\det$ only up to a factor $\det(P)^2$
Canonical form	Jordan form (Section 8.3)	$\operatorname{diag}(1, \dots, 1, -1, \dots, -1, 0, \dots, 0)$ over $\mathbb{R}$ (Section 10.2)

When $P$ is orthogonal, $P^T = P^{-1}$ and the two transformations coincide. This is why symmetric matrices can be studied with either point of view.

Symmetry types. A form is symmetric if $B(\mathbf{u}, \mathbf{v}) = B(\mathbf{v}, \mathbf{u})$ (matrix $A = A^T$) and skew-symmetric if $B(\mathbf{u}, \mathbf{v}) = -B(\mathbf{v}, \mathbf{u})$ ($A = -A^T$). Away from characteristic 2, every form splits uniquely as the sum of the two: $A = \tfrac12(A + A^T) + \tfrac12(A - A^T)$. A skew-symmetric form satisfies $B(\mathbf{v}, \mathbf{v}) = 0$ for every $\mathbf{v}$, so the quadratic form $Q(\mathbf{x}) = B(\mathbf{x}, \mathbf{x})$ only sees the symmetric part. The symmetric form can be recovered from $Q$ by the polarization identity:

$$B(\mathbf{u}, \mathbf{v}) = \tfrac12\big(Q(\mathbf{u} + \mathbf{v}) - Q(\mathbf{u}) - Q(\mathbf{v})\big).$$

Nondegeneracy. The radical of a symmetric form is ${\mathbf{u} : B(\mathbf{u}, \mathbf{v}) = 0 \text{ for all } \mathbf{v}}$, which is the null space of $A$. The form is nondegenerate when the radical is ${\mathbf{0}}$, that is, when $A$ is invertible. A nondegenerate form identifies $V$ with its dual $V^*$ (Section 7.4), by sending $\mathbf{u}$ to the functional $B(\mathbf{u}, \cdot)$. Nondegenerate does not mean positive definite:

Form on $\mathbb{R}^n$	Matrix	Nondegenerate?	Positive definite?	Where it appears
Dot product	$I_n$	yes	yes	Euclidean geometry
Minkowski form	$\operatorname{diag}(-1, 1, 1, 1)$	yes	no (signature $(3, 1)$)	special relativity
Symplectic form	$\begin{bmatrix} 0 & I \ -I & 0 \end{bmatrix}$	yes	no (skew)	Hamiltonian mechanics
$Q(x, y) = x^2$	$\operatorname{diag}(1, 0)$	no (radical = $y$-axis)	no	degenerate conics

Sesquilinear forms. Over $\mathbb{C}$, bilinearity is the wrong requirement for a length. If $B$ were bilinear, then $B(i\mathbf{v}, i\mathbf{v}) = i^2B(\mathbf{v}, \mathbf{v}) = -B(\mathbf{v}, \mathbf{v})$, so no nonzero vector could have positive "squared length" for both $\mathbf{v}$ and $i\mathbf{v}$. The fix is to conjugate scalars in one slot. A sesquilinear form ("one and a half linear") is conjugate-linear in the first argument and linear in the second:

$$h(a\mathbf{u}, \mathbf{v}) = \bar{a},h(\mathbf{u}, \mathbf{v}), \qquad h(\mathbf{u}, a\mathbf{v}) = a,h(\mathbf{u}, \mathbf{v}).$$

Its matrix form is $h(\mathbf{x}, \mathbf{y}) = \mathbf{x}^*H\mathbf{y}$, and a change of basis gives $P^HP$. The form is Hermitian when $h(\mathbf{u}, \mathbf{v}) = \overline{h(\mathbf{v}, \mathbf{u})}$, equivalently $H = H^$. Then $h(\mathbf{v}, \mathbf{v})$ is always real, so positivity makes sense. A positive definite Hermitian form is exactly a complex inner product (Section 9.1).

Convention warning. Physics and most numerical libraries make the first slot conjugate-linear, as here (np.vdot conjugates its first argument). Many mathematics texts conjugate the second slot instead. The theory is identical; only the bookkeeping changes.

import numpy as np

# A bilinear form on R^2 given by a formula, and its matrix A_ij = B(e_i, e_j)
def B(x, y):
    return 2*x[0]*y[0] + 3*x[0]*y[1] - x[1]*y[0] + 4*x[1]*y[1]

E = np.eye(2)
A = np.array([[B(E[i], E[j]) for j in range(2)] for i in range(2)])
print(A)                                      # [[ 2.  3.] [-1.  4.]]
x, y = np.array([1.0, 2.0]), np.array([-3.0, 1.0])
print(B(x, y), x @ A @ y)                     # 11.0 11.0

# Split into symmetric and skew-symmetric parts
S, K = (A + A.T) / 2, (A - A.T) / 2
print(S, K, sep="\n")                         # [[2. 1.] [1. 4.]]  and  [[ 0.  2.] [-2.  0.]]
print(x @ K @ x)                              # 0.0: a skew form vanishes on the diagonal

# Congruence P^T S P preserves the signature (Sylvester), not the eigenvalues
M = np.array([[1.0, 2.0, 0.0], [2.0, -3.0, 1.0], [0.0, 1.0, 0.5]])
P = np.array([[1.0, 4.0, -2.0], [0.0, 1.0, 3.0], [2.0, 0.0, 1.0]])   # any invertible P
M2 = P.T @ M @ P
for X in (M, M2):
    ev = np.linalg.eigvalsh(X)
    print(np.round(ev, 3), (ev > 0).sum(), (ev < 0).sum())
# [-4.018  0.577  1.942] 2 1
# [-43.182   2.744  31.938] 2 1    different eigenvalues, same inertia (2 positive, 1 negative)

# A Hermitian sesquilinear form h(x, y) = x^* H y takes real values on the diagonal
H = np.array([[2.0, 1 - 1j], [1 + 1j, 3.0]])
z = np.array([1 + 2j, -1j])
w = np.array([0.5, 2 - 1j])
h = lambda x, y: np.conj(x) @ H @ y
print(h(z, z))                                # (7+0j): real
print(np.isclose(h(z, w), np.conj(h(w, z))))  # True: conjugate symmetry
print(np.isclose(h(1j * z, z), -1j * h(z, z)))   # True: conjugate-linear in the first slot

Key takeaways

• A bilinear form is linear in each argument separately; in a basis it is $\mathbf{x}^TA\mathbf{y}$ with $A_{ij} = B(\mathbf{e}_i, \mathbf{e}_j)$.
• Form matrices change by congruence $P^TAP$, not by similarity; eigenvalues are not invariant, but rank and signature are (Sylvester's law).
• Every form is the sum of a symmetric part (which determines the quadratic form) and a skew-symmetric part (which vanishes on the diagonal).
• Nondegenerate means $A$ is invertible, which is weaker than positive definite: the Minkowski and symplectic forms are nondegenerate but not definite.
• Over $\mathbb{C}$, use sesquilinear forms $\mathbf{x}^*H\mathbf{y}$; Hermitian ones take real values on the diagonal, and positive definite Hermitian forms are complex inner products.

🧪 Practice

1. Find the matrix of $B(\mathbf{x}, \mathbf{y}) = x_1y_2 - x_2y_1 + 3x_3y_3$ on $\mathbb{R}^3$. Split it into symmetric and skew-symmetric parts and decide whether $B$ is nondegenerate.
2. Show that $B(f, g) = \int_0^1 f(t)g(t),dt + f(0)g(1)$ is a bilinear form on polynomials of degree at most 1. Compute its matrix in the basis ${1, t}$. Is it symmetric?
3. Prove that $\det(P^TAP) = \det(P)^2\det(A)$, and conclude that over $\mathbb{R}$ the sign of the determinant of a symmetric form is basis independent.
4. Interview: Why does a complex inner product need a conjugate in one argument, and what goes wrong if you implement $\langle \mathbf{x}, \mathbf{y} \rangle$ as x @ y for complex arrays? Hint: evaluate $\langle \mathbf{x}, \mathbf{x} \rangle$ for $\mathbf{x} = (1, i)$, then compare with np.vdot.

Tensor Products of Vector Spaces

Why it matters. Bilinear maps are not linear maps, so none of the earlier machinery (kernels, matrices, rank) applies to them directly. The tensor product $V \otimes W$ is a new vector space built so that every bilinear map out of $V \times W$ becomes an ordinary linear map out of $V \otimes W$. The same construction describes composite quantum systems (Section 13.3), multi-index arrays in data science, and the Kronecker products used in numerical linear algebra.

Analogy. A recipe combines an amount of flour with an amount of sugar. The direct sum $V \oplus W$ keeps the two lists side by side: $m + n$ numbers. The tensor product records every pairwise interaction, like a multiplication table: $mn$ numbers, one for each pair of basis vectors.

The defining (universal) property. The tensor product of $V$ and $W$ is a vector space $V \otimes W$ together with a bilinear map $(\mathbf{v}, \mathbf{w}) \mapsto \mathbf{v} \otimes \mathbf{w}$ such that every bilinear map $\beta : V \times W \to U$ factors uniquely through a linear map $L$:

        V x W ----- (v, w) |-> v (x) w ----->  V (x) W
             \                                   |
              \  beta (bilinear)                 |  L (linear, unique)
               \                                 v
                `------------------------------> U

        beta(v, w) = L(v (x) w)   for all v, w

Concrete construction. If $\mathbf{e}_1, \dots, \mathbf{e}_m$ is a basis of $V$ and $\mathbf{f}_1, \dots, \mathbf{f}_n$ a basis of $W$, then the $mn$ symbols $\mathbf{e}_i \otimes \mathbf{f}_j$ form a basis of $V \otimes W$, so

$$\dim(V \otimes W) = \dim V \cdot \dim W.$$

The product of two vectors expands by bilinearity: $\mathbf{v} \otimes \mathbf{w} = \sum_{i,j} v_iw_j,(\mathbf{e}_i \otimes \mathbf{f}_j)$. Its coordinates $v_iw_j$ form the $m \times n$ matrix $\mathbf{v}\mathbf{w}^T$. This identifies $F^m \otimes F^n$ with the space of $m \times n$ matrices, with pure tensors $\mathbf{v} \otimes \mathbf{w}$ corresponding to outer products.

Computation rules. These follow from bilinearity:

• $(\mathbf{v}_1 + \mathbf{v}_2) \otimes \mathbf{w} = \mathbf{v}_1 \otimes \mathbf{w} + \mathbf{v}_2 \otimes \mathbf{w}$, and likewise in the second slot.
• $(a\mathbf{v}) \otimes \mathbf{w} = \mathbf{v} \otimes (a\mathbf{w}) = a(\mathbf{v} \otimes \mathbf{w})$. Scalars move freely between factors.
• $\mathbf{0} \otimes \mathbf{w} = \mathbf{v} \otimes \mathbf{0} = \mathbf{0}$.
• $\mathbf{v} \otimes \mathbf{w} \ne \mathbf{w} \otimes \mathbf{v}$ in general. When $V = W$, they are different elements of $V \otimes V$.

Pure tensors versus general tensors. An element of the form $\mathbf{v} \otimes \mathbf{w}$ is called pure (or simple, or decomposable). Most elements of $V \otimes W$ are not pure; they are sums of pure tensors. Under the matrix identification, pure tensors are exactly the rank-1 matrices, and the smallest number of pure tensors needed to write $T$ is the rank of the matrix $T$. The SVD (Section 11.1) writes any $T$ as a sum of $\operatorname{rank}(T)$ pure tensors $\sigma_k\mathbf{u}_k \otimes \mathbf{v}_k$. In quantum mechanics, non-pure tensors are called entangled states.

	Direct sum $V \oplus W$	Tensor product $V \otimes W$
Elements	pairs $(\mathbf{v}, \mathbf{w})$	sums $\sum_k \mathbf{v}_k \otimes \mathbf{w}_k$
Dimension	$m + n$	$mn$
Basis	${(\mathbf{e}_i, \mathbf{0})} \cup {(\mathbf{0}, \mathbf{f}_j)}$	${\mathbf{e}_i \otimes \mathbf{f}_j}$
Scaling	$a(\mathbf{v}, \mathbf{w}) = (a\mathbf{v}, a\mathbf{w})$	$a(\mathbf{v} \otimes \mathbf{w}) = (a\mathbf{v}) \otimes \mathbf{w}$
Linearizes	pairs of linear maps	bilinear maps
Physics	independent systems with additive degrees of freedom	composite quantum systems

Tensor products of maps. Linear maps $S : V \to V'$ and $R : W \to W'$ induce a linear map $S \otimes R : V \otimes W \to V' \otimes W'$, defined on pure tensors by $(S \otimes R)(\mathbf{v} \otimes \mathbf{w}) = S\mathbf{v} \otimes R\mathbf{w}$. In coordinates, its matrix is the Kronecker product of the two matrices (last topic of this subchapter). Two standard isomorphisms connect tensor products to earlier chapters:

$$\mathcal{L}(V, W) \cong V^* \otimes W, \qquad {\text{bilinear forms on } V \times W} \cong (V \otimes W)^* \cong V^* \otimes W^*.$$

The first says that a rank-1 linear map $\mathbf{x} \mapsto \varphi(\mathbf{x}),\mathbf{w}$ is the pure tensor $\varphi \otimes \mathbf{w}$, so a general matrix is a sum of such maps.

import numpy as np

v = np.array([1.0, 2.0])                      # v in V = R^2
w = np.array([3.0, 0.0, -1.0])                # w in W = R^3

# v (x) w  <->  the 2x3 matrix v w^T, coordinates in the basis e_i (x) f_j
T = np.outer(v, w)
print(T)                                      # [[ 3.  0. -1.] [ 6.  0. -2.]]
print(np.linalg.matrix_rank(T))               # 1: pure tensors are exactly the rank-1 matrices

# A general tensor is a sum of pure tensors; its tensor rank is the matrix rank
U = np.outer([1, 0], [1, 0, 0]) + np.outer([0, 1], [0, 1, 0])   # e1(x)f1 + e2(x)f2
print(np.linalg.matrix_rank(U))               # 2: not of the form v (x) w for any v, w

# Bilinearity: (2v) (x) w = v (x) (2w) = 2 (v (x) w)
print(np.allclose(np.outer(2 * v, w), np.outer(v, 2 * w)))   # True

# Universal property: the bilinear map beta(v, w) = v^T A w factors through v (x) w
# as the LINEAR functional  L(T) = sum_ij A_ij T_ij  (the Frobenius inner product <A, T>)
A = np.array([[1.0, -2.0, 0.5], [4.0, 0.0, 3.0]])
beta = v @ A @ w
L = lambda T: np.sum(A * T)                   # linear in T
print(beta, L(np.outer(v, w)))                # 2.5 2.5

# Linear maps act factor-wise: (S (x) R)(v (x) w) = Sv (x) Rw
S = np.array([[0.0, 1.0], [1.0, 1.0]])
R = np.array([[2.0, 0.0, 0.0], [0.0, 1.0, 1.0], [1.0, 0.0, 1.0]])
lhs = np.kron(S, R) @ np.kron(v, w)           # Kronecker product = matrix of S (x) R
rhs = np.kron(S @ v, R @ w)
print(np.allclose(lhs, rhs))                  # True

Key takeaways

• $V \otimes W$ is the space that turns bilinear maps out of $V \times W$ into linear maps; this universal property defines it up to isomorphism.
• A basis is ${\mathbf{e}_i \otimes \mathbf{f}_j}$, so dimensions multiply, unlike direct sums, where they add.
• $F^m \otimes F^n$ is the space of $m \times n$ matrices, with $\mathbf{v} \otimes \mathbf{w} = \mathbf{v}\mathbf{w}^T$; pure tensors are rank-1 matrices, and most tensors are not pure.
• $S \otimes R$ acts factor by factor, and its matrix is the Kronecker product.

🧪 Practice

1. In $\mathbb{R}^2 \otimes \mathbb{R}^2$, write $(\mathbf{e}_1 + 2\mathbf{e}_2) \otimes (3\mathbf{e}_1 - \mathbf{e}_2)$ in the basis $\mathbf{e}_i \otimes \mathbf{e}_j$.
2. Decide whether $\mathbf{e}_1 \otimes \mathbf{e}_1 + \mathbf{e}_1 \otimes \mathbf{e}_2 + \mathbf{e}_2 \otimes \mathbf{e}_1 + \mathbf{e}_2 \otimes \mathbf{e}_2$ is a pure tensor. If it is, find $\mathbf{v}$ and $\mathbf{w}$ with $T = \mathbf{v} \otimes \mathbf{w}$.
3. Let $P_2$ be the polynomials of degree at most 2. Show that $P_2 \otimes P_2$ is isomorphic to the space of polynomials $p(x, y)$ of degree at most 2 in each variable, via $p \otimes q \mapsto p(x)q(y)$. What is its dimension?
4. Interview: A dataset has $10^3$ users, $10^3$ items, and $10^2$ time steps. Stored densely as a tensor in $\mathbb{R}^{10^3} \otimes \mathbb{R}^{10^3} \otimes \mathbb{R}^{10^2}$ it has $10^8$ entries. How would you store and compute with it if you believe it is approximately low rank? Hint: sums of $r$ pure tensors need $r(m + n + p)$ numbers (a CP decomposition); compare with $mnp$.

Tensors and Index Notation

Why it matters. Physics and engineering describe quantities that need more than one index: the stress inside a material, the curvature of spacetime, the elasticity that links strain to stress. The question that matters is how their components change when the coordinate system changes, because physical laws must not depend on that choice. Index notation, with the Einstein summation convention, makes this bookkeeping almost mechanical, and it is exactly what np.einsum implements.

Analogy. A tensor is like a sentence and its components are a translation. Changing basis is changing language: the words (components) all change, but by fixed rules, so the meaning (the tensor itself) is the same.

Tensors of type $(p, q)$. A tensor of type $(p, q)$ on $V$ is an element of

$$\underbrace{V \otimes \cdots \otimes V}{p} \otimes \underbrace{V^* \otimes \cdots \otimes V^*}{q},$$

or equivalently a multilinear map that eats $q$ vectors and $p$ covectors (linear functionals) and returns a scalar. Its order is $p + q$, the number of indices. Common examples:

Type	Example	Components	Order
$(0, 0)$	scalar (temperature)	$s$	0
$(1, 0)$	vector (velocity)	$v^i$	1
$(0, 1)$	covector (gradient $df$, linear functional)	$a_i$	1
$(1, 1)$	linear operator $V \to V$	$A^i{}_j$	2
$(0, 2)$	bilinear form, metric tensor	$g_{ij}$	2
$(1, 3)$	Riemann curvature tensor	$R^i{}_{jkl}$	4

Terminology clash. In machine learning, "tensor" usually just means a multidimensional array, with no transformation rule attached, and "rank" often means the number of axes. In multilinear algebra, "rank" means the minimum number of pure tensors in a decomposition (previous topic). This guide uses order for the number of indices.

Upper and lower indices. Vector components carry upper indices and covector components carry lower ones. The placement records how the components respond to a change of basis. Suppose the new basis vectors are $\mathbf{e}'_j = P^i{}_j\mathbf{e}_i$ (the columns of $P$). Then:

• Contravariant (upper) components transform with the inverse: $v'^j = (P^{-1})^j{}_i,v^i$. If the basis vectors double in length, the coordinates halve.
• Covariant (lower) components transform like the basis itself: $a'_j = P^i{}_j,a_i$.
• A type-$(p, q)$ tensor gets one factor of $P^{-1}$ per upper index and one factor of $P$ per lower index. For a linear operator, $A' = P^{-1}AP$ (Section 7.2). For a bilinear form, $g' = P^TgP$ (previous topic).

Einstein summation convention. An index that appears once up and once down in the same term is summed over its range, and the summation sign is dropped:

$$w^i = A^i{}_jv^j ;\text{ means }; w^i = \sum_j A^i{}_jv^j, \qquad a_iv^i = \sum_i a_iv^i.$$

An index that is summed is a dummy index and can be renamed freely. An index that is not summed is a free index, and it must appear on both sides of an equation. A full contraction such as $a_iv^i$ has no free indices, so it is a scalar, and the transformation factors $P$ and $P^{-1}$ cancel. That is why it is basis independent.

Standard operations.

Operation	Index form	Matrix form	einsum
Matrix-vector product	$A^i{}_jv^j$	$A\mathbf{v}$	"ij,j->i"
Matrix product	$A^i{}_jB^j{}_k$	$AB$	"ij,jk->ik"
Trace (contraction)	$A^i{}_i$	$\operatorname{tr}A$	"ii->"
Outer product	$v^iw^j$	$\mathbf{v}\mathbf{w}^T$	"i,j->ij"
Lowering an index	$v_i = g_{ij}v^j$	$g\mathbf{v}$	"ij,j->i"
Bilinear form	$g_{ij}u^iv^j$	$\mathbf{u}^Tg\mathbf{v}$	"ij,i,j->"

Metric tensor. An inner product is a symmetric, positive definite $(0, 2)$ tensor $g_{ij}$. It converts vectors to covectors ($v_i = g_{ij}v^j$, lowering an index), and its inverse $g^{ij}$ converts back (raising an index). In an orthonormal basis, $g_{ij} = \delta_{ij}$ and upper and lower components coincide, which is why introductory courses never distinguish them.

Two special symbols. The Kronecker delta $\delta^i{}j$ (1 if $i = j$, else 0) is the identity tensor. The Levi-Civita symbol $\varepsilon{ijk}$ is $+1$ for even permutations of $(1, 2, 3)$, $-1$ for odd ones, and $0$ if an index repeats. With it, $(\mathbf{u} \times \mathbf{w})i = \varepsilon{ijk}u^jw^k$ and $\det M = \varepsilon_{ijk}M_{1i}M_{2j}M_{3k}$ (Section 5.2). Here indices are all written low because the basis is orthonormal.

import numpy as np

rng = np.random.default_rng(0)
A, B = rng.standard_normal((3, 4)), rng.standard_normal((4, 2))
v = rng.standard_normal(4)

# einsum is Einstein summation: repeated indices are summed, free indices remain
print(np.allclose(np.einsum("ij,jk->ik", A, B), A @ B))      # True: C^i_k = A^i_j B^j_k
print(np.allclose(np.einsum("ij,j->i", A, v), A @ v))        # True: w^i = A^i_j v^j
M = rng.standard_normal((4, 4))
print(np.isclose(np.einsum("ii->", M), np.trace(M)))         # True: contraction T^i_i

# Transformation laws. New basis vectors are the columns of P: e'_j = P^i_j e_i
P = np.array([[2.0, 1.0], [1.0, 1.0]])
v = np.array([3.0, 5.0])                      # vector components v^i (contravariant)
a = np.array([1.0, -2.0])                     # covector components a_i (covariant)
v_new = np.linalg.solve(P, v)                 # v'^j = (P^{-1})^j_i v^i
a_new = P.T @ a                               # a'_j = P^i_j a_i
print(v_new, a_new)                           # [-2.  7.] [ 0. -1.]
print(a @ v, a_new @ v_new)                   # -7.0 -7.0: the full contraction a_i v^i is invariant

# Metric tensor: raising/lowering indices. With g = I they coincide; with another g they do not
g = np.array([[2.0, 0.0], [0.0, 3.0]])        # a non-standard inner product on R^2
v_lower = g @ v                               # v_i = g_ij v^j
print(v_lower, v @ g @ v)                     # [ 6. 15.] 93.0: ||v||^2 = g_ij v^i v^j

# Levi-Civita symbol: (u x w)_i = eps_ijk u_j w_k, det(M) = eps_ijk M_1i M_2j M_3k
eps = np.zeros((3, 3, 3))
for (i, j, k), s in {(0, 1, 2): 1, (1, 2, 0): 1, (2, 0, 1): 1,
                     (0, 2, 1): -1, (2, 1, 0): -1, (1, 0, 2): -1}.items():
    eps[i, j, k] = s
u, w = np.array([1.0, 2.0, 3.0]), np.array([4.0, 5.0, 6.0])
print(np.einsum("ijk,j,k->i", eps, u, w), np.cross(u, w))   # [-3.  6. -3.] [-3.  6. -3.]
M3 = rng.standard_normal((3, 3))
print(np.isclose(np.einsum("ijk,i,j,k->", eps, M3[0], M3[1], M3[2]), np.linalg.det(M3)))  # True

Key takeaways

• A type-$(p, q)$ tensor has $p$ upper (contravariant) and $q$ lower (covariant) indices; each index fixes how components transform under a change of basis.
• Einstein convention: a repeated upper-lower index pair is summed; free indices must match on both sides of an equation.
• Full contractions (no free indices) are scalars and are invariant under change of basis.
• The metric tensor raises and lowers indices; in orthonormal coordinates the distinction disappears.
• np.einsum is index notation as code: products, traces, contractions, and transposes in one call.

🧪 Practice

1. Write in index notation, then as an einsum string: (a) $\operatorname{tr}(AB)$; (b) $\mathbf{u}^TA\mathbf{v}$; (c) the Frobenius inner product $\sum_{ij}A_{ij}B_{ij}$.
2. Show that $\operatorname{tr}(A) = A^i{}_i$ is invariant under $A \mapsto P^{-1}AP$ by writing the transformed components in index form and contracting.
3. Using $\varepsilon_{ijk}\varepsilon_{ilm} = \delta_{jl}\delta_{km} - \delta_{jm}\delta_{kl}$, prove $\mathbf{a} \times (\mathbf{b} \times \mathbf{c}) = \mathbf{b}(\mathbf{a} \cdot \mathbf{c}) - \mathbf{c}(\mathbf{a} \cdot \mathbf{b})$.
4. Interview: Given a batch of matrices A of shape (b, n, m) and B of shape (b, m, k), write one einsum call that computes the batched products, and another that computes the trace of each product without forming it. Hint: keep b as a free index in both; for the trace, contract the outer indices as well as m.

Exterior Algebra and Wedge Products

Why it matters. The cross product works only in $\mathbb{R}^3$, and the determinant was introduced in Chapter 5 through formulas whose geometric meaning came later. The wedge product handles both at once, in every dimension. It measures oriented area, volume, and their higher-dimensional analogs, and it gives a definition of the determinant that makes $\det(AB) = \det A\det B$ almost automatic. It is also the algebra behind differential forms and integration on curved spaces (Section 13.4).

Analogy. Think of $\mathbf{u} \wedge \mathbf{w}$ as the parallelogram spanned by $\mathbf{u}$ and $\mathbf{w}$, remembered only by its area, its plane, and its orientation (a direction of circulation). Swapping the two edges reverses the orientation, and a parallelogram with two parallel edges is flat and has area zero.

Definition. The wedge product is a bilinear product that is alternating:

$$\mathbf{v} \wedge \mathbf{v} = 0 \quad \text{for every } \mathbf{v}.$$

Expanding $(\mathbf{u} + \mathbf{w}) \wedge (\mathbf{u} + \mathbf{w}) = 0$ gives anticommutativity, $\mathbf{u} \wedge \mathbf{w} = -\mathbf{w} \wedge \mathbf{u}$. Products of $k$ vectors span the $k$-th exterior power $\Lambda^kV$. If $\mathbf{e}_1, \dots, \mathbf{e}_n$ is a basis of $V$, then

$${\mathbf{e}{i_1} \wedge \cdots \wedge \mathbf{e}{i_k} : i_1 < i_2 < \cdots < i_k} \text{ is a basis of } \Lambda^kV, \qquad \dim\Lambda^kV = \binom{n}{k}.$$

In particular $\Lambda^nV$ is one-dimensional, and $\Lambda^kV = 0$ for $k > n$.

$k$	$\Lambda^k\mathbb{R}^3$ basis	Dimension	Geometric meaning
0	$1$	1	scalars
1	$\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3$	3	oriented lengths (vectors)
2	$\mathbf{e}_1 \wedge \mathbf{e}_2, \mathbf{e}_1 \wedge \mathbf{e}_3, \mathbf{e}_2 \wedge \mathbf{e}_3$	3	oriented areas
3	$\mathbf{e}_1 \wedge \mathbf{e}_2 \wedge \mathbf{e}_3$	1	oriented volumes

Components are minors. For $\mathbf{u}, \mathbf{w} \in \mathbb{R}^n$,

$$\mathbf{u} \wedge \mathbf{w} = \sum_{i < j}(u_iw_j - u_jw_i),\mathbf{e}_i \wedge \mathbf{e}_j.$$

Each coefficient is a $2 \times 2$ minor of the $n \times 2$ matrix $[\mathbf{u}\ \mathbf{w}]$. More generally, the coefficients of $\mathbf{v}_1 \wedge \cdots \wedge \mathbf{v}_k$ are the $k \times k$ minors of $[\mathbf{v}_1 \cdots \mathbf{v}_k]$. In $\mathbb{R}^3$, the three components of $\mathbf{u} \wedge \mathbf{w}$ are, up to order and sign, the components of $\mathbf{u} \times \mathbf{w}$. The cross product is the wedge product, converted back into a vector by the normal direction. That conversion (the Hodge star) only produces a vector when $n = 3$, which is why the cross product exists only in three dimensions.

The determinant, intrinsically. In top degree,

$$\mathbf{v}_1 \wedge \cdots \wedge \mathbf{v}_n = \det[\mathbf{v}_1 \cdots \mathbf{v}_n];\mathbf{e}_1 \wedge \cdots \wedge \mathbf{e}_n.$$

A linear map $T : V \to V$ induces $\Lambda^nT$ on the one-dimensional space $\Lambda^nV$, via $\mathbf{v}_1 \wedge \cdots \wedge \mathbf{v}_n \mapsto T\mathbf{v}_1 \wedge \cdots \wedge T\mathbf{v}_n$. Any linear map on a one-dimensional space is multiplication by a scalar, and that scalar is $\det T$. This definition needs no basis. It also proves the product rule in one line: $\Lambda^n(ST) = \Lambda^nS \circ \Lambda^nT$, so $\det(ST) = \det S,\det T$.

Main facts.

• Independence test: $\mathbf{v}_1, \dots, \mathbf{v}_k$ are linearly independent if and only if $\mathbf{v}_1 \wedge \cdots \wedge \mathbf{v}_k \ne 0$.
• Subspaces: if ${\mathbf{v}_i}$ and ${\mathbf{w}_i}$ span the same $k$-dimensional subspace, their wedges differ only by a nonzero scalar (the determinant of the change of basis). The minors therefore encode the subspace itself; they are its Plücker coordinates.
• Volume: the $k$-dimensional volume of the parallelotope spanned by $\mathbf{v}_1, \dots, \mathbf{v}_k$ is the Euclidean length of the coefficient vector. By the Cauchy-Binet formula, its square equals the Gram determinant $\det(V^TV)$.

import numpy as np
from itertools import combinations

def wedge(*vectors):
    """Components of v1 ^ ... ^ vk in the basis e_I, I = (i1 < ... < ik): the k x k minors."""
    V = np.column_stack(vectors)              # n x k
    n, k = V.shape
    return {I: np.linalg.det(V[list(I), :]) for I in combinations(range(n), k)}

u = np.array([1.0, 2.0, 3.0])
w = np.array([4.0, 5.0, 6.0])
print({I: float(round(c, 6)) for I, c in wedge(u, w).items()})
# {(0, 1): -3.0, (0, 2): -6.0, (1, 2): -3.0}     coefficients of e1^e2, e1^e3, e2^e3
print(np.cross(u, w))                         # [-3.  6. -3.]: the same numbers, reordered as
                                              # (e2^e3, e3^e1, e1^e2) = (-3, -(-6), -3)

# Antisymmetry: w ^ u = -(u ^ w), and u ^ u = 0
print(all(np.isclose(wedge(w, u)[I], -c) for I, c in wedge(u, w).items()))   # True
print(max(abs(c) for c in wedge(u, u).values()))                          # 0.0

# Independence test: v1 ^ ... ^ vk != 0  iff  the vectors are independent
x = u + 2 * w                                 # dependent on u, w
print(max(abs(c) for c in wedge(u, w, x).values()))                       # 0.0

# Top degree: v1 ^ v2 ^ v3 = det[v1 v2 v3] e1^e2^e3
z = np.array([0.0, 1.0, 1.0])
print({I: float(round(c, 6)) for I, c in wedge(u, w, z).items()})          # {(0, 1, 2): 3.0}
print(round(np.linalg.det(np.column_stack([u, w, z])), 6))                # 3.0

# Squared k-volume = sum of squared minors = Gram determinant (Cauchy-Binet)
G = np.column_stack([u, w]).T @ np.column_stack([u, w])
print(sum(c**2 for c in wedge(u, w).values()), round(np.linalg.det(G), 6))   # 54.0 54.0: area = sqrt(54)

Key takeaways

• The wedge product is bilinear and alternating: $\mathbf{v} \wedge \mathbf{v} = 0$ and $\mathbf{u} \wedge \mathbf{w} = -\mathbf{w} \wedge \mathbf{u}$.
• $\Lambda^kV$ has dimension $\binom{n}{k}$; the components of $\mathbf{v}_1 \wedge \cdots \wedge \mathbf{v}_k$ are the $k \times k$ minors of $[\mathbf{v}_1 \cdots \mathbf{v}_k]$.
• A wedge product is nonzero exactly when its factors are independent; its length is the $k$-volume they span.
• In top degree, $\Lambda^nT$ is multiplication by $\det T$: a basis-free definition of the determinant that makes the product rule immediate.
• The cross product is the $n = 3$ special case of $\mathbf{u} \wedge \mathbf{w}$, turned into a vector.

🧪 Practice

1. Expand $(\mathbf{e}_1 + 2\mathbf{e}_2) \wedge (3\mathbf{e}_1 - \mathbf{e}_2)$ in $\Lambda^2\mathbb{R}^2$ and check that the coefficient equals $\det\begin{bmatrix} 1 & 3 \ 2 & -1 \end{bmatrix}$.
2. Compute the area of the parallelogram in $\mathbb{R}^4$ spanned by $(1, 0, 1, 0)$ and $(0, 1, 1, 1)$ in two ways: from the six $2 \times 2$ minors, and from the Gram determinant.
3. Show that $\mathbf{u} \wedge \mathbf{w} = 0$ if and only if $\mathbf{u}$ and $\mathbf{w}$ are parallel, directly from the minor formula.
4. Interview: Explain, without formulas for the determinant, why $\det(AB) = \det(A)\det(B)$ and why $\det A = 0$ exactly when $A$ is singular. Hint: consider how $A$ acts on $\mathbf{e}_1 \wedge \cdots \wedge \mathbf{e}_n$, and use the independence test.

Kronecker Product and Vectorization

Why it matters. Many problems have matrix-valued unknowns: the Sylvester equation $AX + XB = C$ in control theory, the Lyapunov equation for stability, image deblurring $AXB^T = C$, and PDEs discretized on 2D grids. Vectorization turns such a matrix equation into an ordinary system $K\mathbf{x} = \mathbf{c}$, and the Kronecker product describes the matrix $K$. The structure of $K$ also explains why such problems are cheap to solve when they are not treated as a generic $nm \times nm$ system.

Definition. For $A \in F^{m \times n}$ and $B \in F^{p \times q}$, the Kronecker product is the $mp \times nq$ block matrix

$$A \otimes B = \begin{bmatrix} a_{11}B & \cdots & a_{1n}B \ \vdots & \ddots & \vdots \ a_{m1}B & \cdots & a_{mn}B \end{bmatrix}.$$

It is the matrix of the tensor product of maps from the earlier topic, in the basis $\mathbf{e}_i \otimes \mathbf{f}_j$ ordered lexicographically. Vectorization $\operatorname{vec}(X)$ stacks the columns of $X$ into one long vector. The identity that connects them is

$$\boxed{\operatorname{vec}(AXB) = (B^T \otimes A)\operatorname{vec}(X).}$$

To see where it comes from, look at column $j$ of $AXB$: it is $\sum_k b_{kj},A\mathbf{x}_k$, where $\mathbf{x}_k$ is column $k$ of $X$. This is block row $j$ of $(B^T \otimes A)\operatorname{vec}(X)$.

Algebraic rules.

Rule	Statement
Mixed product	$(A \otimes B)(C \otimes D) = AC \otimes BD$ (when the products are defined)
Transpose, inverse	$(A \otimes B)^T = A^T \otimes B^T$, $(A \otimes B)^{-1} = A^{-1} \otimes B^{-1}$
Eigenvalues	the eigenvalues of $A \otimes B$ are all products $\lambda_i\mu_j$, with eigenvectors $\mathbf{x}_i \otimes \mathbf{y}_j$
Kronecker sum	$A \oplus B = A \otimes I_m + I_n \otimes B$ has eigenvalues $\lambda_i + \mu_j$
Determinant, trace	$\det(A \otimes B) = \det(A)^m\det(B)^n$ for $A$ $n \times n$ and $B$ $m \times m$; $\operatorname{tr}(A \otimes B) = \operatorname{tr}A,\operatorname{tr}B$
Rank, singular values	$\operatorname{rank}(A \otimes B) = \operatorname{rank}A \cdot \operatorname{rank}B$; singular values are $\sigma_i(A)\sigma_j(B)$
Not commutative	$A \otimes B \ne B \otimes A$ in general, but they are permutation-similar

The eigenvalue rule follows from the mixed product: $(A \otimes B)(\mathbf{x} \otimes \mathbf{y}) = A\mathbf{x} \otimes B\mathbf{y} = \lambda\mu,(\mathbf{x} \otimes \mathbf{y})$.

Matrix equations as linear systems.

Equation	Vectorized system
$AXB = C$	$(B^T \otimes A)\operatorname{vec}X = \operatorname{vec}C$
Sylvester: $AX + XB = C$	$(I \otimes A + B^T \otimes I)\operatorname{vec}X = \operatorname{vec}C$
Lyapunov: $AX + XA^T = -Q$	$(I \otimes A + A \otimes I)\operatorname{vec}X = -\operatorname{vec}Q$

By the Kronecker sum rule, the Sylvester equation has a unique solution exactly when no eigenvalue of $A$ is the negative of an eigenvalue of $B$.

Never form the big matrix. For $n \times n$ factors, $A \otimes B$ has $n^4$ entries. Solving $(B^T \otimes A)\operatorname{vec}X = \operatorname{vec}C$ densely costs $O(n^6)$, while solving $AXB = C$ as $X = A^{-1}CB^{-1}$ costs $O(n^3)$. The Kronecker form is a tool for reasoning (solvability, eigenvalues, conditioning); the computation should use the factors. The Bartels-Stewart algorithm (scipy.linalg.solve_sylvester) solves Sylvester equations in $O(n^3 + m^3)$ using Schur decompositions of $A$ and $B$.

Grids. The 5-point Laplacian on a $k \times k$ grid is the Kronecker sum $T \otimes I + I \otimes T$ of 1D second-difference matrices $T$. Its eigenvalues are all the sums of pairs of 1D eigenvalues, which is the basis of fast Poisson solvers.

NumPy note. NumPy stores arrays row by row, so X.ravel() stacks rows. Use X.reshape(-1, order="F") for the column-stacking vec, or the identity becomes $\operatorname{vec}{\text{row}}(AXB) = (A \otimes B^T)\operatorname{vec}{\text{row}}(X)$.

import numpy as np
from scipy.linalg import solve_sylvester

A = np.array([[1, 2], [3, 4]])
B = np.array([[0, 1], [1, 0]])
print(np.kron(A, B))
# [[0 1 0 2]
#  [1 0 2 0]
#  [0 3 0 4]
#  [3 0 4 0]]      block (i, j) is a_ij * B

vec = lambda X: X.reshape(-1, order="F")      # stack COLUMNS (NumPy is row-major by default)

rng = np.random.default_rng(1)
A, X, B = rng.standard_normal((3, 4)), rng.standard_normal((4, 5)), rng.standard_normal((5, 2))
print(np.allclose(vec(A @ X @ B), np.kron(B.T, A) @ vec(X)))   # True: vec(AXB) = (B^T (x) A) vec(X)

# Mixed-product rule and the spectrum of a Kronecker product
A = np.array([[1.0, 1.0], [0.0, 2.0]])             # eigenvalues 1, 2
B = np.array([[3.0, 1.0], [0.0, 5.0]])             # eigenvalues 3, 5
print(np.sort(np.linalg.eigvals(np.kron(A, B)).real))          # [ 3.  5.  6. 10.]: all products lambda_i mu_j

# Sylvester equation AX + XB = C  <=>  (I (x) A + B^T (x) I) vec(X) = vec(C)
n, m = 4, 3
A = rng.standard_normal((n, n)) + 4 * np.eye(n)
B = rng.standard_normal((m, m)) + 4 * np.eye(m)
C = rng.standard_normal((n, m))
K = np.kron(np.eye(m), A) + np.kron(B.T, np.eye(n))            # nm x nm: fine for tiny sizes only
X_kron = np.linalg.solve(K, vec(C)).reshape((n, m), order="F")
X_bs = solve_sylvester(A, B, C)               # Bartels-Stewart: O(n^3 + m^3), never forms K
print(np.allclose(X_kron, X_bs), np.allclose(A @ X_bs + X_bs @ B, C))   # True True

# The 2D grid Laplacian is a Kronecker sum of 1D Laplacians
k = 4
T = 2 * np.eye(k) - np.eye(k, k=1) - np.eye(k, k=-1)           # 1D second difference
L2 = np.kron(T, np.eye(k)) + np.kron(np.eye(k), T)             # 16 x 16, acts on vec(U) for a 4x4 grid
lam = np.linalg.eigvalsh(T)
pred = np.sort((lam[:, None] + lam[None, :]).ravel())          # all sums lambda_i + lambda_j
print(np.allclose(np.linalg.eigvalsh(L2), pred))               # True

Key takeaways

• $A \otimes B$ replaces each entry $a_{ij}$ with the block $a_{ij}B$; it is the coordinate matrix of a tensor product of maps.
• $\operatorname{vec}(AXB) = (B^T \otimes A)\operatorname{vec}(X)$ turns linear matrix equations into ordinary linear systems.
• Mixed product $(A \otimes B)(C \otimes D) = AC \otimes BD$; eigenvalues multiply under $\otimes$ and add under the Kronecker sum.
• Use Kronecker structure for analysis, not computation: forming $A \otimes B$ explicitly wastes memory and turns $O(n^3)$ work into $O(n^6)$.
• vec stacks columns; with NumPy's row-major default, pass order="F".

🧪 Practice

1. Compute $I_2 \otimes \begin{bmatrix} 1 & 2 \ 3 & 4 \end{bmatrix}$ and $\begin{bmatrix} 1 & 2 \ 3 & 4 \end{bmatrix} \otimes I_2$ by hand. Which one is block diagonal?
2. Prove $(A \otimes B)^{-1} = A^{-1} \otimes B^{-1}$ using the mixed-product rule.
3. Solve the Lyapunov equation $AX + XA^T = -I$ for $A = \begin{bmatrix} -1 & 1 \ 0 & -2 \end{bmatrix}$ by vectorizing. Verify that $X$ is symmetric positive definite, as Lyapunov theory predicts for a stable $A$.
4. Interview: You must apply $K = A \otimes B$, with $A$ and $B$ both $1000 \times 1000$, to a vector $\mathbf{x}$. How do you do it without forming $K$, and what does it cost? Hint: reshape $\mathbf{x}$ into a matrix $X$ and use the vec identity; compare $O(n^3)$ with $O(n^4)$ memory and time.

<a id="132-positive-matrices-and-graph-theory"></a>

13.2 Positive Matrices and Graph Theory

Matrices whose entries are counts, probabilities, or connection strengths have no negative entries, and this single sign condition gives them a much more rigid spectral theory. This subchapter develops that theory and applies it to graphs, where the eigenvalues and eigenvectors of a few matrices encode connectivity, centrality, and cluster structure.

Nonnegative Matrices

Why it matters. Transition probabilities of a Markov chain, Leontief input-output coefficients (Section 3.5), birth and survival rates in population models, and link counts between web pages cannot be negative. For general matrices, eigenvalues can be complex, the dominant one can be non-unique, and eigenvectors can have entries of any sign. Nonnegativity rules most of that out. The key insight of this topic is that the behavior of a nonnegative matrix is decided by its pattern of zero and nonzero entries, which is a directed graph.

Notation, and a warning. Write $A \ge 0$ if every entry $a_{ij} \ge 0$ (a nonnegative matrix) and $A > 0$ if every entry is strictly positive (a positive matrix). These are entrywise conditions and have nothing to do with definiteness:

	Positive matrix $A > 0$	Positive definite matrix $A \succ 0$
Condition	every entry $a_{ij} > 0$	$\mathbf{x}^TA\mathbf{x} > 0$ for all $\mathbf{x} \ne \mathbf{0}$
Requires symmetry?	no	yes (by convention)
Example	$\begin{bmatrix} 1 & 2 \ 3 & 1 \end{bmatrix}$ (eigenvalues $1 \pm \sqrt{6}$, one negative)	$\begin{bmatrix} 2 & -1 \ -1 & 2 \end{bmatrix}$ (has negative entries)
Theory	Perron-Frobenius (next topic)	spectral theorem, Cholesky (Chapter 10)

The graph of a matrix. An $n \times n$ matrix $A \ge 0$ defines a directed graph on vertices $1, \dots, n$, with an edge $i \to j$ whenever $a_{ij} > 0$. Matrix powers count weighted walks:

$$(A^k){ij} = \sum{i_1, \dots, i_{k-1}} a_{ii_1}a_{i_1i_2}\cdots a_{i_{k-1}j},$$

and every term is nonnegative, so $(A^k)_{ij} > 0$ exactly when the graph has a walk of length $k$ from $i$ to $j$. No cancellation can occur, which is why the pattern alone decides positivity.

Irreducible matrices. $A \ge 0$ is irreducible if its graph is strongly connected: every vertex can reach every other vertex by some walk. Equivalent conditions:

• $(I + A)^{n-1} > 0$.
• There is no permutation matrix $P$ that puts $A$ in block upper triangular form $P^TAP = \begin{bmatrix} B & C \ 0 & D \end{bmatrix}$ with square diagonal blocks.

A reducible matrix has a set of states that, once entered, can never be left. The problem then splits into smaller problems, one per strongly connected component, as in the block triangular systems of Section 4.5.

Primitive matrices and the period. $A \ge 0$ is primitive if $A^k > 0$ for some $k$: for a long enough length, there is a walk of exactly that length between every pair of vertices. Primitive implies irreducible, but not conversely. The obstruction is periodicity. The period $h$ of an irreducible matrix is the greatest common divisor of the lengths of all its closed walks. Then:

• $h = 1$ (aperiodic) $\iff$ primitive. A single self-loop ($a_{ii} > 0$) is enough to force $h = 1$.
• $h > 1$: the vertices split into $h$ classes visited in cyclic order, like the hands of a clock, so $A^k$ is never positive.
• If $A$ is primitive, then $A^k > 0$ already for $k = (n-1)^2 + 1$ (Wielandt's bound), so primitivity can be tested with finitely many multiplications.

    irreducible, period 3      primitive (period 1)       reducible

                                  +--+
                                  v  |  self-loop
    1 ----> 2                     1 ----> 2               1 ----> 2
    ^       |                     ^       |
    |       v                     |       v               2 can never
    +------ 3                     +------ 3               return to 1

    closed walks: 3, 6, 9, ...    closed walks: 1, 3, 4, ...
    gcd = 3                       gcd = 1

Spectral radius bounds. The spectral radius $\rho(A) = \max\lvert\lambda_i\rvert$ (Section 11.3) of a nonnegative matrix is easy to bound:

• $\min_i \sum_j a_{ij} \le \rho(A) \le \max_i \sum_j a_{ij}$ (row sums), and likewise for column sums. In particular, a (row-)stochastic matrix (nonnegative, row sums 1) has $\rho = 1$.
• Monotonicity: $0 \le A \le B$ entrywise implies $\rho(A) \le \rho(B)$. Increasing any entry can never decrease the spectral radius.
• Neumann series: if $\rho(A) < 1$, then $(I - A)^{-1} = I + A + A^2 + \cdots \ge 0$. This is why a productive Leontief economy has a nonnegative production vector for every nonnegative demand.

import numpy as np
from math import gcd

def is_irreducible(A):
    n = len(A)
    R = np.linalg.matrix_power(np.eye(n) + (A > 0), n - 1)   # (I + A)^(n-1) > 0  iff  strongly connected
    return bool((R > 0).all())

def is_primitive(A):
    n = len(A)
    Ak = np.linalg.matrix_power((A > 0).astype(float), (n - 1) ** 2 + 1)   # Wielandt's exponent bound
    return bool((Ak > 0).all())

def period(A):
    """gcd of the lengths of all closed walks through vertex 0 (for an irreducible A)."""
    n, g = len(A), 0
    P, M = np.eye(n), (A > 0).astype(float)
    for k in range(1, 2 * n + 1):
        P = np.minimum(P @ M, 1)              # keep only the 0/1 pattern, avoid overflow
        if P[0, 0] > 0:
            g = gcd(g, k)
    return g

cycle = np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]])          # 1 -> 2 -> 3 -> 1
cycle_plus = cycle + np.diag([1, 0, 0])                       # add a self-loop at vertex 1
reducible = np.array([[1, 1], [0, 1]])                        # 2 cannot reach 1
for name, A in [("cycle", cycle), ("cycle+loop", cycle_plus), ("reducible", reducible)]:
    print(name, is_irreducible(A), is_primitive(A), period(A) if is_irreducible(A) else "-")
# cycle True False 3
# cycle+loop True True 1
# reducible False False -

# Row-sum bounds on the spectral radius, and monotonicity
A = np.array([[0.2, 0.3, 0.1], [0.4, 0.1, 0.2], [0.1, 0.5, 0.3]])
rho = max(abs(np.linalg.eigvals(A)))
rows = np.round(A.sum(axis=1), 10)
print(rows.min(), round(rho, 4), rows.max())                       # 0.6 0.7186 0.9
print(max(abs(np.linalg.eigvals(A + 0.05))) > rho)                 # True: larger entries, larger rho

# Leontief model (Section 3.5): rho(A) < 1  =>  (I - A)^{-1} = I + A + A^2 + ... >= 0
d = np.array([100.0, 50.0, 80.0])                            # external demand
x = np.linalg.solve(np.eye(3) - A, d)                        # production needed
series = sum(np.linalg.matrix_power(A, k) for k in range(200)) @ d
print(np.round(x, 2), np.allclose(x, series))                # [255.74 240.98 322.95] True

Key takeaways

• "Positive matrix" means positive entries, not positive definite; the two theories are unrelated.
• For $A \ge 0$, $(A^k)_{ij} > 0$ exactly when there is a walk of length $k$ from $i$ to $j$; the zero pattern is a directed graph.
• Irreducible means strongly connected; primitive means some power is positive, which for irreducible matrices is the same as being aperiodic.
• $\rho(A)$ lies between the smallest and largest row sums and grows monotonically with the entries.
• If $\rho(A) < 1$, then $(I - A)^{-1} = \sum_k A^k$ is nonnegative.

🧪 Practice

1. Draw the graph of $A = \begin{bmatrix} 0 & 1 & 0 \ 0 & 0 & 1 \ 1 & 1 & 0 \end{bmatrix}$. Is $A$ irreducible? What is its period? Find the smallest $k$ with $A^k > 0$.
2. Show that if $A \ge 0$ is irreducible and has at least one positive diagonal entry, then $A$ is primitive.
3. Prove the upper row-sum bound $\rho(A) \le \max_i\sum_j a_{ij}$ using the induced $\infty$-norm (Section 11.3).
4. Interview: Given a directed graph with a million nodes stored as a sparse adjacency matrix, how would you decide whether it is strongly connected? Why is computing $(I + A)^{n-1}$ a bad idea? Hint: two breadth-first searches, on $A$ and on $A^T$, cost $O(n + m)$ and never fill in the matrix; Tarjan's algorithm also finds all strongly connected components.

Perron-Frobenius Theorem

Why it matters. Section 8.4 computed Markov steady states, PageRank scores, and population growth rates as dominant eigenvectors, and assumed that such an eigenvector exists, is unique, has nonnegative entries (so it can be read as probabilities or population shares), and is the limit of power iteration. The Perron-Frobenius theorem is the guarantee behind all of those assumptions, and it states exactly when they can fail.

Analogy. Picture repeatedly redistributing water through a network of pipes that only ever carry water forward (all entries nonnegative). If every tank eventually feeds every other tank, the proportions of water in the tanks settle to a single pattern, whatever the starting amounts. The total grows or shrinks by a fixed factor per step: that factor is the Perron root, and the pattern is the Perron vector.

The theorem for positive matrices (Perron, 1907). Let $A > 0$ be $n \times n$. Then:

1. $r = \rho(A)$ is a positive real eigenvalue of $A$ (the Perron root).
2. $r$ is algebraically simple: a simple root of the characteristic polynomial.
3. There is an eigenvector $\mathbf{p} > 0$ with $A\mathbf{p} = r\mathbf{p}$ (the Perron vector), and every nonnegative eigenvector of $A$ is a multiple of $\mathbf{p}$.
4. Every other eigenvalue satisfies $\lvert\lambda\rvert < r$ (strict dominance).

The theorem for irreducible matrices (Frobenius, 1912). Let $A \ge 0$ be irreducible with period $h$. Then items 1-3 still hold, but item 4 is replaced by:

4'. There are exactly $h$ eigenvalues of modulus $r$, namely $r, r\omega, r\omega^2, \dots, r\omega^{h-1}$ with $\omega = e^{2\pi i/h}$, and the whole spectrum is invariant under rotation by $2\pi/h$.

So strict dominance holds exactly when $A$ is primitive ($h = 1$). For a general (reducible) $A \ge 0$, only a weak form survives: $\rho(A)$ is an eigenvalue with some nonnegative eigenvector, which need not be unique or positive.

Class	$\rho(A)$ an eigenvalue?	Simple?	Positive eigenvector?	Strictly dominant?	Power iteration converges?
positive / primitive	yes	yes	yes, unique	yes	yes
irreducible, period $h > 1$	yes	yes	yes, unique	no: $h$ eigenvalues on the circle	no (oscillates with period $h$)
reducible	yes	not necessarily	nonnegative only	not necessarily	depends on the start

Proof idea (positive case). Define, for $\mathbf{x} \ge 0$, $\mathbf{x} \ne \mathbf{0}$, the number $f(\mathbf{x}) = \min_{i:,x_i > 0}(A\mathbf{x})_i/x_i$: the largest $t$ with $A\mathbf{x} \ge t\mathbf{x}$. Maximize $f$ over the compact set of nonnegative vectors with entries summing to 1. At a maximizer $\mathbf{p}$, if $A\mathbf{p} - r\mathbf{p}$ were nonzero, it would be nonnegative, so $A(A\mathbf{p} - r\mathbf{p}) > 0$ because $A > 0$. Then $A\mathbf{p}$ would give a larger value of $f$, a contradiction. So $A\mathbf{p} = r\mathbf{p}$, and $\mathbf{p} = A\mathbf{p}/r > 0$. Comparing $\lvert A\mathbf{x} \rvert \le A\lvert\mathbf{x}\rvert$ for any eigenvector $\mathbf{x}$ shows $\lvert\lambda\rvert \le r$, with equality only when the entries of $\mathbf{x}$ share a common phase. The max-min construction gives the useful Collatz-Wielandt bounds: for any $\mathbf{x} > 0$,

$$\min_i \frac{(A\mathbf{x})_i}{x_i} ;\le; \rho(A) ;\le; \max_i \frac{(A\mathbf{x})_i}{x_i}.$$

Consequences.

• Power iteration (Section 12.3) on a primitive matrix, started from any positive vector, converges to the Perron vector at rate $\lvert\lambda_2\rvert/r$.
• Markov chains: a column-stochastic matrix has $r = 1$ (column sums 1), and the left eigenvector $\mathbf{1}^T$ is positive. If the chain is irreducible, the stationary distribution is unique and positive. If it is also aperiodic, every starting distribution converges to it.
• PageRank: the raw link matrix of the web is reducible (dangling pages, closed subgraphs). The Google matrix $G = \alpha S + (1 - \alpha)\tfrac{1}{n}\mathbf{1}\mathbf{1}^T$ is positive, so the ranking exists, is unique, and power iteration converges with $\lvert\lambda_2\rvert \le \alpha$.
• Population models: for a Leslie matrix, $r$ is the long-run growth factor per period, and $\mathbf{p}$ is the stable age distribution.

import numpy as np

# Perron vector of a positive matrix by power iteration (Section 12.3)
A = np.array([[2.0, 1.0, 1.0], [1.0, 3.0, 1.0], [1.0, 1.0, 1.0]])
x = np.ones(3)
for _ in range(50):
    x = A @ x
    x /= x.sum()                              # normalize to sum 1 (it stays positive)
lam = (A @ x)[0] / x[0]
ev = np.linalg.eigvals(A)
print(np.round(x, 4), round(lam, 4))          # [0.3111 0.4516 0.2373] 4.2143: positive vector, rho(A)
print(np.round(sorted(abs(ev), reverse=True), 4))   # [4.2143 1.4608 0.3249]: rho is strictly dominant

# Collatz-Wielandt: min_i (Ax)_i / x_i <= rho <= max_i (Ax)_i / x_i for any positive x
y = np.array([1.0, 2.0, 1.0])
r = (A @ y) / y
print(r.min(), r.max())                       # 4.0 5.0: brackets rho = 4.2143

# Irreducible but periodic: h = 3 eigenvalues on the circle |z| = rho, evenly spaced
C = np.array([[0, 2, 0], [0, 0, 1], [4, 0, 0]], dtype=float)  # cycle 1 -> 2 -> 3 -> 1 with weights
ev = np.linalg.eigvals(C)
print(np.round(abs(ev), 4), np.round(np.angle(ev, deg=True) % 360, 1))
# [2. 2. 2.] [120. 240.   0.]     spectrum invariant under rotation by 360/3 degrees

# PageRank (Section 8.4): the damping term makes the Google matrix positive
links = {0: [1, 2], 1: [2], 2: [0], 3: [2]}   # page -> pages it links to; nobody links to 3
n, alpha = 4, 0.85
S = np.zeros((n, n))
for j, outs in links.items():
    S[outs, j] = 1 / len(outs)                # column-stochastic link matrix
G = alpha * S + (1 - alpha) / n               # every entry > 0
print((G > 0).all(), np.allclose(G.sum(axis=0), 1))           # True True
w, V = np.linalg.eig(G)
k = np.argmax(w.real)
pr = np.abs(V[:, k].real); pr /= pr.sum()
print(round(w[k].real, 6), np.round(pr, 4))   # 1.0 [0.3725 0.1958 0.3941 0.0375]
print(np.round(np.sort(abs(w))[::-1], 4))     # [1.    0.601 0.601 0.   ]: |lambda_2| <= 0.85 sets the convergence rate

# Leslie population model: the Perron root is the long-run growth rate per period
L = np.array([[0.0, 1.5, 1.2], [0.6, 0.0, 0.0], [0.0, 0.8, 0.0]])  # fertilities; survival rates
w, V = np.linalg.eig(L)
k = np.argmax(abs(w))
stable = np.abs(V[:, k].real); stable /= stable.sum()
print(round(w[k].real, 4), np.round(stable, 3))   # 1.1785 [0.539 0.275 0.186]: grows about 18% per period

Key takeaways

• For $A > 0$, $\rho(A)$ is a simple, strictly dominant eigenvalue with a positive eigenvector, the only nonnegative eigenvector up to scale.
• Irreducibility keeps simplicity and positivity, but a period $h > 1$ puts $h$ eigenvalues on the circle $\lvert z \rvert = \rho$, spaced by $2\pi/h$.
• Strict dominance, and therefore convergence of power iteration, requires primitivity.
• Collatz-Wielandt: any positive test vector brackets $\rho(A)$ between its smallest and largest ratios $(A\mathbf{x})_i/x_i$.
• PageRank's damping and the self-loops of lazy random walks are there to make the matrix primitive.

🧪 Practice

1. For $A = \begin{bmatrix} 1 & 2 \ 3 & 2 \end{bmatrix}$, find the Perron root and a positive Perron vector by hand. Check that the other eigenvalue is smaller in modulus.
2. Show that the permutation matrix of an $n$-cycle is irreducible with period $n$, and that its eigenvalues are exactly the $n$-th roots of unity.
3. Use the Collatz-Wielandt bounds with $\mathbf{x} = (1, 1, 1)$ and then with $\mathbf{x} = (1, 2, 1)$ to bracket $\rho(A)$ for the $3 \times 3$ matrix in the code. Which test vector gives the tighter bracket, and why?
4. Interview: Your PageRank power iteration on a crawl converges very slowly when you set $\alpha = 0.99$. Explain why, and what trade-off the choice of $\alpha$ represents. Hint: $\lvert\lambda_2(G)\rvert \le \alpha$, so the error shrinks roughly like $\alpha^k$; a small $\alpha$ converges fast but pulls the ranking toward uniform.

Adjacency and Incidence Matrices

Why it matters. A graph is a combinatorial object: vertices and edges. Writing it as a matrix makes the whole of linear algebra available. Matrix powers count walks, eigenvalues detect bipartiteness and bound degrees, and the four fundamental subspaces of the incidence matrix (Section 6.6) are exactly the conservation laws of electrical networks (Section 3.5). This topic introduces the two basic encodings.

Adjacency matrix. For a graph on vertices $1, \dots, n$, the adjacency matrix has $A_{ij} = 1$ if there is an edge from $i$ to $j$ and $0$ otherwise (or the edge weight, for weighted graphs). For an undirected graph, $A$ is symmetric, so it has real eigenvalues and an orthonormal eigenbasis (Section 10.1), and the degree of vertex $i$ is the row sum $d_i = \sum_jA_{ij}$.

Walk counting. As in the first topic of this subchapter, $(A^k)_{ij}$ is the number of walks of length $k$ from $i$ to $j$. This yields simple formulas for an undirected simple graph:

• $\operatorname{tr}(A) = 0$ (no loops), $\operatorname{tr}(A^2) = 2m$ (each edge walked there and back from both ends).
• $\operatorname{tr}(A^3) = 6 \times (\text{number of triangles})$: each triangle has 3 starting vertices and 2 directions.
• Since $\operatorname{tr}(A^k) = \sum_i\lambda_i^k$, these counts are determined by the spectrum.

What the adjacency spectrum tells you.

Spectral fact	Graph meaning
$\bar{d} \le \lambda_{\max} \le d_{\max}$	$\lambda_{\max}$ lies between the average and maximum degree; equals $d$ for a $d$-regular graph
spectrum symmetric about 0	bipartite (for connected graphs, iff $-\lambda_{\max}$ is an eigenvalue)
Perron vector of $A$	eigenvector centrality: important vertices are those linked to important vertices
number of distinct eigenvalues	greater than the diameter

Cospectral graphs exist (non-isomorphic graphs with the same eigenvalues), so the spectrum does not determine the graph completely.

Incidence matrix. Give each edge an arbitrary orientation. The oriented incidence matrix $B$ is $n \times m$, with one column per edge: for edge $e = (i \to j)$, $B_{ie} = -1$, $B_{je} = +1$, and the rest of the column is zero. Its transpose and itself act as discrete calculus:

• $B^T\mathbf{x}$ assigns to each edge the difference $x_j - x_i$ of a vertex potential $\mathbf{x}$: a discrete gradient.
• $B\mathbf{y}$ assigns to each vertex the net inflow of an edge flow $\mathbf{y}$: a discrete divergence (up to sign).

The four subspaces of $B$. For a graph with $n$ vertices, $m$ edges, and $c$ connected components:

Subspace	Description	Dimension	Physical law
$N(B^T)$	potentials constant on each component	$c$	potentials are defined only up to a constant
$C(B^T)$	edge vectors that are potential differences	$n - c$	Kirchhoff's voltage law: they sum to zero around every loop
$N(B)$	flows conserved at every vertex: the cycle space	$m - n + c$	Kirchhoff's current law
$C(B)$	net vertex inflows produced by some flow	$n - c$	total supply on each component is zero

So $\operatorname{rank}B = n - c$, and the number of independent loops in a connected graph is $m - n + 1$, which is the number of loop equations a circuit analysis needs. The product $BB^T$ does not depend on the chosen orientations: it equals $D - A$, the graph Laplacian of the next topic. The unoriented incidence matrix (all entries $+1$) satisfies $|B||B|^T = D + A$ instead.

import numpy as np

#   0 --- 1
#    \   /
#     \ /
#      2 --- 3 --- 4
edges = [(0, 1), (0, 2), (1, 2), (2, 3), (3, 4)]
n, m = 5, len(edges)

A = np.zeros((n, n), dtype=int)
for i, j in edges:
    A[i, j] = A[j, i] = 1                     # undirected: symmetric
print(A.sum(axis=1))                          # [2 2 3 2 1]: degrees

A2, A3 = A @ A, A @ A @ A
print(A2[0, 3], A3[0, 0])                     # 1 2: one walk 0-2-3; two closed 3-walks at 0 (both directions round the triangle)
print(np.trace(A3) // 6)                      # 1 triangle: each is counted 3 vertices x 2 directions

lam = np.linalg.eigvalsh(A)
print(np.round(lam, 4))                       # [-1.6751 -1.     -0.5392  1.      2.2143]: not symmetric about 0: not bipartite
print(2 * m / n, round(lam.max(), 4), A.sum(axis=1).max())   # 2.0 2.2143 3: average degree <= lambda_max <= max degree

# Oriented incidence matrix: column e = (i, j) has -1 at the tail i and +1 at the head j
Bm = np.zeros((n, m), dtype=int)
for e, (i, j) in enumerate(edges):
    Bm[i, e], Bm[j, e] = -1, 1
print(Bm)
# [[-1 -1  0  0  0]
#  [ 1  0 -1  0  0]
#  [ 0  1  1 -1  0]
#  [ 0  0  0  1 -1]
#  [ 0  0  0  0  1]]

x = np.array([0.0, 1.0, 3.0, 6.0, 10.0])      # a potential (value) at each vertex
print(Bm.T @ x)                               # [1. 3. 2. 3. 4.]: B^T is the gradient (difference along each edge)

print(np.linalg.matrix_rank(Bm))              # 4 = n - (number of components)
# Cycle space = null(B): flows that are conserved at every vertex (Kirchhoff's current law)
_, s, Vt = np.linalg.svd(Bm.astype(float))
cycles = Vt[np.sum(s > 1e-10):]               # basis of null(B)
print(len(cycles), np.round(cycles / cycles[0, 0], 4) + 0.0)   # 1 [[ 1. -1.  1.  0.  0.]]: the loop 0 -> 1 -> 2 -> 0

D = np.diag(A.sum(axis=1))
print(np.array_equal(Bm @ Bm.T, D - A))       # True: B B^T is the graph Laplacian (next topic)

Key takeaways

• The adjacency matrix $A$ is symmetric for undirected graphs; $(A^k)_{ij}$ counts walks of length $k$, and $\operatorname{tr}(A^3)/6$ counts triangles.
• $\lambda_{\max}(A)$ lies between the average and maximum degree; a spectrum symmetric about 0 signals a bipartite graph.
• The oriented incidence matrix $B$ is a discrete divergence, and $B^T$ is a discrete gradient.
• $\operatorname{rank}B = n - c$; the null space of $B$ is the cycle space, of dimension $m - n + c$.
• $BB^T = D - A$ regardless of the orientation chosen.

🧪 Practice

1. Write $A$ and $B$ for the 4-cycle $1 - 2 - 3 - 4 - 1$. Compute $A^2$ and interpret its diagonal and off-diagonal entries.
2. Show that the eigenvalues of the complete graph $K_n$ are $n - 1$ (once) and $-1$ ($n - 1$ times), using $A = J - I$ where $J$ is the all-ones matrix.
3. Prove that if a graph is bipartite, its adjacency spectrum is symmetric about 0. Hint: order the vertices by side, so that $A = \begin{bmatrix} 0 & C \ C^T & 0 \end{bmatrix}$, and flip the sign of half of an eigenvector.
4. Interview: How would you count the triangles in a large undirected graph using linear algebra, and what is the cost of the naive approach? Hint: $\operatorname{tr}(A^3)/6$; computing $A^3$ densely is $O(n^3)$, but $\sum_{ij}(A^2 \circ A)_{ij}/6$ with sparse products is far cheaper.

Graph Laplacian

Why it matters. The adjacency matrix records which vertices are connected. The Laplacian measures how much a function on the vertices varies across those connections. That turns graph questions into eigenvalue questions: how many pieces the graph has, how well connected it is, how heat or information spreads, and where to cut it (next topic). It is the discrete counterpart of the Laplace operator $\nabla^2$ of physics, which is where its name comes from.

Analogy. Imagine each edge is a spring connecting two beads, and $x_i$ is the vertical displacement of bead $i$. The potential energy stored in the springs is $\tfrac12\sum_{\text{edges}}(x_i - x_j)^2$. The Laplacian is the matrix of that energy.

Definition. For an undirected graph with adjacency matrix $A$ and degree matrix $D = \operatorname{diag}(d_1, \dots, d_n)$,

$$L = D - A = BB^T.$$

The key identity, which follows from $\mathbf{x}^TL\mathbf{x} = \mathbf{x}^TBB^T\mathbf{x} = \lVert B^T\mathbf{x} \rVert^2$, is

$$\boxed{\mathbf{x}^TL\mathbf{x} = \sum_{{i,j} \in E}(x_i - x_j)^2.}$$

This sum is small when $\mathbf{x}$ is "smooth" on the graph, that is, when neighboring vertices have similar values. For weighted graphs, $A_{ij} = w_{ij}$ and the sum becomes $\sum w_{ij}(x_i - x_j)^2$.

Spectral properties. Order the eigenvalues as $0 = \lambda_1 \le \lambda_2 \le \cdots \le \lambda_n$.

• $L$ is symmetric positive semidefinite (it is $BB^T$), so all $\lambda_i \ge 0$.
• $L\mathbf{1} = \mathbf{0}$ because each row sums to $d_i - d_i = 0$, so $\lambda_1 = 0$.
• The multiplicity of eigenvalue 0 equals the number of connected components. If $\mathbf{x}^TL\mathbf{x} = 0$, then $x_i = x_j$ across every edge, so $\mathbf{x}$ is constant on each component. The null space is spanned by the component indicator vectors.
• $\lambda_2 > 0$ if and only if the graph is connected. Fiedler called $\lambda_2$ the algebraic connectivity: it is small when the graph has a bottleneck, and its eigenvector, the Fiedler vector, changes sign across that bottleneck.
• $\lambda_n \le 2d_{\max}$.

Matrix-tree theorem (Kirchhoff). The number of spanning trees of a connected graph equals any cofactor of $L$ (delete row $i$ and column $i$, take the determinant), and also

$$\tau(G) = \frac{1}{n}\lambda_2\lambda_3\cdots\lambda_n.$$

For the complete graph $K_n$, $L = nI - J$ has eigenvalues $0$ and $n$ (with multiplicity $n - 1$), which gives Cayley's formula $\tau(K_n) = n^{n-2}$.

Diffusion. Heat flowing along edges from hotter to cooler neighbors obeys $\dot{\mathbf{x}} = -L\mathbf{x}$, the graph heat equation. Expanding in the eigenbasis (Section 8.4) gives $\mathbf{x}(t) = \sum_k e^{-\lambda_kt}(\mathbf{u}_k^T\mathbf{x}_0)\mathbf{u}_k$. Every mode decays except the constant one, so the temperatures approach their average, and the slowest decay rate is $\lambda_2$. Consensus algorithms in distributed systems and sensor networks are built on exactly this equation.

Path graph and the second derivative. For a path $1 - 2 - \cdots - n$, row $i$ of $L$ is $(-1, 2, -1)$, the second-difference stencil $-(x_{i-1} - 2x_i + x_{i+1}) \approx -h^2x''$. Its eigenvalues are $2 - 2\cos(k\pi/n)$ for $k = 0, \dots, n - 1$, with cosine eigenvectors: the graph's own Fourier modes.

Normalized Laplacians. When degrees vary widely, high-degree vertices dominate $L$. Two rescalings are common:

Matrix	Definition	Symmetric?	Eigenvalues	Use
combinatorial $L$	$D - A$	yes	$[0, 2d_{\max}]$	regular-ish graphs, physics, RatioCut
symmetric normalized $L_{\text{sym}}$	$I - D^{-1/2}AD^{-1/2}$	yes	$[0, 2]$	spectral clustering, Normalized Cut
random walk $L_{\text{rw}}$	$I - D^{-1}A$	no (similar to $L_{\text{sym}}$)	$[0, 2]$	random walks: $D^{-1}A$ is the transition matrix

import numpy as np

def laplacian(n, edges):
    A = np.zeros((n, n))
    for i, j in edges:
        A[i, j] = A[j, i] = 1
    return np.diag(A.sum(axis=1)) - A, A

# Two triangles {0,1,2} and {3,4,5} joined by the single "bridge" edge 2-3
edges = [(0, 1), (0, 2), (1, 2), (3, 4), (3, 5), (4, 5), (2, 3)]
L, A = laplacian(6, edges)
x = np.array([1.0, 2.0, 0.0, -1.0, 3.0, 1.0])
print(x @ L @ x, sum((x[i] - x[j]) ** 2 for i, j in edges))   # 31.0 31.0: x^T L x = sum over edges of (x_i - x_j)^2

lam, U = np.linalg.eigh(L)
print(np.round(lam, 4) + 0.0)                 # [0.     0.4384 3.     3.     3.     4.5616]
print(np.round(np.abs(U[:, 0]) * np.sqrt(6), 4))   # [1. 1. 1. 1. 1. 1.]: the constant vector spans the null space
print(np.round(U[:, 1], 3))                   # [-0.465 -0.465 -0.261  0.261  0.465  0.465] (up to sign): Fiedler vector,
                                              # one sign per triangle, smallest magnitude at the bridge

# Removing the bridge disconnects the graph: the eigenvalue 0 becomes double
L2, _ = laplacian(6, edges[:-1])
print(np.round(np.linalg.eigvalsh(L2), 4) + 0.0)   # [0. 0. 3. 3. 3. 3.]: two components

# Matrix-tree theorem: number of spanning trees = any cofactor = (product of nonzero eigenvalues) / n
LK4, _ = laplacian(4, [(i, j) for i in range(4) for j in range(i + 1, 4)])
print(round(np.linalg.det(LK4[1:, 1:])), round(np.prod(np.linalg.eigvalsh(LK4)[1:]) / 4))   # 16 16 = 4^(4-2) (Cayley)
print(round(np.linalg.det(L[1:, 1:])))        # 9: one spanning tree per choice of an edge dropped from each triangle

# Heat diffusion dx/dt = -L x: every component converges to the average, at rate lambda_2
x0 = np.array([6.0, 0.0, 0.0, 0.0, 0.0, 0.0])
for t in (0.5, 2.0, 10.0):
    xt = U @ (np.exp(-lam * t) * (U.T @ x0))  # solve via the eigenbasis (Section 8.2)
    print(t, np.round(xt, 3))
# 0.5 [2.842 1.504 1.287 0.267 0.05  0.05 ]
# 2.0 [1.548 1.533 1.3   0.695 0.462 0.462]
# 10.0 [1.016 1.016 1.009 0.991 0.984 0.984]    -> average 1.0; the slowest mode exp(-0.4384 t) is the Fiedler vector

# Normalized Laplacian I - D^{-1/2} A D^{-1/2}: eigenvalues in [0, 2]
d = A.sum(axis=1)
Lsym = np.eye(6) - A / np.sqrt(np.outer(d, d))
print(np.round(np.linalg.eigvalsh(Lsym), 4) + 0.0)   # [0.     0.2047 1.1667 1.5    1.5    1.6287]

Key takeaways

• $L = D - A = BB^T$, and $\mathbf{x}^TL\mathbf{x} = \sum_{\text{edges}}(x_i - x_j)^2$ measures how much $\mathbf{x}$ varies across edges.
• $L$ is positive semidefinite, and the multiplicity of eigenvalue 0 is the number of connected components.
• $\lambda_2$ (algebraic connectivity) is small when the graph has a bottleneck; the Fiedler vector changes sign across it.
• Matrix-tree theorem: the number of spanning trees is any cofactor of $L$, or $\frac1n\prod_{k \ge 2}\lambda_k$.
• The heat equation $\dot{\mathbf{x}} = -L\mathbf{x}$ drives every vertex to the average, at rate $\lambda_2$.

🧪 Practice

1. Compute $L$ and its eigenvalues for the star graph with one center and three leaves. How many spanning trees does it have? Check with the matrix-tree theorem.
2. Prove that the eigenvalue 0 of $L$ has multiplicity equal to the number of connected components, using $\mathbf{x}^TL\mathbf{x} = \sum(x_i - x_j)^2$.
3. Show that $L_{\text{rw}} = D^{-1/2}L_{\text{sym}}D^{1/2}$, so the two normalized Laplacians have the same eigenvalues. How are their eigenvectors related?
4. Interview: You are designing a peer-to-peer network where nodes average their values with neighbors to reach consensus. Which single number predicts how fast consensus is reached, and how would you add a few links to speed it up? Hint: the convergence rate is governed by $\lambda_2(L)$; add edges between vertices whose Fiedler-vector entries are far apart (opposite signs).

Spectral Clustering

Why it matters. k-means clusters points around centroids, so it only finds round, convex groups. Real clusters can be rings, spirals, or communities in a social network with no coordinates at all. Spectral clustering builds a similarity graph, looks for a cut that separates it into well-connected pieces with few edges between them, and solves a relaxed version of that cut problem with the eigenvectors of the Laplacian. It is one of the clearest examples of an NP-hard combinatorial problem made tractable by turning it into an eigenvalue problem.

Analogy. Shake the spring network from the previous topic and watch the slowest vibration. Tightly connected groups of beads move together, and the weak springs between groups stretch. The slowest modes (smallest Laplacian eigenvalues) reveal the groups.

From cuts to eigenvectors. To split the vertices into $S$ and $\bar{S}$, encode the split as $x_i = +1$ on $S$ and $-1$ on $\bar{S}$. Then

$$\mathbf{x}^TL\mathbf{x} = \sum_{\text{edges}}(x_i - x_j)^2 = 4\cdot\operatorname{cut}(S, \bar{S}),$$

where $\operatorname{cut}$ counts (or sums the weights of) edges crossing between the two sides. Minimizing the cut alone is useless, because it prefers to chop off a single vertex. Balanced objectives fix this:

$$\operatorname{RatioCut} = \operatorname{cut}(S, \bar{S})\Big(\frac{1}{\lvert S \rvert} + \frac{1}{\lvert \bar{S} \rvert}\Big), \qquad \operatorname{Ncut} = \operatorname{cut}(S, \bar{S})\Big(\frac{1}{\operatorname{vol}S} + \frac{1}{\operatorname{vol}\bar{S}}\Big),$$

where $\operatorname{vol}S = \sum_{i \in S}d_i$. Minimizing either over discrete labelings is NP-hard. The relaxation lets $\mathbf{x}$ take any real values, subject to $\mathbf{x} \perp \mathbf{1}$ (balance) and $\lVert\mathbf{x}\rVert = \sqrt{n}$. Then, by the Courant-Fischer theorem (Section 10.3),

$$\min_{\mathbf{x} \perp \mathbf{1},,\mathbf{x} \ne \mathbf{0}}\frac{\mathbf{x}^TL\mathbf{x}}{\mathbf{x}^T\mathbf{x}} = \lambda_2,$$

attained at the Fiedler vector. Rounding by sign, $S = {i : u_{2,i} > 0}$, recovers a partition. RatioCut relaxes to $L$, and Ncut relaxes to the generalized problem $L\mathbf{u} = \lambda D\mathbf{u}$, which is $L_{\text{rw}}$ or $L_{\text{sym}}$. For $k$ clusters, the first $k$ eigenvectors play the same role.

The algorithm (Ng-Jordan-Weiss version).

input: points x_1..x_N, number of clusters k, scale sigma

1. similarity graph   W_ij = exp(-||x_i - x_j||^2 / (2 sigma^2))
                      (often kept only for k-nearest neighbors, then symmetrized)
2. Laplacian          L_sym = I - D^{-1/2} W D^{-1/2}
3. embedding          U = [u_1 ... u_k]: eigenvectors of the k smallest eigenvalues
                      normalize each ROW of U to unit length -> point y_i in R^k
4. cluster            run k-means on y_1..y_N; point i gets the label of y_i

In the ideal case of $k$ disconnected components, the bottom eigenvectors are (scaled) component indicators, so the rows $\mathbf{y}_i$ collapse to $k$ distinct points and k-means cannot fail. Small perturbations move the points only slightly, which is why the method is robust (Davis-Kahan theorem).

Choosing $k$: the eigengap heuristic. If the graph has $k$ clear clusters, $\lambda_1, \dots, \lambda_k$ are close to 0 and $\lambda_{k+1}$ is noticeably larger. Pick the $k$ that maximizes $\lambda_{k+1} - \lambda_k$.

Practical issues.

Issue	Effect	Remedy
scale $\sigma$ too small	graph falls apart into many tiny components	use a kNN graph, or a per-point scale $\sigma_i$ (distance to the 7th neighbor)
scale $\sigma$ too large	everything is similar to everything; clusters merge	tune on the eigengap
cost of a dense eigendecomposition	$O(N^3)$ time, $O(N^2)$ memory	sparse kNN graph plus Lanczos (eigsh, Section 12.3) for $k$ eigenvectors
unbalanced degrees	$L$ favors cutting off low-degree vertices	use $L_{\text{sym}}$ or $L_{\text{rw}}$ (Ncut)

import numpy as np

rng = np.random.default_rng(0)

# Two concentric rings: no straight line (or k-means centroid pair) separates them
n = 200
r = np.r_[np.full(n, 1.0), np.full(n, 3.0)] + 0.1 * rng.standard_normal(2 * n)
t = rng.uniform(0, 2 * np.pi, 2 * n)
X = np.c_[r * np.cos(t), r * np.sin(t)]
truth = np.r_[np.zeros(n, int), np.ones(n, int)]

def kmeans(X, k, iters=100):
    C = [X[0]]                                # farthest-point initialization (deterministic)
    for _ in range(k - 1):
        C.append(X[np.argmax(np.min([((X - c) ** 2).sum(1) for c in C], axis=0))])
    C = np.array(C)
    for _ in range(iters):
        labels = np.argmin(((X[:, None, :] - C[None]) ** 2).sum(-1), axis=1)
        C = np.array([X[labels == j].mean(axis=0) for j in range(k)])
    return labels

def accuracy(labels):
    a = (labels == truth).mean()
    return max(a, 1 - a)                      # cluster names are arbitrary

print(accuracy(kmeans(X, 2)))                 # 0.54: k-means cuts both rings in half (chance level)

# 1. Similarity graph: Gaussian kernel restricted to each point's 10 nearest neighbors
D2 = ((X[:, None, :] - X[None]) ** 2).sum(-1)
sigma = 0.5
W = np.exp(-D2 / (2 * sigma**2))
knn = np.argsort(D2, axis=1)[:, 1:11]
mask = np.zeros_like(W, dtype=bool)
mask[np.arange(len(X))[:, None], knn] = True
W = np.where(mask | mask.T, W, 0.0)           # symmetrize the kNN graph

# 2. Normalized Laplacian L_sym = I - D^{-1/2} W D^{-1/2}
d = W.sum(axis=1)
Lsym = np.eye(len(X)) - W / np.sqrt(np.outer(d, d))

# 3. Bottom k eigenvectors, rows normalized to unit length (Ng-Jordan-Weiss)
lam, U = np.linalg.eigh(Lsym)
print(np.round(lam[:4], 4) + 0.0)             # [0.     0.     0.0037 0.005 ]: eigenvalue 0 twice -> k = 2
k = 2
Y = U[:, :k] / np.linalg.norm(U[:, :k], axis=1, keepdims=True)

# 4. k-means in the embedded space
print(accuracy(kmeans(Y, k)))                 # 1.0

Key takeaways

• For a $\pm 1$ labeling, $\mathbf{x}^TL\mathbf{x} = 4\cdot\operatorname{cut}$; balanced cut objectives (RatioCut, Ncut) are NP-hard to minimize exactly.
• Relaxing to real vectors turns the problem into a Rayleigh quotient minimization, solved by the bottom eigenvectors of a Laplacian.
• The algorithm: similarity graph, normalized Laplacian, bottom-$k$ eigenvectors as new coordinates, then k-means.
• It finds non-convex clusters that k-means on raw coordinates cannot.
• Choose $k$ with the eigengap; choose $\sigma$ or the neighbor count carefully, and use sparse eigensolvers at scale.

🧪 Practice

1. For the two-triangle graph in the Laplacian topic, partition the vertices by the sign of the Fiedler vector and compute the RatioCut of that partition.
2. Rerun the code with sigma = 0.05 and sigma = 5, and with 3 or 50 neighbors instead of 10. Report what happens to the bottom eigenvalues and to the accuracy, and explain each change.
3. Show that $\mathbf{x}^TL\mathbf{x} = 4\cdot\operatorname{cut}(S, \bar{S})$ for $x_i = \pm 1$, and that $\sum_ix_i = 0$ forces $\lvert S \rvert = \lvert \bar{S} \rvert$.
4. Interview: Spectral clustering works perfectly on 2,000 points but runs out of memory on 2 million. Which steps are the bottleneck, and how would you scale it? Hint: the dense $N \times N$ affinity matrix and the full eigendecomposition; use a sparse kNN graph (approximate nearest neighbors), Lanczos for only $k$ eigenvectors, or Nyström landmark approximations.

<a id="133-linear-algebra-in-applied-domains"></a>

13.3 Linear Algebra in Applied Domains

This subchapter visits six fields where linear algebra is the main computational tool, not a supporting one. Each topic shows which earlier concept carries the application and what extra structure the domain adds.

Computer Graphics and Homogeneous Coordinates

Why it matters. Every frame of a video game or CAD view moves millions of vertices through a chain of transformations: place the model in the world, view it from a camera, project it onto the screen. GPUs are built to do one thing very fast, multiply $4 \times 4$ matrices by 4-vectors. The difficulty is that translation, the most basic motion, is not linear in $\mathbb{R}^3$, since it moves the origin. Homogeneous coordinates (introduced in 2D in Section 2.1) fix that by adding one coordinate, which turns translations, and even perspective projection, into matrix multiplications.

Analogy. Translation is like trying to slide a sheet of paper while only being allowed to rotate and stretch it around a pin at the origin. Lift the paper into 3D, stand it on the plane $w = 1$, and a shear of 3D space (which is linear and keeps the origin fixed) slides the paper sideways. Homogeneous coordinates use the same trick one dimension up.

Points and directions. Represent a point $(x, y, z)$ as $(x, y, z, 1)$ and a direction (displacement, velocity, normal) as $(x, y, z, 0)$. More generally, $(x, y, z, w)$ with $w \ne 0$ represents the point $(x/w, y/w, z/w)$, so all nonzero multiples of a 4-vector are the same point. This is the projective space $\mathbb{RP}^3$.

Affine transformations as $4 \times 4$ matrices. A linear part $M$ (rotation, scaling, shear) followed by a translation $\mathbf{t}$ becomes

$$\begin{bmatrix} M & \mathbf{t} \ \mathbf{0}^T & 1 \end{bmatrix}\begin{bmatrix} \mathbf{x} \ 1 \end{bmatrix} = \begin{bmatrix} M\mathbf{x} + \mathbf{t} \ 1 \end{bmatrix}, \qquad \begin{bmatrix} M & \mathbf{t} \ \mathbf{0}^T & 1 \end{bmatrix}\begin{bmatrix} \mathbf{d} \ 0 \end{bmatrix} = \begin{bmatrix} M\mathbf{d} \ 0 \end{bmatrix}.$$

Directions are automatically unaffected by translation. Composition is matrix multiplication, so an entire chain collapses into one matrix that is applied to every vertex.

Transform	Matrix (block form)	Inverse
translation by $\mathbf{t}$	$\begin{bmatrix} I & \mathbf{t} \ \mathbf{0}^T & 1 \end{bmatrix}$	translation by $-\mathbf{t}$
rotation $R$ (orthogonal, $\det R = 1$)	$\begin{bmatrix} R & \mathbf{0} \ \mathbf{0}^T & 1 \end{bmatrix}$	$R^T$ in place of $R$
scaling $S = \operatorname{diag}(s_x, s_y, s_z)$	$\begin{bmatrix} S & \mathbf{0} \ \mathbf{0}^T & 1 \end{bmatrix}$	$S^{-1}$
rigid motion (rotate, then translate)	$\begin{bmatrix} R & \mathbf{t} \ \mathbf{0}^T & 1 \end{bmatrix}$	$\begin{bmatrix} R^T & -R^T\mathbf{t} \ \mathbf{0}^T & 1 \end{bmatrix}$

Order and pivots. Matrices act right to left, and they do not commute: "rotate then translate" ($TR$) differs from "translate then rotate" ($RT$). To rotate about a point $\mathbf{c}$ instead of the origin, conjugate: $T(\mathbf{c}),R,T(-\mathbf{c})$ moves $\mathbf{c}$ to the origin, rotates, and moves it back. This is a change of basis (Section 7.2) in the affine setting.

Normals need the inverse transpose. A surface normal $\mathbf{n}$ is defined by being perpendicular to every tangent vector $\mathbf{t}$. If tangents transform by $M$, normals must transform by some $N$ with $(M\mathbf{t})^T(N\mathbf{n}) = \mathbf{t}^T\mathbf{n}$ for all $\mathbf{t}$, which forces $N = (M^{-1})^T$. For rotations, $N = M$. For non-uniform scales and shears it differs, and using $M$ produces visibly wrong lighting. Normals are really covectors (Section 13.1), which is why they transform differently.

Perspective projection. Objects twice as far away should look half as large: a point $(x, y, z)$ seen by a camera looking down the $-z$ axis appears at $(fx/(-z), fy/(-z))$. Division is not linear, but homogeneous coordinates absorb it: a matrix puts $-z$ into the $w$ slot, and the hardware's perspective divide by $w$ does the rest. A real pipeline uses a matrix that also maps the visible depth range to $[-1, 1]$ for the depth buffer.

  model coords --[model M]--> world --[view V]--> camera --[projection P]--> clip
       (x, y, z, 1)                                                          (x', y', z', w')
                                                                                  |
                                                       divide by w' (perspective) v
  pixels <--[viewport: scale + shift]-- normalized device coordinates in [-1, 1]^3

  one matrix per vertex:  clip = P V M p

Rotations are often stored as unit quaternions, which interpolate smoothly and avoid the gimbal lock of Euler angles, then converted to the rotation block $R$ for rendering.

import numpy as np

def translate(tx, ty, tz):
    M = np.eye(4); M[:3, 3] = [tx, ty, tz]; return M

def scale(sx, sy, sz):
    return np.diag([sx, sy, sz, 1.0])

def rotate_z(deg):
    c, s = np.cos(np.radians(deg)), np.sin(np.radians(deg))
    M = np.eye(4); M[:2, :2] = [[c, -s], [s, c]]; return M

p = np.array([1.0, 0.0, 0.0, 1.0])            # a point: w = 1
d = np.array([1.0, 0.0, 0.0, 0.0])            # a direction: w = 0
T = translate(5, 0, 0)
print(T @ p, T @ d)                           # [6. 0. 0. 1.] [1. 0. 0. 0.]: directions ignore translation

# Order matters: matrices apply right to left
R = rotate_z(90)
print(np.round(T @ R @ p, 6), np.round(R @ T @ p, 6))   # [5. 1. 0. 1.] [0. 6. 0. 1.]

# Rotate a unit square by 90 degrees about its own center c = (0.5, 0.5), not the origin
square = np.array([[0, 0, 0, 1], [1, 0, 0, 1], [1, 1, 0, 1], [0, 1, 0, 1]], float).T   # columns are points
M = translate(0.5, 0.5, 0) @ R @ translate(-0.5, -0.5, 0)   # move c to origin, rotate, move back
print(np.round(M @ square, 6)[:2])            # [[1. 1. 0. 0.] [0. 1. 1. 0.]]: same square, corners cycled

# Normals transform by the inverse transpose of the linear part
S = scale(2, 1, 1)[:3, :3]                    # non-uniform scale
tangent = np.array([1.0, -1.0, 0.0])          # along the surface x + y = 1
normal = np.array([1.0, 1.0, 0.0])            # perpendicular to it
print((S @ tangent) @ (S @ normal))           # 3.0: naive transform breaks perpendicularity
print((S @ tangent) @ (np.linalg.inv(S).T @ normal))     # 0.0: inverse transpose keeps it

# Perspective projection: camera at the origin looking down -z, image plane at distance f
def perspective(f):
    P = np.zeros((4, 4))
    P[0, 0] = P[1, 1] = f
    P[2, 2] = 1.0                             # keep depth (a real pipeline maps it to [-1, 1])
    P[3, 2] = -1.0                            # w' = -z: the divide by w is the perspective
    return P

pts = np.array([[1, 1, -2, 1], [1, 1, -4, 1], [2, 2, -4, 1]], float).T
clip = perspective(1.0) @ pts
ndc = clip[:2] / clip[3]                      # homogeneous divide
print(ndc.T)                                  # [[0.5 0.5] [0.25 0.25] [0.5 0.5]]: twice as far looks half as big

Key takeaways

• Homogeneous coordinates append $w$: points have $w = 1$, directions $w = 0$, and $(x, y, z, w) \sim (x/w, y/w, z/w)$.
• Translations become linear maps in one higher dimension, so every affine transform, and any chain of them, is one $4 \times 4$ matrix.
• Order matters (right to left); rotation about a point $\mathbf{c}$ is $T(\mathbf{c})RT(-\mathbf{c})$.
• Normals transform by the inverse transpose of the linear part.
• Perspective is a matrix that writes depth into $w$, followed by the divide by $w$.

🧪 Practice

1. Write the $3 \times 3$ homogeneous matrix (2D) that scales by 2 about the point $(1, 1)$. Check that $(1, 1)$ is fixed and $(2, 1)$ maps to $(3, 1)$.
2. Derive the inverse of the rigid motion $\begin{bmatrix} R & \mathbf{t} \ \mathbf{0}^T & 1 \end{bmatrix}$ and verify it by multiplication.
3. Show that under the perspective matrix of the code, all points on the line through the origin in direction $(a, b, c)$ with $c < 0$ map to the same image point. What does this mean physically?
4. Interview: After applying a non-uniform scale to a model, its lighting looks wrong even though the geometry is correct. What is the bug, and what is the fix? Hint: the code transforms normals with the model matrix; normals must stay perpendicular to tangents, which requires $(M^{-1})^T$ (and renormalization).

Cryptography and Coding Theory

Why it matters. Digital communication must survive two adversaries: noise, which flips bits at random, and eavesdroppers, who try to read or alter messages. Coding theory fights noise by adding structured redundancy, and the best-understood codes are subspaces of $\mathbb{F}_q^n$, vector spaces over finite fields (Section 1.2). Cryptography has a more complicated relationship with linearity: linear ciphers are easy to break, yet the leading post-quantum public-key systems are built on linear algebra with deliberately added noise.

Analogy. A linear code is a dictionary of valid words chosen so that any two differ in many letters. A smudged word is corrected to the nearest valid one. Linearity means the dictionary is a subspace, so it can be described by a small matrix instead of a list of every word.

Linear codes. A linear $[n, k, d]$ code over $\mathbb{F}_2$ is a $k$-dimensional subspace $C \subseteq \mathbb{F}_2^n$:

• A generator matrix $G$ ($k \times n$) has rows forming a basis of $C$. Encoding a message $\mathbf{m} \in \mathbb{F}_2^k$ is $\mathbf{c} = \mathbf{m}G$.
• A parity-check matrix $H$ ($(n - k) \times n$) has $C$ as its null space: $\mathbf{c} \in C \iff H\mathbf{c} = \mathbf{0}$. If $G = [I_k \mid P]$ (systematic form), then $H = [P^T \mid I_{n-k}]$ works over $\mathbb{F}_2$, since $GH^T = P + P = 0$.
• The Hamming distance $d(\mathbf{x}, \mathbf{y})$ counts the positions where $\mathbf{x}$ and $\mathbf{y}$ differ. Because $C$ is a subspace, $d(\mathbf{x}, \mathbf{y}) = \operatorname{wt}(\mathbf{x} - \mathbf{y})$ with $\mathbf{x} - \mathbf{y} \in C$, so the minimum distance $d$ is the smallest weight of a nonzero codeword.
• A code with minimum distance $d$ detects up to $d - 1$ errors and corrects up to $\lfloor(d-1)/2\rfloor$.

Syndrome decoding. If $\mathbf{c}$ is sent and $\mathbf{r} = \mathbf{c} + \mathbf{e}$ is received, then

$$H\mathbf{r} = H\mathbf{c} + H\mathbf{e} = H\mathbf{e}.$$

The syndrome $H\mathbf{r}$ depends only on the error, not on the message. For the Hamming $[7, 4, 3]$ code, the columns of $H$ are the seven distinct nonzero vectors of $\mathbb{F}_2^3$. A single error in position $j$ gives syndrome equal to column $j$, which pinpoints the bit to flip. The decoder never needs a list of all 16 codewords.

Code	$[n, k, d]$	Rate $k/n$	Corrects	Used in
repetition (3x)	$[3, 1, 3]$	0.33	1 error	teaching, simple telemetry
Hamming	$[7, 4, 3]$	0.57	1 error	ECC memory (extended variants)
extended Golay	$[24, 12, 8]$	0.5	3 errors	Voyager probes
Reed-Solomon (over $\mathbb{F}_{256}$)	$[255, 223, 33]$	0.87	16 symbol errors	CDs, QR codes, deep-space links

Linear ciphers and why they fail. The Hill cipher encrypts blocks of letters (as numbers mod 26) by $\mathbf{c} = K\mathbf{p} \pmod{26}$ with an invertible key matrix $K$. Decryption is $K^{-1}\mathbf{c}$, which requires $\gcd(\det K, 26) = 1$ so that the determinant can be inverted mod 26. Its weakness is linearity itself. An attacker who knows a few plaintext-ciphertext pairs stacks them into matrices with $C = KP$ and solves $K = CP^{-1}$. This is Gaussian elimination mod 26. Modern block ciphers such as AES alternate linear mixing layers (matrix multiplication over $\mathbb{F}_{2^8}$) with nonlinear substitution boxes, precisely so that no system of linear equations describes them.

Linear algebra with noise: lattice cryptography. In the Learning With Errors (LWE) problem, you are given a random matrix $A \in \mathbb{Z}_q^{m \times n}$ and

$$\mathbf{b} = A\mathbf{s} + \mathbf{e} \pmod q,$$

where $\mathbf{s}$ is secret and $\mathbf{e}$ has small random entries. Without $\mathbf{e}$, elimination recovers $\mathbf{s}$ in polynomial time. With $\mathbf{e}$, the problem is believed to be hard even for quantum computers, because elimination amplifies the noise until it swamps the answer. The NIST post-quantum standard ML-KEM (Kyber) is built on a structured variant of LWE.

import numpy as np

# Hamming(7,4) code over F_2: generator G = [I | P], parity check H = [P^T | I]
P = np.array([[1, 1, 0], [1, 0, 1], [0, 1, 1], [1, 1, 1]])
G = np.hstack([np.eye(4, dtype=int), P])      # 4 x 7
H = np.hstack([P.T, np.eye(3, dtype=int)])    # 3 x 7
print(np.all((G @ H.T) % 2 == 0))             # True: every row of G (so every codeword) satisfies H c = 0

msg = np.array([1, 0, 1, 1])
c = msg @ G % 2                               # encode: 4 data bits -> 7 bits
print(c)                                      # [1 0 1 1 0 1 0]

r = c.copy(); r[2] ^= 1                       # the channel flips bit 2
s = H @ r % 2                                 # syndrome depends only on the error, not on c
print(s, H[:, 2])                             # [0 1 1] [0 1 1]: syndrome = column 2 of H
pos = next(j for j in range(7) if np.array_equal(H[:, j], s))
r[pos] ^= 1
print(pos, np.array_equal(r, c), r[:4])       # 2 True [1 0 1 1]: corrected, message recovered

# Minimum distance = smallest weight of a nonzero codeword (the code is a subspace)
words = [np.array([(k >> i) & 1 for i in range(4)]) @ G % 2 for k in range(1, 16)]
print(min(w.sum() for w in words))            # 3: detects 2 errors, corrects 1

# Hill cipher: encrypt pairs of letters with an invertible 2x2 key matrix mod 26
def to_nums(s): return np.array([ord(ch) - 65 for ch in s]).reshape(-1, 2).T   # columns = letter pairs
def to_text(M): return "".join(chr(int(v) + 65) for v in M.T.ravel())
def inv_mod26(K):
    det = int(round(np.linalg.det(K))) % 26
    adj = np.array([[K[1, 1], -K[0, 1]], [-K[1, 0], K[0, 0]]])
    return (pow(det, -1, 26) * adj) % 26      # needs gcd(det, 26) = 1

K = np.array([[3, 3], [2, 5]])                # det = 9, invertible mod 26
plain = "HELPMEPLEASE"
cipher = K @ to_nums(plain) % 26
print(to_text(cipher), to_text(inv_mod26(K) @ cipher % 26))   # HIATWSAHMIOE HELPMEPLEASE

# Known-plaintext attack: two known pairs give the key by solving a linear system mod 26
Pk, Ck = to_nums("HELP"), K @ to_nums("HELP") % 26
K_found = Ck @ inv_mod26(Pk) % 26             # K = C P^{-1} (mod 26), if P is invertible mod 26
print(K_found)                                # [[3 3] [2 5]]: key recovered

Key takeaways

• A linear code is a subspace of $\mathbb{F}_q^n$, encoded by a generator matrix $G$ and checked by a parity-check matrix $H$ with $C = N(H)$.
• Minimum distance equals the minimum nonzero codeword weight; distance $d$ corrects $\lfloor(d-1)/2\rfloor$ errors.
• The syndrome $H\mathbf{r} = H\mathbf{e}$ depends only on the error; for Hamming codes it directly names the flipped bit.
• Purely linear ciphers fall to known-plaintext attacks by solving a linear system modulo $q$; secure ciphers mix in nonlinearity.
• LWE adds small noise to a linear system, which turns easy elimination into a problem believed to be hard even for quantum computers.

🧪 Practice

1. Encode the message $(0, 1, 1, 0)$ with the Hamming code from the example, flip bit 5, and decode by computing the syndrome by hand.
2. Show that the Hamming $[7, 4, 3]$ code is perfect: the radius-1 balls around its 16 codewords contain $16 \times 8 = 128 = 2^7$ vectors, so they tile $\mathbb{F}_2^7$ exactly.
3. Which of $\begin{bmatrix} 2 & 3 \ 1 & 4 \end{bmatrix}$ and $\begin{bmatrix} 3 & 1 \ 4 & 2 \end{bmatrix}$ is a valid Hill key mod 26? For the valid one, compute $K^{-1} \bmod 26$.
4. Interview: Why does adding a small error vector make $\mathbf{b} = A\mathbf{s} + \mathbf{e} \pmod q$ hard to solve, when $\mathbf{b} = A\mathbf{s}$ is solved by Gaussian elimination? Hint: follow what row operations with multipliers of size up to $q$ do to the entries of $\mathbf{e}$.

Linear Programming and the Simplex Method

Why it matters. Airlines schedule crews, refineries blend crude oil, and logistics networks route shipments by solving linear programs with millions of variables. A linear program optimizes a linear objective subject to linear inequality constraints. The simplex method, one of the most important algorithms of the 20th century, solves it by moving along the edges of a polyhedron, and each move is one Gauss-Jordan pivot (Section 3.3).

Analogy. Picture the feasible region as a gem with flat faces. Gravity is the objective, and you want the lowest point. It is always at a corner, or along a whole edge or face that includes corners. The simplex method starts at one corner and repeatedly rolls along an edge to a neighboring corner that is lower, until no neighbor is lower.

Standard form. In inequality form,

$$\max; \mathbf{c}^T\mathbf{x} \quad \text{subject to} \quad A\mathbf{x} \le \mathbf{b},; \mathbf{x} \ge \mathbf{0}.$$

Each constraint $\mathbf{a}_i^T\mathbf{x} \le b_i$ is a half-space, so the feasible set is a convex polyhedron. Adding a slack variable $s_i \ge 0$ per row gives equalities $A\mathbf{x} + \mathbf{s} = \mathbf{b}$, a linear system with more unknowns than equations.

Geometry: optima occur at vertices. A linear function on a polyhedron has no interior maximum. Its level sets are parallel hyperplanes, and you can always slide toward a better one until you hit the boundary. If an optimum exists, some vertex attains it. Algebraically, a vertex is a basic feasible solution: choose $m$ columns of $[A \mid I]$ forming an invertible basis matrix $B$, set the other variables to zero, and solve $\mathbf{x}_B = B^{-1}\mathbf{b}$. The solution must also satisfy $\mathbf{x}_B \ge \mathbf{0}$.

     y
     |
   9 +
     |
   6 +=========@ (2,6)   <- optimum: 3x + 5y = 36
     |          \
     |   feasible\        constraints:
     |   region   \         x <= 4           (vertical line)
     |             @ (4,3)   2y <= 12         (y <= 6)
     |             |         3x + 2y <= 18    (slanted edge)
     @-------------@-----> x
   (0,0)         (4,0)

The simplex method. Each iteration does three things:

1. Pricing. Compute the reduced costs $\bar{\mathbf{c}}_N = \mathbf{c}_N - N^TB^{-T}\mathbf{c}_B$ of the nonbasic variables, which measure how fast the objective improves when each one is increased from zero. If none is positive (for maximization), the current vertex is optimal.
2. Ratio test. Increase one improving variable as far as possible while keeping $\mathbf{x}_B \ge \mathbf{0}$. The first basic variable to hit zero leaves the basis.
3. Pivot. Swap the two variables and update the tableau by one Gauss-Jordan elimination step on the pivot entry.

Production codes never form $B^{-1}$. They maintain an LU factorization of $B$ (Section 4.4) and update it after each swap, which is a rank-1 change.

Duality. Every LP has a dual. For the problem above it is

$$\min; \mathbf{b}^T\mathbf{y} \quad \text{subject to} \quad A^T\mathbf{y} \ge \mathbf{c},; \mathbf{y} \ge \mathbf{0}.$$

• Weak duality: $\mathbf{c}^T\mathbf{x} \le \mathbf{b}^T\mathbf{y}$ for any feasible pair, since $\mathbf{c}^T\mathbf{x} \le (A^T\mathbf{y})^T\mathbf{x} = \mathbf{y}^T(A\mathbf{x}) \le \mathbf{y}^T\mathbf{b}$.
• Strong duality: if either problem has an optimum, both do, and the optimal values are equal.
• Shadow prices: the optimal $y_i$ is the rate at which the optimal objective increases per unit increase of $b_i$. A constraint that is slack at the optimum has $y_i = 0$ (complementary slackness).

Complexity.

Method	Worst case	In practice	Notes
simplex (Dantzig, 1947)	exponential (Klee-Minty cubes)	typically a small multiple of $m$ pivots	excellent warm starts; exact vertex solutions
ellipsoid (Khachiyan, 1979)	polynomial	slow	first proof that LP is in P
interior point (Karmarkar, 1984)	polynomial	very fast on large sparse problems	each step solves a linear system (Newton's method)

import numpy as np
from scipy.optimize import linprog

# max 3x + 5y  s.t.  x <= 4,  2y <= 12,  3x + 2y <= 18,  x, y >= 0
c = np.array([3.0, 5.0])
A = np.array([[1.0, 0.0], [0.0, 2.0], [3.0, 2.0]])
b = np.array([4.0, 12.0, 18.0])

def simplex(c, A, b):
    """Tableau simplex for max c^T x, Ax <= b, x >= 0 with b >= 0 (the origin is feasible)."""
    m, n = A.shape
    T = np.zeros((m + 1, n + m + 1))
    T[:m, :n], T[:m, n:n + m], T[:m, -1] = A, np.eye(m), b    # slack variables make a basis
    T[-1, :n] = -c                            # objective row: z - c^T x = 0
    basis = list(range(n, n + m))
    while True:
        j = np.argmin(T[-1, :-1])             # entering variable: most negative reduced cost
        if T[-1, j] >= -1e-12:
            break                             # no improving direction: optimal
        col, ratios = T[:m, j], np.full(m, np.inf)
        pos = col > 1e-12
        ratios[pos] = T[:m, -1][pos] / col[pos]
        i = np.argmin(ratios)                 # leaving variable: ratio test keeps x >= 0
        if ratios[i] == np.inf:
            raise ValueError("unbounded")
        T[i] /= T[i, j]                       # pivot: one Gauss-Jordan step (Section 3.3)
        for r in range(m + 1):
            if r != i:
                T[r] -= T[r, j] * T[i]
        basis[i] = j
        x = np.zeros(n + m); x[basis] = T[:m, -1]
        print("vertex", x[:n], "objective", T[-1, -1])
    x = np.zeros(n + m); x[basis] = T[:m, -1]
    return x[:n], T[-1, -1], T[-1, n:n + m]   # last part: dual values (shadow prices)

x, z, y = simplex(c, A, b)
# vertex [0. 6.] objective 30.0                 walks (0,0) -> (0,6) -> (2,6)
# vertex [2. 6.] objective 36.0
print(x, z, y)                                # [2. 6.] 36.0 [0.  1.5 1. ]

# Library solver (it minimizes, so negate c); HiGHS also reports the duals
res = linprog(-c, A_ub=A, b_ub=b, bounds=[(0, None)] * 2, method="highs")
print(res.x, -res.fun, -res.ineqlin.marginals)   # [2. 6.] 36.0 [0.  1.5 1. ]

# Strong duality: min b^T y s.t. A^T y >= c, y >= 0 has the same optimal value
dual = linprog(b, A_ub=-A.T, b_ub=-c, bounds=[(0, None)] * 3, method="highs")
print(dual.x, dual.fun)                       # [0.  1.5 1. ] 36.0

Key takeaways

• An LP optimizes a linear objective over a polyhedron; if an optimum exists, one is attained at a vertex.
• Vertices are basic feasible solutions $\mathbf{x}_B = B^{-1}\mathbf{b} \ge \mathbf{0}$ for a choice of $m$ basic columns.
• Simplex moves between adjacent vertices: reduced costs choose the entering variable, the ratio test chooses the leaving one, and a Gauss-Jordan pivot updates the tableau.
• Every LP has a dual with the same optimal value; the dual variables are shadow prices of the constraints.
• Simplex is exponential in the worst case but fast in practice; interior-point methods are polynomial and dominate on very large problems.

🧪 Practice

1. Solve $\max;x + y$ subject to $x + 2y \le 4$, $3x + y \le 6$, $x, y \ge 0$ graphically by listing all vertices, then check with linprog.
2. In the example, increase the right-hand side of $2y \le 12$ to $13$ and resolve. Does the objective rise by the shadow price $y_2 = 1.5$? Why does this prediction eventually break down if you keep increasing it?
3. Write the dual of $\min;\mathbf{c}^T\mathbf{x}$ subject to $A\mathbf{x} = \mathbf{b}$, $\mathbf{x} \ge \mathbf{0}$, and prove weak duality for this pair.
4. Interview: Your LP solver reports "infeasible" for a scheduling model with thousands of constraints. How would you find which constraints conflict? Hint: Farkas' lemma gives a certificate of infeasibility (a dual ray); solvers can compute an irreducible infeasible subset (IIS), or you can add penalized slack variables to every constraint and minimize total violation.

Signal Processing and Discrete Fourier Transform

Why it matters. Audio, images, radio, and sensor data are sampled signals, vectors in $\mathbb{C}^n$. The discrete Fourier transform (DFT) is a change of basis from samples in time to amplitudes of pure frequencies. In the frequency basis, the most common operation in signal processing, filtering by convolution, becomes a diagonal matrix. The Fast Fourier Transform (FFT) computes the change of basis in $O(n\log n)$ instead of $O(n^2)$, and it is behind MP3, JPEG, Wi-Fi, MRI, and fast polynomial multiplication.

Analogy. A chord played on a piano is one sound wave, but your ear hears separate notes. The DFT does what the ear does: it re-expresses the same signal in the "note" basis, where each coefficient says how much of each frequency is present.

The DFT matrix. Let $\omega = e^{-2\pi i/n}$, a primitive $n$-th root of unity. The DFT of $\mathbf{x} \in \mathbb{C}^n$ is $\hat{\mathbf{x}} = F\mathbf{x}$ with

$$F_{jk} = \omega^{jk}, \qquad \hat{x}j = \sum{k=0}^{n-1}x_k,e^{-2\pi ijk/n}, \qquad j, k = 0, \dots, n - 1.$$

Column $k$ of $\bar{F}$ samples the complex sinusoid $e^{2\pi ijk/n}$, which oscillates $k$ times over the window.

Orthogonality. The columns of $F$ are orthogonal, each with squared norm $n$, because for $j \ne l$ the sum $\sum_k\omega^{(j - l)k}$ is a full geometric series of roots of unity and vanishes. So

$$F^F = nI, \qquad F^{-1} = \tfrac1nF^, \qquad \tfrac{1}{\sqrt{n}}F \text{ is unitary (Section 9.4)}.$$

Consequences: the inverse DFT is the same computation with conjugated $\omega$, and Parseval's identity $\lVert\mathbf{x}\rVert^2 = \tfrac1n\lVert\hat{\mathbf{x}}\rVert^2$ says the change of basis preserves energy. This is the finite version of Fourier series as orthogonal projection (Section 9.5).

Circulant matrices and convolution. A circulant matrix $C$ has each column equal to the previous one shifted down cyclically. It is completely determined by its first column $\mathbf{c}$, and $C\mathbf{x} = \mathbf{c} \circledast \mathbf{x}$ is the circular convolution $(\mathbf{c} \circledast \mathbf{x})t = \sum_mc_mx{(t - m) \bmod n}$. The key theorem is that every circulant matrix is diagonalized by the DFT:

$$C = F^{-1}\operatorname{diag}(F\mathbf{c}),F.$$

The Fourier vectors are eigenvectors of every circulant matrix (they are eigenvectors of the cyclic shift, and a circulant is a polynomial in the shift), and the eigenvalues are the DFT of the first column. Reading the formula right to left gives the convolution theorem:

$$\mathbf{c} \circledast \mathbf{x} = F^{-1}\big((F\mathbf{c}) \odot (F\mathbf{x})\big).$$

Convolution in time is pointwise multiplication in frequency. A filter is therefore a diagonal matrix in the Fourier basis, and its eigenvalues $F\mathbf{c}$ are the frequency response.

The FFT as a matrix factorization. For $n = 2^p$, splitting $\mathbf{x}$ into even- and odd-indexed entries gives

$$F_n = \begin{bmatrix} I & D \ I & -D \end{bmatrix}\begin{bmatrix} F_{n/2} & 0 \ 0 & F_{n/2} \end{bmatrix}P, \qquad D = \operatorname{diag}(1, \omega, \dots, \omega^{n/2 - 1}),$$

where $P$ is the even-odd permutation. Applying the factorization recursively writes $F_n$ as a product of $\log_2n$ sparse "butterfly" matrices, each costing $O(n)$ to apply, for $O(n\log n)$ in total. For $n = 10^6$ this is about 50,000 times fewer operations than the dense product.

Operation on signals	Matrix in time basis	Matrix in Fourier basis
cyclic shift by one sample	permutation $S$	$\operatorname{diag}(\omega^{-j})$ (phase ramp)
filtering (circular convolution)	circulant $C$	$\operatorname{diag}(F\mathbf{c})$
ideal low-pass filter	dense circulant	0/1 diagonal: projection onto low frequencies
energy $\lVert\mathbf{x}\rVert^2$	$\mathbf{x}^*\mathbf{x}$	$\tfrac1n\hat{\mathbf{x}}^*\hat{\mathbf{x}}$

In practice, convolution is usually linear (not wrapping around). Zero-padding both signals to length at least $n_1 + n_2 - 1$ makes circular and linear convolution agree. Sampling a signal too slowly makes high frequencies indistinguishable from low ones (aliasing): $e^{2\pi i(k + n)t/n}$ and $e^{2\pi ikt/n}$ agree at every sample point.

import numpy as np

n = 8
j, k = np.meshgrid(np.arange(n), np.arange(n), indexing="ij")
F = np.exp(-2j * np.pi * j * k / n)           # DFT matrix F_jk = w^{jk}, w = e^{-2 pi i / n}

x = np.random.default_rng(0).standard_normal(n)
print(np.allclose(F @ x, np.fft.fft(x)))      # True: np.fft.fft is multiplication by F, in O(n log n)
print(np.allclose(F.conj().T @ F, n * np.eye(n)))   # True: F / sqrt(n) is unitary
print(np.isclose(np.sum(x**2), np.sum(abs(F @ x)**2) / n))   # True: Parseval (energy is preserved)

# Every circulant matrix is diagonalized by F: C = F^{-1} diag(F c) F
c = np.array([2.0, -1.0, 0, 0, 0, 0, 0, -1.0])     # first column: periodic second difference
C = np.array([np.roll(c, s) for s in range(n)]).T  # column s is c shifted down by s
Lam = F @ C @ np.linalg.inv(F)
print(np.allclose(Lam, np.diag(np.diag(Lam))), np.allclose(np.diag(Lam), np.fft.fft(c)))   # True True

# Convolution theorem: circular convolution = pointwise product of spectra
h = np.array([0.25, 0.5, 0.25, 0, 0, 0, 0, 0])     # a smoothing kernel
direct = np.array([sum(h[m] * x[(t - m) % n] for m in range(n)) for t in range(n)])
via_fft = np.fft.ifft(np.fft.fft(h) * np.fft.fft(x)).real
print(np.allclose(direct, via_fft))           # True

# Denoising: project onto the low-frequency subspace
N = 256
t = np.arange(N) / N
clean = np.sin(2 * np.pi * 3 * t) + 0.5 * np.sin(2 * np.pi * 7 * t)
noisy = clean + 0.4 * np.random.default_rng(1).standard_normal(N)
X = np.fft.fft(noisy)
freqs = np.fft.fftfreq(N, d=1 / N)            # integer frequencies 0, 1, ..., -1
X[np.abs(freqs) > 10] = 0                     # keep 21 of 256 basis vectors
denoised = np.fft.ifft(X).real
rms = lambda e: np.sqrt(np.mean(e**2))
print(round(rms(noisy - clean), 3), round(rms(denoised - clean), 3))   # 0.367 0.069: about 5x less error

Key takeaways

• The DFT is the change of basis $\hat{\mathbf{x}} = F\mathbf{x}$, $F_{jk} = \omega^{jk}$, onto sampled complex sinusoids; $F/\sqrt{n}$ is unitary and preserves energy (Parseval).
• Every circulant matrix (circular convolution) is diagonalized by $F$, with eigenvalues $F\mathbf{c}$: convolution in time is pointwise multiplication in frequency.
• Filtering is a diagonal operation in the Fourier basis; denoising by discarding high frequencies is an orthogonal projection.
• The FFT factors $F$ into $\log_2n$ sparse butterfly stages, so the transform costs $O(n\log n)$.
• Zero-pad to get linear convolution from circular; sample fast enough to avoid aliasing.

🧪 Practice

1. Write $F_4$ explicitly (entries $\pm1, \pm i$) and compute the DFT of $(1, 0, 0, 0)$ and of $(1, 1, 1, 1)$. Interpret both results.
2. Prove that the vector $\mathbf{f}j = (\omega^{jk}){k}$ is an eigenvector of the cyclic shift matrix, and find its eigenvalue.
3. Multiply the polynomials $1 + 2x + 3x^2$ and $4 + 5x$ by zero-padding the coefficient vectors to length 4 and using np.fft. Compare with np.convolve.
4. Interview: You need to convolve a signal of length $10^6$ with a filter of length $10^4$. Compare the cost of direct convolution and FFT-based convolution, and say how you would handle a signal that streams in indefinitely. Hint: $O(nm)$ versus $O(n\log n)$; for streams, use overlap-add or overlap-save with blocks of moderate FFT size.

Machine Learning Foundations

Why it matters. Under the vocabulary of machine learning (features, weights, embeddings, layers, attention) are matrices and the operations of this guide. Datasets are matrices, models are compositions of linear maps and simple nonlinearities, training solves least-squares-like problems with gradients that are matrix products, and GPUs are fast because they multiply matrices. Knowing the linear algebra explains why models fail (collinearity, ill conditioning, rank collapse) and why the standard fixes work.

Analogy. If a neural network were a factory, the matrices are the machines that mix and route materials, and the nonlinearities are the quality-control gates between them. Without the gates, a long line of mixing machines could be replaced by a single machine.

The dictionary.

Machine learning concept	Linear algebra concept	Where in this guide
dataset ($n$ samples, $d$ features)	matrix $X \in \mathbb{R}^{n \times d}$	Section 4.1
linear regression	least squares, normal equations, QR	Section 9.5
collinear features	near-dependent columns, large $\kappa(X)$	Sections 6.3, 11.3
ridge regularization	$(X^TX + \lambda I)^{-1}X^T\mathbf{y}$: shifts eigenvalues up by $\lambda$	Section 10.3
PCA, dimensionality reduction	SVD of the centered data matrix	Section 11.2
embeddings, similarity search	vectors in $\mathbb{R}^d$, cosine of angles	Section 2.3
dense layer	affine map $W\mathbf{x} + \mathbf{b}$	Chapter 7
kernel methods	Gram matrix $K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j) \succeq 0$	Section 10.3
attention	row-stochastic matrix times value matrix	this topic, Section 13.2
LoRA fine-tuning	low-rank update $W + BA$ with $B \in \mathbb{R}^{d \times r}$, $r \ll d$	Section 11.2 (Eckart-Young)

Regression and conditioning. Least squares minimizes $\lVert X\mathbf{w} - \mathbf{y} \rVert^2$. When two features are nearly collinear, $X^TX$ has a tiny eigenvalue, and the solution can put huge weights of opposite sign on the twin features: they cancel on the training data but generalize badly. Ridge regression minimizes $\lVert X\mathbf{w} - \mathbf{y} \rVert^2 + \lambda\lVert\mathbf{w}\rVert^2$, with solution $(X^TX + \lambda I)^{-1}X^T\mathbf{y}$. Every eigenvalue of $X^TX$ is raised by $\lambda$, which bounds the condition number by $(\sigma_{\max}^2 + \lambda)/\lambda$. In SVD terms, ridge multiplies each component by $\sigma_i^2/(\sigma_i^2 + \lambda)$, so directions with small singular values are damped smoothly instead of amplified.

Gradients are matrix products. The gradient of $f(\mathbf{w}) = \lVert X\mathbf{w} - \mathbf{y} \rVert^2$ is $2X^T(X\mathbf{w} - \mathbf{y})$. Gradient descent needs only products with $X$ and $X^T$, never an inverse, which is how models with billions of parameters are trained. The step size must be below $2/L$, where $L$ is the largest eigenvalue of the Hessian $2X^TX$. The number of iterations needed grows with its condition number, the same phenomenon as for iterative solvers (Section 12.2). Backpropagation applies the chain rule layer by layer, and each step multiplies by a transposed weight matrix.

Why networks need nonlinearities. A composition of linear maps is linear: $W_3W_2W_1\mathbf{x} = (W_3W_2W_1)\mathbf{x}$. A deep network without activation functions can represent only what one matrix can, and its "depth" adds nothing but parameters. An elementwise nonlinearity such as $\operatorname{ReLU}(z) = \max(z, 0)$ breaks this collapse. Each layer is then piecewise linear, and a network of them can approximate any continuous function on a bounded region.

Attention. A transformer layer maps a sequence of token embeddings (rows of $E$) through three learned linear projections $Q = EW_Q$, $K = EW_K$, $V = EW_V$, then computes

$$\operatorname{Attention}(Q, K, V) = \operatorname{softmax}\Big(\frac{QK^T}{\sqrt{d_k}}\Big)V.$$

$QK^T$ is a matrix of dot products, a similarity score for every pair of tokens. The row-wise softmax turns it into a row-stochastic matrix, and multiplying by $V$ replaces each token by a weighted average of value vectors. The scale $1/\sqrt{d_k}$ keeps the dot products of random $d_k$-dimensional vectors at unit variance, so the softmax does not saturate.

import numpy as np

rng = np.random.default_rng(0)

# Data matrix: n samples (rows) x d features (columns); two features nearly collinear
n, d = 50, 3
X = rng.standard_normal((n, d))
X[:, 2] = X[:, 1] + 1e-3 * rng.standard_normal(n)   # feature 2 is almost a copy of feature 1
w_true = np.array([1.0, 2.0, 0.0])
y = X @ w_true + 0.1 * rng.standard_normal(n)

# Least squares (Section 9.5) is unstable here; ridge regularization fixes the conditioning
print(f"{np.linalg.cond(X.T @ X):.1e}")      # 3.3e+06: nearly singular normal equations
w_ls = np.linalg.lstsq(X, y, rcond=None)[0]
lam = 0.1
w_ridge = np.linalg.solve(X.T @ X + lam * np.eye(d), X.T @ y)
print(np.round(w_ls, 2), np.round(w_ridge, 2))   # [ 1.    6.27 -4.26] [1.   1.01 1.  ]
# least squares: huge opposite weights on the twin features; ridge: splits the weight 2 evenly
print(f"{np.linalg.cond(X.T @ X + lam * np.eye(d)):.1e}")   # 7.3e+02: adding lam I lifts the small eigenvalue

# Gradient descent on f(w) = ||Xw - y||^2 / n + lam ||w||^2 / n uses only matrix-vector products
w = np.zeros(d)
step = 1 / np.linalg.eigvalsh(2 * (X.T @ X + lam * np.eye(d)) / n).max()   # 1 / Lipschitz constant
for _ in range(20000):
    grad = 2 * (X.T @ (X @ w - y) + lam * w) / n
    w -= step * grad
print(np.allclose(w, w_ridge, atol=1e-4))     # True: same minimizer

# Without a nonlinearity, a deep network is just one matrix
W1, W2, W3 = rng.standard_normal((16, 3)), rng.standard_normal((16, 16)), rng.standard_normal((1, 16))
x = rng.standard_normal(3)
deep_linear = W3 @ (W2 @ (W1 @ x))
print(np.allclose(deep_linear, (W3 @ W2 @ W1) @ x))     # True: three "layers" collapse to one 1 x 3 matrix
relu = lambda z: np.maximum(z, 0)
print(np.round(W3 @ relu(W2 @ relu(W1 @ x)), 4))   # [4.2202]: a genuinely nonlinear function of x

# Scaled dot-product attention: softmax(Q K^T / sqrt(d_k)) V
tokens, d_model, d_k = 4, 8, 4
E = rng.standard_normal((tokens, d_model))    # one embedding per row
Wq, Wk, Wv = (rng.standard_normal((d_model, d_k)) for _ in range(3))
Q, K, V = E @ Wq, E @ Wk, E @ Wv              # three linear projections
scores = Q @ K.T / np.sqrt(d_k)               # tokens x tokens matrix of dot products
A = np.exp(scores - scores.max(axis=1, keepdims=True))
A /= A.sum(axis=1, keepdims=True)             # each row is a probability vector (row-stochastic)
out = A @ V                                   # each output is a weighted average of value vectors
print(A.shape, np.allclose(A.sum(axis=1), 1), out.shape)   # (4, 4) True (4, 4)

Key takeaways

• Datasets are matrices; regression is least squares; PCA is the SVD; embeddings are vectors compared by cosine similarity.
• Collinear features make $X^TX$ ill conditioned; ridge adds $\lambda I$, raising every eigenvalue and damping small-singular-value directions.
• Training uses gradients such as $2X^T(X\mathbf{w} - \mathbf{y})$, which need only matrix-vector products; convergence speed depends on conditioning.
• Without nonlinearities, a deep network collapses to a single matrix.
• Attention is $\operatorname{softmax}(QK^T/\sqrt{d_k})V$: pairwise dot products, a row-stochastic weight matrix, and a weighted average of values.

🧪 Practice

1. Derive the ridge solution by setting the gradient of $\lVert X\mathbf{w} - \mathbf{y} \rVert^2 + \lambda\lVert\mathbf{w}\rVert^2$ to zero. Why is $X^TX + \lambda I$ invertible for every $\lambda > 0$, even when $X^TX$ is singular?
2. Using the SVD $X = U\Sigma V^T$, show that the ridge prediction is $X\mathbf{w}_\lambda = \sum_i\frac{\sigma_i^2}{\sigma_i^2 + \lambda}(\mathbf{u}_i^T\mathbf{y})\mathbf{u}_i$. What happens as $\lambda \to 0$ and as $\lambda \to \infty$?
3. A dense layer maps $\mathbb{R}^{4096} \to \mathbb{R}^{4096}$. How many parameters does a rank-8 LoRA update $BA$ have, compared with a full update of $W$? What subspace can the update change?
4. Interview: Attention costs $O(n^2d)$ for a sequence of length $n$. Which matrix causes the quadratic cost, and name one linear-algebra idea used to reduce it. Hint: the $n \times n$ score matrix $QK^T$; low-rank approximations (Linformer), kernel feature maps that reorder the product as $Q(K^TV)$, or computing it in blocks without storing it (FlashAttention).

Quantum State Spaces

Why it matters. Quantum mechanics is linear algebra over $\mathbb{C}$, with a probability rule attached. States are vectors, physical evolution is unitary matrices, measurable quantities are Hermitian matrices, and composite systems are tensor products (Section 13.1). Quantum computing algorithms are sequences of matrix multiplications on vectors whose dimension doubles with every qubit. With the right linear algebra, the "strange" phenomena of superposition, interference, and entanglement become concrete facts about vectors.

Analogy. A classical bit is a light switch: on or off. A qubit is more like an arrow in a two-dimensional complex space, pointing partly along "off" and partly along "on." Measuring the qubit forces the arrow onto one axis, at random, with probabilities given by the squared lengths of its components along each axis.

The postulates in linear algebra terms.

Physical notion	Mathematical object	Key property
state of a system	unit vector $\lvert\psi\rangle$ in a complex Hilbert space $\mathcal{H}$	$\lvert\psi\rangle$ and $e^{i\theta}\lvert\psi\rangle$ are the same state
evolution (gates)	unitary matrix $U$	preserves norms, hence total probability
observable	Hermitian matrix $A$	real eigenvalues = possible outcomes; orthogonal eigenvectors
measurement probability (Born rule)	$\lvert\langle\mathbf{e}_k\vert\psi\rangle\rvert^2$	probabilities of all outcomes sum to 1
expected value	$\langle\psi\vert A\vert\psi\rangle$	a Hermitian form, so it is real
composite system	tensor product $\mathcal{H}_A \otimes \mathcal{H}_B$	dimensions multiply

Dirac notation. A ket $\lvert\psi\rangle$ is a column vector, and a bra $\langle\psi\rvert$ is its conjugate transpose, a row vector. Then $\langle\phi\vert\psi\rangle = \boldsymbol{\phi}^*\boldsymbol{\psi}$ is the inner product, conjugate-linear in the first slot as in Section 13.1, and $\lvert\psi\rangle\langle\psi\rvert$ is the rank-1 orthogonal projection onto $\lvert\psi\rangle$.

Qubits and gates. A qubit lives in $\mathbb{C}^2$ with basis $\lvert0\rangle = (1, 0)$ and $\lvert1\rangle = (0, 1)$. A general state is $\alpha\lvert0\rangle + \beta\lvert1\rangle$ with $\lvert\alpha\rvert^2 + \lvert\beta\rvert^2 = 1$; measuring gives 0 with probability $\lvert\alpha\rvert^2$ and 1 with probability $\lvert\beta\rvert^2$. Common gates:

$$X = \begin{bmatrix} 0 & 1 \ 1 & 0 \end{bmatrix}, \qquad Z = \begin{bmatrix} 1 & 0 \ 0 & -1 \end{bmatrix}, \qquad H = \frac{1}{\sqrt{2}}\begin{bmatrix} 1 & 1 \ 1 & -1 \end{bmatrix}.$$

The Hadamard gate $H$ creates equal superpositions: $H\lvert0\rangle = (\lvert0\rangle + \lvert1\rangle)/\sqrt{2}$ and $H\lvert1\rangle = (\lvert0\rangle - \lvert1\rangle)/\sqrt{2}$. These two states give identical measurement statistics in the standard basis. They differ by a relative phase, which another $H$ converts into a certain outcome. This is interference: amplitudes, unlike probabilities, can cancel.

Many qubits and entanglement. An $n$-qubit system lives in $(\mathbb{C}^2)^{\otimes n} = \mathbb{C}^{2^n}$. Simulating 50 qubits classically means storing $2^{50} \approx 10^{15}$ complex amplitudes, which is where the potential power of quantum computers comes from. A two-qubit state $\sum_{ij}c_{ij}\lvert i\rangle \otimes \lvert j\rangle$ is a product state if it equals $\lvert a\rangle \otimes \lvert b\rangle$. By the tensor product topic, this happens exactly when the $2 \times 2$ coefficient matrix $C = [c_{ij}]$ has rank 1. Otherwise the state is entangled. The SVD of $C$ is the Schmidt decomposition $\lvert\psi\rangle = \sum_k\sigma_k\lvert u_k\rangle \otimes \lvert v_k\rangle$, and the number of nonzero singular values (the Schmidt rank) measures entanglement. The Bell state $(\lvert00\rangle + \lvert11\rangle)/\sqrt{2}$ has $C = I/\sqrt{2}$, rank 2: measuring either qubit instantly determines the other.

Mixed states. When a state is only known statistically, it is described by a density matrix $\rho = \sum_ip_i\lvert\psi_i\rangle\langle\psi_i\rvert$, which is Hermitian, positive semidefinite, and has trace 1. Expected values become $\operatorname{tr}(\rho A)$.

import numpy as np

ket0, ket1 = np.array([1, 0], complex), np.array([0, 1], complex)
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)  # Hadamard gate
X = np.array([[0, 1], [1, 0]])                # NOT gate (Pauli X)
Z = np.diag([1, -1])                          # Pauli Z: the observable "which basis state?"

psi = H @ ket0                                # (|0> + |1>) / sqrt(2)
print(np.round(psi, 4), np.round(np.abs(psi) ** 2, 4))   # [0.7071+0.j 0.7071+0.j] [0.5 0.5]: Born rule
print(np.allclose(H.conj().T @ H, np.eye(2))) # True: gates are unitary, so norms (total probability) are preserved
print(np.round(np.vdot(psi, Z @ psi).real, 4), np.round(np.vdot(ket0, Z @ ket0).real, 4))   # 0.0 1.0: <psi|Z|psi>

# Global phase is unobservable; relative phase is not
phased = np.exp(1j * 0.7) * psi
minus = H @ ket1                              # (|0> - |1>) / sqrt(2): same probabilities as psi in this basis
print(np.allclose(np.abs(phased) ** 2, np.abs(psi) ** 2), np.round(np.abs(H @ psi) ** 2, 4), np.round(np.abs(H @ minus) ** 2, 4))
# True [1. 0.] [0. 1.]: after another H, the relative phase decides the outcome with certainty

# Two qubits live in C^2 (x) C^2 = C^4, basis |00>, |01>, |10>, |11>
CNOT = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
bell = CNOT @ np.kron(H, np.eye(2)) @ np.kron(ket0, ket0)   # Hadamard on qubit 1, then CNOT
print(np.round(bell.real, 4))                 # [0.7071 0.     0.     0.7071]: (|00> + |11>) / sqrt(2)

# Entanglement test: reshape the 4 amplitudes into a 2x2 matrix; product state <=> rank 1
def schmidt(state):
    return np.round(np.linalg.svd(state.reshape(2, 2), compute_uv=False), 4)   # Schmidt coefficients
product = np.kron(H @ ket0, ket1)             # |+> (x) |1>
print(schmidt(product), schmidt(bell))        # [1. 0.] [0.7071 0.7071]: rank 1 vs rank 2 (entangled)

# Measuring both qubits of the Bell state gives 00 or 11, never 01 or 10: perfectly correlated
print(np.round(np.abs(bell) ** 2, 4))         # [0.5 0.  0.  0.5]

Key takeaways

• States are unit vectors in $\mathbb{C}^n$ up to global phase; gates are unitary; observables are Hermitian, with eigenvalues as outcomes.
• Born rule: outcome $k$ has probability $\lvert\langle\mathbf{e}_k\vert\psi\rangle\rvert^2$; expected values are $\langle\psi\vert A\vert\psi\rangle$.
• Relative phases are invisible to one measurement but cause interference after further gates.
• Composite systems are tensor products, so dimension grows as $2^n$; product states are rank-1 coefficient matrices.
• Entanglement is measured by the SVD of the coefficient matrix (Schmidt decomposition); density matrices are PSD with trace 1.

🧪 Practice

1. Verify that $H$, $X$, and $Z$ are unitary and Hermitian. Compute $HXH$ and $HZH$, and interpret the results.
2. For $\lvert\psi\rangle = \tfrac{1}{\sqrt{3}}\lvert0\rangle + \sqrt{\tfrac{2}{3}},i\lvert1\rangle$, compute the probabilities of measuring 0 and 1 and the expected value $\langle\psi\vert Z\vert\psi\rangle$.
3. Decide whether $\tfrac12(\lvert00\rangle + \lvert01\rangle + \lvert10\rangle + \lvert11\rangle)$ and $\tfrac{1}{\sqrt{2}}(\lvert01\rangle - \lvert10\rangle)$ are entangled by computing the rank of their coefficient matrices.
4. Interview: Why can't a quantum computer simply copy an unknown qubit (the no-cloning theorem)? Give the linear-algebra argument. Hint: suppose a unitary $U$ maps $\lvert\psi\rangle \otimes \lvert0\rangle \mapsto \lvert\psi\rangle \otimes \lvert\psi\rangle$ for all $\lvert\psi\rangle$; apply it to a superposition and compare with linearity (or compare inner products before and after).

<a id="134-connections-to-further-mathematics"></a>

13.4 Connections to Further Mathematics

Linear algebra is the entry point to much of modern mathematics. This subchapter previews five directions that start by relaxing one assumption of this guide (scalars from a field, finite dimension, a flat space) or by applying linear algebra to symmetry. In each case, it shows what survives and what breaks.

Modules over Rings

Why it matters. Every vector space in this guide had scalars from a field, where every nonzero scalar has an inverse. Many natural problems have scalars that cannot be divided: integer linear systems (Diophantine equations), abelian groups, lattices in cryptography, and even a single linear operator $T$ viewed as acting through polynomials $p(T)$. A module is a "vector space over a ring." Studying modules explains which theorems of linear algebra really needed division, and it gives a single theorem that contains both the classification of finite abelian groups and the Jordan canonical form.

Analogy. Working in a field is like paying with a currency that can be split into any fraction. Working over the integers is like paying only in whole coins. Many transactions still work, but some amounts cannot be made from the coins available, and the familiar rules about change need rethinking.

Definition. Let $R$ be a commutative ring with 1. An $R$-module $M$ is an abelian group with a scalar multiplication $R \times M \to M$ satisfying the same axioms as a vector space (Section 6.1): distributivity in both arguments, $(rs)m = r(sm)$, and $1m = m$. Examples:

• If $R$ is a field, $R$-modules are exactly vector spaces.
• $\mathbb{Z}$-modules are exactly abelian groups, with $nm = m + \cdots + m$.
• $R^n$, with componentwise operations, is the free module of rank $n$.
• $\mathbb{Z}/n\mathbb{Z}$ is a $\mathbb{Z}$-module in which $n \cdot m = 0$ for every $m$.
• A vector space $V$ with an operator $T$ is an $F[x]$-module, where a polynomial acts by $p(x) \cdot \mathbf{v} = p(T)\mathbf{v}$.

What breaks.

Vector spaces (over a field)	Modules over a ring
every space has a basis	many modules have none: $\mathbb{Z}/6$ has no nonempty independent set over $\mathbb{Z}$
$a\mathbf{v} = \mathbf{0}$ with $a \ne 0$ implies $\mathbf{v} = \mathbf{0}$	torsion: $3 \cdot [2] = [0]$ in $\mathbb{Z}/6$
every independent set extends to a basis	${2} \subset \mathbb{Z}$ is independent, but no basis of $\mathbb{Z}$ contains it
every subspace has a complement	$2\mathbb{Z} \subset \mathbb{Z}$ has no complementary submodule
dimension classifies spaces	rank alone is not enough (torsion must be described too)
$A$ invertible $\iff \det A \ne 0$	$A$ invertible over $R$ $\iff \det A$ is a unit of $R$ ($\pm1$ over $\mathbb{Z}$)

The structure theorem. Over a principal ideal domain (PID) such as $\mathbb{Z}$ or $F[x]$, much is recovered. Every finitely generated module decomposes as

$$M \cong R^r \oplus R/(d_1) \oplus R/(d_2) \oplus \cdots \oplus R/(d_k), \qquad d_1 \mid d_2 \mid \cdots \mid d_k,$$

with a free part of rank $r$ and a torsion part described by the invariant factors $d_i$, both unique.

Smith normal form: the computation. A finitely generated module is presented by generators and a matrix of relations. Over a PID, any integer matrix $A$ can be brought to diagonal form by unimodular row and column operations (integer matrices with determinant $\pm1$, invertible over $\mathbb{Z}$):

$$UAV = \operatorname{diag}(d_1, d_2, \dots, d_k, 0, \dots, 0), \qquad d_1 \mid d_2 \mid \cdots.$$

The algorithm is Gaussian elimination in which division is replaced by the Euclidean algorithm: to clear an entry, repeatedly subtract integer multiples until a gcd sits in the pivot position. The $d_i$ are the invariant factors, so the module is read off directly. Applications:

• Integer linear systems. $A\mathbf{x} = \mathbf{b}$ has an integer solution if and only if the transformed right side $U\mathbf{b}$ is divisible entrywise by the $d_i$. Rank conditions over $\mathbb{Q}$ are not enough: $6x + 10y = 1$ has rational solutions but no integer ones.
• Finite abelian groups. Every finite abelian group is $\mathbb{Z}/d_1 \oplus \cdots \oplus \mathbb{Z}/d_k$.
• Canonical forms. For $V$ as an $F[x]$-module via $T$, the decomposition gives the rational canonical form (Section 8.3), with the minimal polynomial as $d_k$. Over $\mathbb{C}$, splitting each $d_i$ into powers $(x - \lambda)^m$ gives the Jordan canonical form, one block per factor.
• Topology. Homology groups of a space are computed as the Smith normal forms of integer boundary matrices.

from sympy import Matrix, ZZ, gcd
from sympy.matrices.normalforms import smith_normal_form, invariant_factors

# Z-module (abelian group) with generators a, b, c and relations given by the ROWS of R:
#   2a + 4b + 4c = 0,  -6a + 6b + 12c = 0,  10a - 4b - 16c = 0
R = Matrix([[2, 4, 4], [-6, 6, 12], [10, -4, -16]])
print(smith_normal_form(R, domain=ZZ))        # Matrix([[2, 0, 0], [0, 6, 0], [0, 0, 12]])
print(invariant_factors(R, domain=ZZ))        # (2, 6, 12): each divides the next
# So the group Z^3 / (row space of R) is  Z/2 + Z/6 + Z/12,  of order |det R| = 144
print(abs(R.det()))                           # 144

# A relation matrix with a zero invariant factor has a free part
R2 = Matrix([[2, 4], [3, 6]])                 # rank 1 over Q
print(smith_normal_form(R2, domain=ZZ))       # Matrix([[1, 0], [0, 0]]): Z/1 + Z = Z (free of rank 1)

# Integer linear systems: 6x + 10y + 15z = 1 has integer solutions iff the gcd of coefficients divides 1
print(gcd(gcd(6, 10), 15))                    # 1: solvable over Z (e.g. x = 1, y = 1, z = -1)
print(6 * 1 + 10 * 1 + 15 * (-1))             # 1
# but 6x + 10y = 1 has none: gcd(6, 10) = 2 does not divide 1, although it is solvable over Q

# Torsion breaks familiar facts: in Z/6 (a Z-module), 3 * [2] = [0] with 3 != 0 and [2] != 0
print((3 * 2) % 6)                            # 0: the set {[2]} is linearly DEPENDENT over Z

Key takeaways

• A module is a vector space whose scalars form a ring; $\mathbb{Z}$-modules are abelian groups, and an operator $T$ makes $V$ an $F[x]$-module.
• Without division, bases may fail to exist, torsion appears, and independent sets need not extend to bases.
• Over a PID, finitely generated modules are $R^r \oplus R/(d_1) \oplus \cdots \oplus R/(d_k)$ with $d_1 \mid \cdots \mid d_k$.
• The Smith normal form computes the $d_i$ using elimination with the Euclidean algorithm and unimodular operations.
• One theorem yields the classification of finite abelian groups, integer solvability of linear systems, and the rational and Jordan canonical forms.

🧪 Practice

1. Show that $\mathbb{Q}$ is not a free $\mathbb{Z}$-module. Hint: any two rationals are dependent over $\mathbb{Z}$, so a basis would have one element.
2. Compute the Smith normal form of $\begin{bmatrix} 4 & 6 \ 6 & 4 \end{bmatrix}$ by hand. Which abelian group has this relation matrix, and what is its order?
3. Decide whether $2x + 4y = 3$ and $3x + 5y = 7$, $x + y = 2$ have integer solutions, by comparing gcd and Smith normal form arguments with the rational solutions.
4. Interview: Explain how the Jordan canonical form of an operator $T$ on $\mathbb{C}^n$ follows from the structure theorem for modules over a PID. Hint: make $\mathbb{C}^n$ a $\mathbb{C}[x]$-module via $x \cdot \mathbf{v} = T\mathbf{v}$; it is all torsion; split each invariant factor into prime powers $(x - \lambda)^m$, and find the matrix of $T$ on $\mathbb{C}[x]/((x - \lambda)^m)$ in the basis $(x - \lambda)^{m-1}, \dots, 1$.

Group Representations

Why it matters. Symmetry simplifies problems: a molecule's symmetry predicts which vibrations are infrared-active, the rotational symmetry of an atom explains its energy levels, and the symmetry of a ring of $n$ points is why the DFT diagonalizes circulant matrices (Section 13.3). Representation theory makes symmetry computable by realizing the elements of an abstract group as matrices. Decomposing that action into its smallest invariant pieces is block diagonalization, the same idea as eigendecomposition, done for a whole group of matrices at once.

Analogy. A group is a script of moves (rotate, reflect, swap). A representation is a cast of actors (matrices) performing the script on a stage (a vector space). Different casts can perform the same script, and each performance breaks into independent scenes (irreducible subspaces) that never interact.

Definition. A representation of a group $G$ on a vector space $V$ is a homomorphism $\rho : G \to GL(V)$: each $g \in G$ becomes an invertible matrix $\rho(g)$, with $\rho(gh) = \rho(g)\rho(h)$ and $\rho(e) = I$. Its dimension is $\dim V$. Examples:

• Trivial: $\rho(g) = 1$ on $\mathbb{C}$.
• Rotations: the cyclic group $\mathbb{Z}_n$ acting on $\mathbb{R}^2$ by $\rho(k) = $ rotation by $2\pi k/n$.
• Characters of $\mathbb{Z}_n$: $\rho_j(k) = \omega^{jk}$ on $\mathbb{C}$, for $j = 0, \dots, n - 1$. These are the Fourier modes.
• Permutations: $S_n$ acting on $\mathbb{C}^n$ by permutation matrices. Sign: $\rho(\sigma) = \operatorname{sgn}(\sigma) = \det(\text{permutation matrix})$.

Invariant subspaces and irreducibility. A subspace $W$ is invariant if $\rho(g)W \subseteq W$ for all $g$, a common invariant subspace for a whole family of matrices (Section 8.3). A representation with no invariant subspaces other than $0$ and $V$ is irreducible. For permutations of $\mathbb{C}^3$, the line spanned by $(1, 1, 1)$ is invariant, since permuting coordinates does not change it, and so is the plane $x_1 + x_2 + x_3 = 0$.

Maschke's theorem. For a finite group over $\mathbb{C}$ (or $\mathbb{R}$), every representation is a direct sum of irreducible ones. The proof is an averaging trick. Start from any inner product and define $\langle\mathbf{u}, \mathbf{v}\rangle_G = \frac{1}{\lvert G \rvert}\sum_g\langle\rho(g)\mathbf{u}, \rho(g)\mathbf{v}\rangle$. Every $\rho(g)$ is unitary for this inner product, so the orthogonal complement of an invariant subspace is again invariant (Section 9.2). Split repeatedly until only irreducible pieces remain. In a suitable basis, every $\rho(g)$ is then block diagonal with the same block pattern.

Schur's lemma. If $\rho$ and $\sigma$ are irreducible and $TA = BT$ for all group elements (with $A = \rho(g)$, $B = \sigma(g)$), then $T = 0$ or $T$ is invertible. Over $\mathbb{C}$, a matrix that commutes with every $\rho(g)$ of an irreducible $\rho$ is a scalar multiple of $I$. This is the representation-theoretic version of "commuting matrices share eigenvectors."

Characters. The character of $\rho$ is $\chi_\rho(g) = \operatorname{tr}\rho(g)$. Because trace is a similarity invariant (Section 8.2), $\chi$ is constant on conjugacy classes, and isomorphic representations have equal characters. With the inner product $\langle\chi, \psi\rangle = \frac{1}{\lvert G \rvert}\sum_g\overline{\chi(g)}\psi(g)$:

• Characters of distinct irreducibles are orthonormal.
• A representation is irreducible if and only if $\langle\chi, \chi\rangle = 1$.
• The multiplicity of irreducible $\pi$ in $\rho$ is $\langle\chi_\pi, \chi_\rho\rangle$.
• There are as many irreducibles as conjugacy classes, and $\sum_\pi(\dim\pi)^2 = \lvert G \rvert$.

The character table of $S_3$ (classes: identity, 3 transpositions, 2 three-cycles):

Irreducible	$e$	$(12)$	$(123)$	Dimension
trivial	1	1	1	1
sign	1	$-1$	1	1
standard	2	0	$-1$	2
permutation rep on $\mathbb{C}^3$	3	1	0	3 = trivial + standard

Fourier analysis is representation theory. For the abelian group $\mathbb{Z}_n$, all irreducibles are one-dimensional, namely the $n$ characters $k \mapsto \omega^{jk}$. Decomposing the regular representation (cyclic shifts on $\mathbb{C}^n$) into irreducibles is exactly the DFT, and Schur's lemma explains why every shift-invariant (circulant) operator becomes diagonal. For non-abelian groups, the same program yields block-diagonal "Fourier transforms," used in chemistry, crystallography, and quantum mechanics.

import numpy as np
from itertools import permutations

# The permutation representation of S_3 on C^3: rho(g) e_i = e_{g(i)}
G = list(permutations(range(3)))              # 6 group elements, as tuples g = (g(0), g(1), g(2))
def rho(g):
    P = np.zeros((3, 3))
    for i in range(3):
        P[g[i], i] = 1
    return P
compose = lambda g, h: tuple(g[h[i]] for i in range(3))   # (gh)(i) = g(h(i))

print(all(np.array_equal(rho(compose(g, h)), rho(g) @ rho(h)) for g in G for h in G))   # True: a homomorphism

# Averaging over the group projects onto the vectors fixed by every group element
P_triv = sum(rho(g) for g in G) / len(G)
print(np.round(P_triv, 4)[0])                 # [0.3333 0.3333 0.3333]: all entries 1/3, projection onto span{(1, 1, 1)}
print(np.allclose(P_triv @ P_triv, P_triv))   # True: it is a projection
# Its complement {x : x1 + x2 + x3 = 0} is also invariant: C^3 = trivial (+) 2-dim "standard" rep

# Characters: chi(g) = trace rho(g); depends only on the conjugacy class (cycle type)
chi = np.array([np.trace(rho(g)) for g in G])
print({g: int(c) for g, c in zip(G, chi)})
# {(0, 1, 2): 3, (0, 2, 1): 1, (1, 0, 2): 1, (1, 2, 0): 0, (2, 0, 1): 0, (2, 1, 0): 1}
# identity: 3 (fixes all points), transpositions: 1, 3-cycles: 0
inner = lambda a, b: np.vdot(a, b).real / len(G)   # <a, b> = (1/|G|) sum conj(a(g)) b(g)
print(inner(chi, chi))                        # 2.0: a sum of 2 distinct irreducibles (1 would mean irreducible)

# Character table of S_3: trivial, sign, standard (2-dim)
def sign(g):
    return round(np.linalg.det(rho(g)))
triv = np.ones(6)
sgn = np.array([sign(g) for g in G], float)
std = chi - triv                              # the character of the 2-dim complement
for name, c in [("trivial", triv), ("sign", sgn), ("standard", std)]:
    print(name, inner(c, c), inner(c, chi))   # each irreducible has norm 1; multiplicity in chi
# trivial 1.0 1.0
# sign 1.0 0.0
# standard 1.0 1.0
print(1**2 + 1**2 + 2**2 == len(G))           # True: sum of squared dimensions = |G|

Key takeaways

• A representation is a homomorphism $\rho : G \to GL(V)$; it turns symmetry into matrices.
• Irreducible representations have no common invariant subspaces; by Maschke, finite-group representations split into irreducibles, so all $\rho(g)$ block-diagonalize at once.
• Averaging over $G$ makes $\rho$ unitary and projects onto invariant vectors.
• Characters $\chi(g) = \operatorname{tr}\rho(g)$ are class functions; their orthogonality decides irreducibility and multiplicities.
• The DFT is the decomposition of the regular representation of $\mathbb{Z}_n$; Schur's lemma is why it diagonalizes circulants.

🧪 Practice

1. Check that $k \mapsto \begin{bmatrix} \cos(2\pi k/4) & -\sin(2\pi k/4) \ \sin(2\pi k/4) & \cos(2\pi k/4) \end{bmatrix}$ is a representation of $\mathbb{Z}_4$. Is it irreducible over $\mathbb{R}$? Over $\mathbb{C}$? Hint: look for common eigenvectors.
2. Using the character table, decompose the representation of $S_3$ on $\mathbb{C}^3 \otimes \mathbb{C}^3$ (character $\chi^2$) into irreducibles.
3. Prove that for a one-dimensional representation, $\rho(g)$ is a root of unity of order dividing the order of $g$.
4. Interview: A $6 \times 6$ symmetric matrix $M$ describes coupled oscillators arranged on a hexagon with all couplings equal, so $M$ commutes with the rotation of the hexagon. How can you find its eigenvectors without solving a degree-6 characteristic polynomial? Hint: $M$ is circulant; Schur's lemma and the characters of $\mathbb{Z}_6$ say the Fourier vectors are eigenvectors, and the eigenvalues are the DFT of its first row.

Infinite-Dimensional Spaces and Hilbert Spaces

Why it matters. Functions are vectors (Section 6.1): signals, wavefunctions, temperature profiles, probability densities. Most spaces of functions are infinite-dimensional, and many results of this guide assumed finite dimension (rank-nullity, existence of eigenvalues, "injective implies surjective"). Hilbert spaces are the infinite-dimensional spaces where geometry still works: orthogonal projection, orthonormal bases, and least squares all carry over, provided the space is complete. Fourier series (Section 9.5), quantum mechanics, PDE theory, and kernel methods in machine learning all take place in Hilbert spaces.

Analogy. Infinite dimensions are like an infinite hotel. Some finite-hotel common sense survives (every room has a number), and some fails dramatically: a fully occupied hotel can take in a new guest by moving everyone up one room, which is an injective map that is not surjective.

Why algebra alone is not enough. Every vector space has a Hamel basis, but in infinite dimensions it is uncountable and cannot be written down (its existence needs the axiom of choice). A Hamel basis allows only finite linear combinations, and the useful expansions, such as Fourier series $\sum_{n=1}^\infty c_n\mathbf{e}_n$, are infinite. Infinite sums require a notion of convergence, so the space needs a norm, and limits that should exist must actually exist:

• A Banach space is a complete normed space: every Cauchy sequence converges.
• A Hilbert space is a complete inner product space.

Main examples.

Space	Elements	Inner product	Role
$\ell^2$	sequences with $\sum\lvert x_n\rvert^2 < \infty$	$\sum\overline{x_n}y_n$	the model separable Hilbert space
$L^2[a, b]$	(classes of) functions with $\int\lvert f\rvert^2 < \infty$	$\int\overline{f}g$	signals, Fourier series, quantum wavefunctions
Sobolev $H^1$	$L^2$ functions with $L^2$ derivatives	$\int(\overline{f}g + \overline{f'}g')$	weak solutions of PDEs, finite elements
$C[a, b]$ with the sup norm	continuous functions	none (a Banach space)	uniform approximation
$C[a, b]$ with the $L^2$ norm	continuous functions	$\int\overline{f}g$	not complete: its completion is $L^2$

Orthonormal bases. An orthonormal set ${\mathbf{e}_n}$ in a Hilbert space $H$ is a basis (complete orthonormal system) if its finite linear combinations are dense. Then every $\mathbf{x} \in H$ satisfies

$$\mathbf{x} = \sum_{n=1}^\infty\langle\mathbf{e}_n, \mathbf{x}\rangle\mathbf{e}n, \qquad \lVert\mathbf{x}\rVert^2 = \sum{n=1}^\infty\lvert\langle\mathbf{e}_n, \mathbf{x}\rangle\rvert^2 \quad \text{(Parseval)}.$$

For any orthonormal set, complete or not, Bessel's inequality $\sum\lvert\langle\mathbf{e}_n, \mathbf{x}\rangle\rvert^2 \le \lVert\mathbf{x}\rVert^2$ holds. The partial sum $S_N$ is the orthogonal projection onto $\operatorname{span}{\mathbf{e}_1, \dots, \mathbf{e}_N}$, the best approximation from that subspace, exactly as in Chapter 9. The functions $e^{inx}/\sqrt{2\pi}$, $n \in \mathbb{Z}$, form an orthonormal basis of $L^2[-\pi, \pi]$. The Legendre polynomials, obtained by Gram-Schmidt on $1, x, x^2, \dots$ (Section 9.3), form one for $L^2[-1, 1]$. Every infinite-dimensional separable Hilbert space has a countable orthonormal basis, and the coefficient map makes it isometrically isomorphic to $\ell^2$, the infinite-dimensional analog of "every $n$-dimensional space is $F^n$."

What still works.

• Projection theorem: for every closed subspace $W$, $H = W \oplus W^\perp$, and the closest point in $W$ exists and is unique. Least squares works in infinite dimensions.
• Riesz representation: every continuous linear functional is $\varphi(\mathbf{x}) = \langle\mathbf{y}, \mathbf{x}\rangle$ for a unique $\mathbf{y}$, so $H \cong H^*$ (Section 7.4).
• Cauchy-Schwarz, the triangle inequality, and Pythagoras hold unchanged.

What fails.

Finite dimensions	Infinite dimensions
every subspace is closed	subspaces can be dense but proper (polynomials in $L^2[0, 1]$)
all norms are equivalent	different norms give different convergence ($L^2$ vs uniform, see the code)
the closed unit ball is compact	the orthonormal sequence $\mathbf{e}_n$ has no convergent subsequence ($\lVert\mathbf{e}_n - \mathbf{e}_m\rVert = \sqrt{2}$)
linear maps are continuous	differentiation is linear but discontinuous on $L^2$ (next topic)
injective $\iff$ surjective for $T : V \to V$	the shift $(x_1, x_2, \dots) \mapsto (0, x_1, x_2, \dots)$ is injective, not surjective

import numpy as np

# f(x) = x on [-pi, pi], expanded in the orthonormal family sin(nx)/sqrt(pi) of L^2[-pi, pi]
x = np.linspace(-np.pi, np.pi, 200001)
f = x
inner = lambda g, h: np.trapezoid(g * h, x)   # <g, h> = integral of g h over [-pi, pi]
e = lambda n: np.sin(n * x) / np.sqrt(np.pi)  # unit vectors: <e_n, e_m> = delta_nm

norm2 = inner(f, f)
print(round(norm2, 6), round(2 * np.pi**3 / 3, 6))   # 20.670851 20.670851: ||f||^2 = 2 pi^3 / 3
print(round(inner(f, e(3)), 6), round(2 * np.sqrt(np.pi) / 3, 6))   # 1.181636 1.181636: c_3, numeric vs exact

# Partial sums S_N = sum_{n <= N} c_n e_n are orthogonal projections onto span{e_1..e_N}
for N in (1, 5, 25, 125):
    c = np.array([2 * np.sqrt(np.pi) * (-1) ** (n + 1) / n for n in range(1, N + 1)])   # exact c_n
    SN = sum(c[n - 1] * e(n) for n in range(1, N + 1))
    print(N, round(inner(f - SN, f - SN), 5), round(norm2 - np.sum(c**2), 5), round(np.max(np.abs(f - SN)), 3))
# 1 8.10448 8.10448 3.142
# 5 2.27857 2.27857 3.142
# 25 0.49274 0.49274 3.142
# 125 0.10013 0.10013 3.142
# ||f - S_N||^2 = ||f||^2 - sum c_n^2 (Pythagoras) -> 0, yet the max error stays pi (at x = +-pi):
# convergence in the L^2 norm, not uniform convergence

# Parseval: sum_n c_n^2 = ||f||^2, i.e. 4 pi * sum 1/n^2 = 2 pi^3 / 3, i.e. sum 1/n^2 = pi^2 / 6
print(round(4 * np.pi * np.pi**2 / 6, 6))     # 20.670851

Key takeaways

• Infinite expansions need a norm and completeness: Banach spaces are complete normed spaces, and Hilbert spaces are complete inner product spaces.
• $\ell^2$ and $L^2$ are the central examples; all separable infinite-dimensional Hilbert spaces are isometric to $\ell^2$.
• Orthonormal bases give expansions $\sum\langle\mathbf{e}_n, \mathbf{x}\rangle\mathbf{e}_n$ with Parseval's identity; partial sums are orthogonal projections.
• Projections onto closed subspaces and the Riesz representation theorem carry over from finite dimensions.
• Compactness of the unit ball, equivalence of norms, continuity of linear maps, and "injective iff surjective" all fail.

🧪 Practice

1. Show that $\mathbf{x} = (1, \tfrac12, \tfrac13, \dots)$ is in $\ell^2$ but not in $\ell^1$, and compute $\lVert\mathbf{x}\rVert_2$ exactly.
2. Find a sequence of continuous functions on $[0, 1]$ that is Cauchy in the $L^2$ norm but whose limit is discontinuous. Conclude that $C[0, 1]$ with the $L^2$ norm is not complete.
3. Apply Gram-Schmidt in $L^2[-1, 1]$ to $1, x, x^2$, and check that the results are multiples of the Legendre polynomials $1$, $x$, $\tfrac12(3x^2 - 1)$.
4. Interview: Why do Fourier partial sums of a square wave converge in $L^2$ but always overshoot by about 9% near the jump (the Gibbs phenomenon)? What does that say about the difference between norms in infinite dimensions? Hint: $L^2$ convergence controls the integral of the squared error, which can go to zero while the maximum error does not; the overshoot narrows but never shrinks.

Operator Theory Preview

Why it matters. Differential and integral equations are linear equations $T\mathbf{u} = \mathbf{f}$ in which $T$ acts on functions: $T = -\frac{d^2}{dx^2}$ for heat and vibration, $T = -\frac{\hbar^2}{2m}\nabla^2 + V$ (the Schrödinger operator) in quantum mechanics, and integral operators in potential theory and machine learning kernels. Operator theory asks the questions of Chapters 7, 8, and 10 for such $T$: when is it invertible, what is its spectrum, is there an orthonormal eigenbasis? Every numerical method for these equations replaces $T$ by a matrix, so knowing which properties survive that replacement is essential for trusting the computed answer.

Analogy. A matrix is a finite machine with a list of eigenvalues on its label. An operator is a machine with infinitely many settings. Some operators (compact ones) behave like very long matrices whose eigenvalues fade toward zero. Others have "eigenvalues" smeared into a continuum, or settings that blow up any input.

Bounded operators. A linear map $T : H \to H$ on a Hilbert space is bounded if $\lVert T \rVert = \sup_{\lVert\mathbf{x}\rVert = 1}\lVert T\mathbf{x}\rVert < \infty$ (the induced norm of Section 11.3). For linear maps, bounded is equivalent to continuous. Every bounded $T$ has an adjoint $T^$ defined by $\langle T^\mathbf{y}, \mathbf{x}\rangle = \langle\mathbf{y}, T\mathbf{x}\rangle$, the analog of the conjugate transpose, and $T$ is self-adjoint if $T^* = T$. Differentiation is unbounded: on $L^2[0, 1]$, the unit vectors $\sqrt{2}\sin(k\pi x)$ have derivatives of norm $k\pi \to \infty$. Unbounded operators can only be defined on a dense subspace (their domain), and boundary conditions are part of that domain, so they are part of the operator.

Spectrum versus eigenvalues. The spectrum of $T$ is

$$\sigma(T) = {\lambda \in \mathbb{C} : T - \lambda I \text{ has no bounded inverse}}.$$

In finite dimensions, $T - \lambda I$ fails to be invertible exactly when it fails to be injective, so the spectrum is the set of eigenvalues. In infinite dimensions, $T - \lambda I$ can be injective but not surjective, or surjective onto a dense subspace only, so the spectrum can contain points that are not eigenvalues:

Operator	Eigenvalues	Spectrum	Lesson
matrix $A \in \mathbb{C}^{n \times n}$	$\lambda_1, \dots, \lambda_n$	the same	finite dimensions: spectrum = eigenvalues
multiplication $(M\mathbf{u})(x) = x,u(x)$ on $L^2[0, 1]$	none	$[0, 1]$ (continuous spectrum)	self-adjoint without a single eigenvector
right shift $S$ on $\ell^2$	none	closed unit disk	truncations (nilpotent matrices, spectrum ${0}$) misrepresent it
left shift $S^*$ on $\ell^2$	the open unit disk	closed unit disk	$(1, \lambda, \lambda^2, \dots)$ is an eigenvector for $\lvert\lambda\rvert < 1$
compact self-adjoint $K$	real $\lambda_k \to 0$	eigenvalues together with ${0}$	behaves like a symmetric matrix

The spectrum of a bounded operator is always a nonempty compact set, and $\max{\lvert\lambda\rvert : \lambda \in \sigma(T)} = \lim\lVert T^k\rVert^{1/k}$ (Gelfand's formula, which extends the spectral radius of Section 11.3).

Compact operators: matrices in disguise. $K$ is compact if it maps bounded sets into sets with compact closure; equivalently, in a Hilbert space, $K$ is a norm limit of finite-rank operators. Integral operators $(K\mathbf{u})(x) = \int k(x, y)u(y),dy$ with square-integrable kernels are compact. The spectral theorem for compact self-adjoint operators is the direct generalization of Section 10.1:

$$K\mathbf{u} = \sum_{k=1}^\infty\lambda_k\langle\mathbf{e}_k, \mathbf{u}\rangle\mathbf{e}_k, \qquad \lambda_k \in \mathbb{R},; \lambda_k \to 0,$$

with ${\mathbf{e}_k}$ orthonormal (and complete after adding a basis of $N(K)$). Truncating the sum is the infinite-dimensional Eckart-Young theorem (Section 11.2). Two consequences:

• Fredholm alternative: for $\lambda \ne 0$, $K\mathbf{u} - \lambda\mathbf{u} = \mathbf{f}$ has a unique solution for every $\mathbf{f}$ exactly when the homogeneous equation has only the zero solution. This is the finite-dimensional "injective iff surjective" rescued for compact perturbations of the identity.
• Differential operators through their inverses: $-\frac{d^2}{dx^2}$ with boundary conditions is unbounded, but its inverse is an integral operator with a Green's function kernel, which is compact and self-adjoint. Its eigenfunctions (sines, cosines, Legendre and Bessel functions in Sturm-Liouville theory) therefore form orthonormal bases, which explains why separation of variables works.

Discretization. Numerical methods replace $T$ by an $n \times n$ matrix $T_n$. For compact operators, the eigenvalues of good discretizations converge to the true ones, as the code shows. For non-compact operators, truncation can create spurious eigenvalues (spectral pollution) or miss continuous spectrum, so finite-section results must be interpreted with care.

import numpy as np

# Integral operator (Ku)(x) = int_0^1 min(x, y) u(y) dy on L^2[0, 1]: compact and self-adjoint.
# It inverts -u'' with u(0) = 0, u'(1) = 0; exact eigenvalues 1 / ((k - 1/2)^2 pi^2), k = 1, 2, ...
exact = 1 / ((np.arange(1, 6) - 0.5) ** 2 * np.pi**2)
for n in (10, 100, 1000):
    y = (np.arange(n) + 0.5) / n              # midpoint rule nodes, weight 1/n each
    Kn = np.minimum.outer(y, y) / n           # n x n symmetric matrix approximating K
    lam = np.sort(np.linalg.eigvalsh(Kn))[::-1][:5]
    print(n, np.round(lam, 5))
print("exact", np.round(exact, 5))
# 10 [0.40612 0.04587 0.01707 0.00916 0.00593]
# 100 [0.40529 0.04504 0.01622 0.00828 0.00501]
# 1000 [0.40528 0.04503 0.01621 0.00827 0.005  ]
# exact [0.40528 0.04503 0.01621 0.00827 0.005  ]    matrices converge to the operator; eigenvalues -> 0

# Bounded vs unbounded: differentiation on the unit-norm functions sin(k pi x) * sqrt(2)
xs = np.linspace(0, 1, 100001)
for k in (1, 10, 100):
    u = np.sqrt(2) * np.sin(k * np.pi * xs)
    du = np.sqrt(2) * k * np.pi * np.cos(k * np.pi * xs)
    print(k, round(np.sqrt(np.trapezoid(u**2, xs)), 4), round(np.sqrt(np.trapezoid(du**2, xs)), 1))
# 1 1.0 3.1
# 10 1.0 31.4
# 100 1.0 314.2     ||u|| = 1 but ||u'|| = k pi: d/dx is unbounded (not continuous) on L^2

# Spectrum without eigenvalues: truncations of the right shift on l^2 are nilpotent Jordan blocks
n = 6
S = np.eye(n, k=-1)                           # (x1, x2, ...) -> (0, x1, x2, ...)
print(np.linalg.eigvals(S).real, np.linalg.norm(S, 2))   # [0. 0. 0. 0. 0. 0.] 1.0: all eigenvalues 0, yet ||S|| = 1
# The true shift on l^2 has NO eigenvalues; its spectrum is the closed unit disk

Key takeaways

• Bounded operators are the continuous linear maps; differentiation is unbounded and needs a domain, which includes the boundary conditions.
• The spectrum is where $T - \lambda I$ has no bounded inverse; in infinite dimensions it can contain non-eigenvalues (continuous spectrum), even for self-adjoint operators.
• Compact operators are norm limits of finite-rank ones; compact self-adjoint operators have real eigenvalues tending to 0 and an orthonormal eigenbasis.
• Inverses of differential operators are compact integral operators, which is why Sturm-Liouville eigenfunctions form bases.
• Discretizations converge spectrally for compact operators, but truncations of non-compact operators can mislead.

🧪 Practice

1. Show that the left shift $S^*$ on $\ell^2$ has every $\lambda$ with $\lvert\lambda\rvert < 1$ as an eigenvalue, and that the right shift $S$ has no eigenvalues at all.
2. Verify that $u_k(x) = \sqrt{2}\sin((k - \tfrac12)\pi x)$ satisfies $\int_0^1\min(x, y)u_k(y),dy = \frac{u_k(x)}{(k - \frac12)^2\pi^2}$. Which boundary conditions does $u_k$ satisfy?
3. Prove that the multiplication operator $(M\mathbf{u})(x) = x,u(x)$ on $L^2[0, 1]$ is bounded, self-adjoint, and has no eigenvalues. Hint: $x,u(x) = \lambda u(x)$ forces $u = 0$ almost everywhere.
4. Interview: You discretize a PDE eigenvalue problem with finite differences and observe that the smallest eigenvalues converge as the grid is refined while the largest ones grow without bound. Is this a bug? Hint: the differential operator is unbounded, so its eigenvalues go to infinity; only the low end of the spectrum is resolved by a given grid, and the matrix norm grows like $1/h^2$, which also drives the condition number.

Linear Algebra in Differential Geometry

Why it matters. Surfaces, curved spacetime in general relativity, robot configuration spaces, and the parameter spaces of statistical models are manifolds: spaces that look like $\mathbb{R}^n$ near every point but may be curved globally. Differential geometry studies them by attaching a vector space to every point and doing linear algebra there. The derivative becomes a linear map, an inner product varying from point to point measures lengths and angles, and the exterior algebra of Section 13.1 provides the objects that can be integrated. Calculus on curved spaces is linear algebra applied one point at a time.

Analogy. The Earth is round, but a city map is flat and works perfectly within the city. A manifold is a surface covered by overlapping city maps (charts), and the tangent space at a point is the flat map that matches the surface best right there.

Tangent spaces. At each point $p$ of an $n$-dimensional manifold $M$ there is an $n$-dimensional vector space $T_pM$, the tangent space, consisting of the velocity vectors of curves through $p$. For a surface in $\mathbb{R}^3$ parametrized by $X(u, v)$, the tangent space is the column space of the Jacobian $J = [\partial_uX\ \partial_vX]$, a plane through $p$. Tangent spaces at different points are different vector spaces. Comparing vectors at different points requires extra structure (a connection), and the failure of that comparison to be path independent is curvature.

The derivative is a linear map. A smooth map $F : M \to N$ has at each point a differential $dF_p : T_pM \to T_{F(p)}N$, the best linear approximation of $F$ near $p$. In coordinates it is the Jacobian matrix. The chain rule $d(G \circ F)p = dG{F(p)} \circ dF_p$ is matrix multiplication (Section 7.2). The inverse function theorem says: if $dF_p$ is invertible (a rank condition, Section 6.6), then $F$ is locally invertible. Rank conditions on the Jacobian also decide when a level set ${F = c}$ is a smooth manifold.

Covectors and the gradient. For a function $f : M \to \mathbb{R}$, the differential $df_p : T_pM \to \mathbb{R}$ is a linear functional, an element of the cotangent space $T_p^*M$ (the dual space of Section 7.4). It is naturally a covector, not a vector. Turning it into the gradient vector requires an inner product, $\nabla f = g^{-1},df$ (raising an index, Section 13.1). In a coordinate system that is not orthonormal, the gradient is not the vector of partial derivatives. This is why "steepest descent" depends on the metric, and why natural-gradient methods in machine learning use the Fisher information matrix as a metric.

Riemannian metrics. A Riemannian metric assigns to each point an inner product on $T_pM$: in coordinates, a symmetric positive definite matrix $g_{ij}(p)$ that varies smoothly. For a surface $X(u, v)$ in $\mathbb{R}^3$, the induced metric (the first fundamental form) is the Gram matrix

$$g = J^TJ = \begin{bmatrix} \partial_uX \cdot \partial_uX & \partial_uX \cdot \partial_vX \ \partial_vX \cdot \partial_uX & \partial_vX \cdot \partial_vX \end{bmatrix}.$$

It computes everything intrinsic to the surface:

Quantity	Formula	Linear algebra behind it
length of a tangent vector $\mathbf{v}$	$\sqrt{\mathbf{v}^Tg\mathbf{v}}$	induced norm (Section 9.1)
angle between tangent vectors	$\cos\theta = \mathbf{v}^Tg\mathbf{w}/(\lVert\mathbf{v}\rVert_g\lVert\mathbf{w}\rVert_g)$	inner product
length of a curve $\gamma(t)$	$\int\sqrt{\dot\gamma^Tg,\dot\gamma},dt$	norms of velocities
area element	$\sqrt{\det g},du,dv$	Gram determinant = squared area of $\partial_uX \wedge \partial_vX$ (Section 13.1)
gradient	$g^{-1},df$	raising an index

Differential forms. A $k$-form assigns to each point an alternating $k$-linear map on $T_pM$, an element of $\Lambda^k(T_p^*M)$. Forms are what integration is defined for: a $k$-form eats the $k$ tangent vectors spanning an infinitesimal parallelotope and returns its signed weight. The change of variables formula $\int f,d\mathbf{x} = \int f(\Phi(\mathbf{u}))\lvert\det D\Phi\rvert,d\mathbf{u}$ (Section 5.3) is the statement that $n$-forms scale by the determinant, which is the top exterior power (Section 13.1). The generalized Stokes theorem $\int_M d\omega = \int_{\partial M}\omega$ unifies the fundamental theorem of calculus with the theorems of Green, Gauss, and Stokes.

Curvature. The Riemann curvature tensor $R^i{}_{jkl}$ is a type-$(1, 3)$ tensor that measures how much a vector rotates when carried around a small loop. For surfaces it reduces to one number at each point, the Gaussian curvature $K = \det(\text{shape operator})$: the product of the principal curvatures, which are the eigenvalues of a symmetric matrix. Gauss's Theorema Egregium says $K$ depends only on $g$, so a sphere ($K > 0$) cannot be flattened onto a plane ($K = 0$) without distortion, which is why every world map distorts something.

import numpy as np

# Unit sphere parametrized by colatitude theta and longitude phi
def X(theta, phi):
    return np.array([np.sin(theta) * np.cos(phi), np.sin(theta) * np.sin(phi), np.cos(theta)])

def jacobian(theta, phi, h=1e-6):             # 3 x 2: columns span the tangent plane T_p M
    return np.column_stack([(X(theta + h, phi) - X(theta - h, phi)) / (2 * h),
                            (X(theta, phi + h) - X(theta, phi - h)) / (2 * h)])

th, ph = np.pi / 3, np.pi / 4
J = jacobian(th, ph)
p = X(th, ph)
print(np.round(J.T @ p, 8) + 0.0)             # [0. 0.]: tangent vectors are orthogonal to the radius
g = J.T @ J                                   # metric tensor (first fundamental form) in these coordinates
print(np.round(g, 6) + 0.0)                   # [[1. 0.] [0. 0.75]] = diag(1, sin^2 theta)

# Length of a tangent vector given in coordinates (d theta, d phi) = (0, 1): sqrt(v^T g v)
v = np.array([0.0, 1.0])
print(round(np.sqrt(v @ g @ v), 6), round(np.sin(th), 6))   # 0.866025 0.866025: a circle of latitude is shorter

# Area = integral of sqrt(det g) d theta d phi  (Section 5.3: Jacobian change of variables)
n = 400
thetas = (np.arange(n) + 0.5) * np.pi / n
area = sum(np.sqrt(np.linalg.det(jacobian(t, 0.0).T @ jacobian(t, 0.0))) for t in thetas) * (np.pi / n) * 2 * np.pi
print(round(area, 4), round(4 * np.pi, 4))    # 12.5664 12.5664

# The differential of f(x, y, z) = z at p is a covector; the gradient needs the metric to become a vector
df = np.array([0.0, 0.0, 1.0]) @ J            # components of df in (d theta, d phi)
grad = np.linalg.solve(g, df)                 # raise the index: grad^i = g^{ij} (df)_j
print(np.round(df, 6) + 0.0, np.round(J @ grad, 6))   # [-0.866025  0.] [-0.306186 -0.306186  0.75]
# df = (-sin theta, 0); J @ grad is the gradient as a 3D tangent vector
print(np.round(np.array([0, 0, 1.0]) - p * p[2], 6))   # [-0.306186 -0.306186  0.75]: e_z minus its normal part

Key takeaways

• Each point of a manifold has a tangent space $T_pM$, a vector space spanned locally by the columns of the Jacobian of a parametrization.
• The differential $dF_p$ is a linear map; the chain rule is matrix multiplication, and rank conditions control local invertibility.
• $df$ is a covector; the gradient $g^{-1}df$ requires a metric and is not just the vector of partial derivatives.
• A Riemannian metric is a smoothly varying SPD matrix $g = J^TJ$ (for embedded surfaces) that gives lengths, angles, and the area element $\sqrt{\det g}$.
• Differential forms are alternating multilinear maps; integration and change of variables are exterior algebra and determinants applied pointwise.

🧪 Practice

1. For polar coordinates $X(r, \theta) = (r\cos\theta, r\sin\theta)$, compute $J$, $g = J^TJ$, and $\sqrt{\det g}$. Recover the area element $r,dr,d\theta$.
2. Compute the metric of the cylinder $X(\theta, z) = (\cos\theta, \sin\theta, z)$ and show that it equals the identity. Explain why a cylinder can be unrolled flat without distortion while a sphere cannot.
3. For $f(x, y) = x$ on $\mathbb{R}^2$ with the metric $g = \operatorname{diag}(1, 4)$, compute $df$ and the gradient $g^{-1}df$. In which direction does $f$ increase fastest when lengths are measured with $g$?
4. Interview: In optimization over probability distributions, natural gradient descent replaces $\nabla_\theta L$ by $F(\theta)^{-1}\nabla_\theta L$, where $F$ is the Fisher information matrix. Explain this in the language of this topic, and why it makes the method invariant to reparametrization. Hint: $\nabla_\theta L$ is a covector $dL$; $F$ is a Riemannian metric on parameter space; raising the index with $F^{-1}$ gives a tangent vector whose meaning does not depend on the coordinates $\theta$.