Take the matrix $\begin{bmatrix}2 & 1\\ 1 & 2\end{bmatrix}$ and multiply it into the vector $(1, 0)$, then multiply the result by the same matrix again, twenty times in all. Two things happen. The vector swings around until it points along $(1, 1)$, and after that every multiplication makes it three times longer, to many decimal places, so by the end it is about 2.5 billion times longer than it started. Start from almost any other vector and the same thing happens.

That small experiment contains most of this section. The direction the vector settles into is an eigenvector, the factor of three is its eigenvalue, and the procedure itself is power iteration, the algorithm that estimates the spectral norm inside spectral normalization. Swap "twenty multiplications" for "twenty steps of a recurrent network" and you have the reason RNN gradients explode or vanish.

The rest of the book keeps asking questions of this kind. Why do gradients through fifty steps of an RNN die (§3.3)? Why can LoRA fine-tune a $4096 \times 4096$ matrix by training only 65,536 numbers (§14.5)? Why does the Muon optimizer throw away the singular values of its update (§8.6)? What exactly does PCA keep (§1.6)? Each one is a question about what a matrix does to directions. We will answer them with a handful of ideas (eigenvalues, the singular value decomposition, norms, and conditioning), building each from a $2 \times 2$ picture you can draw, to a proof, to NumPy, to the weight matrices of a trained Transformer.

A matrix is a map you can draw

A matrix $\mathbf{A} \in \R^{2 \times 2}$ sends each vector $\mathbf{x}$ to $\mathbf{A}\mathbf{x}$, and its columns tell you where the two basis vectors go. For our matrix, $\mathbf{e}_1 = (1, 0)$ lands on the first column $(2, 1)$ and $\mathbf{e}_2 = (0, 1)$ on the second column $(1, 2)$. Every other vector follows by linearity, $\mathbf{A}(x_1\mathbf{e}_1 + x_2\mathbf{e}_2) = x_1\mathbf{A}\mathbf{e}_1 + x_2\mathbf{A}\mathbf{e}_2$, so the whole map is fixed by two arrows.

A good way to see a matrix is to watch what it does to the unit circle. A linear map sends the circle to an ellipse (flattened to a line segment if the matrix is singular), and the unit square to a parallelogram whose area is $|\det\mathbf{A}|$. For our matrix the determinant is $2 \cdot 2 - 1 \cdot 1 = 3$, so areas triple. A negative determinant means the map also flips orientation, like a mirror.

The widget below lets you set the four entries and watch the circle deform. It also draws two special pairs of directions, which the rest of this section explains. The eigenvectors are the directions that the map only stretches, without turning. The singular vectors are the directions on the circle that land on the axes of the ellipse.

The shear leaves exactly one direction unturned, the $x$ axis, and the ellipse's axes point somewhere else entirely. The rotation turns every direction, so it has no real eigenvectors at all, yet its "ellipse" is the circle itself, with every direction an axis. Eigenvectors and the axes of the ellipse are different objects, and they agree only for special matrices. We take them one at a time.

Eigenvectors, the directions a matrix only stretches

An eigenvector of a square matrix $\mathbf{A} \in \R^{n \times n}$ is a nonzero vector $\mathbf{q}$ that $\mathbf{A}$ maps to a multiple of itself, and the multiple $\lambda$ is its eigenvalue,

In plain words, along an eigenvector the matrix acts like the single number $\lambda$. It stretches (if $|\lambda| > 1$), shrinks (if $|\lambda| < 1$), or flips (if $\lambda < 0$), and never turns.

How do we find them? Rewrite the equation as $(\mathbf{A} - \lambda\mathbf{I})\mathbf{q} = \mathbf{0}$. A nonzero solution exists exactly when $\mathbf{A} - \lambda\mathbf{I}$ is singular, that is, when $\det(\mathbf{A} - \lambda\mathbf{I}) = 0$. For a $2 \times 2$ matrix $\begin{bmatrix}a & b\\ c & d\end{bmatrix}$ this determinant is $(a - \lambda)(d - \lambda) - bc$, which expands to

where $\operatorname{tr}\mathbf{A} = a + d$ is the trace. Comparing with $(\lambda - \lambda_1)(\lambda - \lambda_2) = \lambda^2 - (\lambda_1 + \lambda_2)\lambda + \lambda_1\lambda_2$ shows that the eigenvalues add up to the trace and multiply to the determinant. Both facts hold for $n \times n$ matrices too.

A worked example. For $\mathbf{A} = \begin{bmatrix}2 & 1\\ 1 & 2\end{bmatrix}$ the trace is 4 and the determinant is 3, so $\lambda^2 - 4\lambda + 3 = (\lambda - 3)(\lambda - 1) = 0$ and the eigenvalues are 3 and 1. For $\lambda = 3$, $\mathbf{A} - 3\mathbf{I} = \begin{bmatrix}-1 & 1\\ 1 & -1\end{bmatrix}$, whose null space is spanned by $(1, 1)$, so $\mathbf{q}_1 = (1, 1)/\sqrt2$. For $\lambda = 1$, $\mathbf{A} - \mathbf{I}$ has every entry equal to 1, and $\mathbf{q}_2 = (1, -1)/\sqrt2$. Now the opening experiment makes sense. The starting vector $(1, 0)$ is $(\mathbf{q}_1 + \mathbf{q}_2)/\sqrt2$. After $k$ multiplications it is $(3^k\mathbf{q}_1 + 1^k\mathbf{q}_2)/\sqrt2$, so the $\mathbf{q}_1$ part wins by a factor of $3^k$, and the tangent of the angle to $\mathbf{q}_1$ is exactly $(1/3)^k$.

Three other $2 \times 2$ matrices show everything that can go wrong.

• A rotation by 90 degrees, $\begin{bmatrix}0 & -1\\ 1 & 0\end{bmatrix}$, has trace 0 and determinant 1, so $\lambda^2 + 1 = 0$ and $\lambda = \pm i$. No real direction survives a quarter turn unturned. The eigenvalues are complex, of absolute value 1, because a rotation neither stretches nor shrinks.
• A shear, $\begin{bmatrix}1 & 1\\ 0 & 1\end{bmatrix}$, has $(\lambda - 1)^2 = 0$, a double eigenvalue of 1, but $\mathbf{A} - \mathbf{I} = \begin{bmatrix}0 & 1\\ 0 & 0\end{bmatrix}$ has only a one-dimensional null space, the $x$ axis. There is one eigenvector where we needed two. Such a matrix is called defective.
• A triangular matrix, $\begin{bmatrix}3 & 0\\ 4 & 5\end{bmatrix}$, has eigenvalues 3 and 5 (for a triangular matrix they are the diagonal entries). Its eigenvectors are $(1, -2)/\sqrt5$ and $(0, 1)$, which meet at about 27 degrees instead of 90. Real eigenvalues, but skewed eigenvectors. We will meet this matrix again.

Diagonalization and powers

If $\mathbf{A}$ has $n$ linearly independent eigenvectors, collect them as the columns of $\mathbf{Q}$ and the eigenvalues in the diagonal matrix $\boldsymbol\Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_n)$. The $n$ equations $\mathbf{A}\mathbf{q}_i = \lambda_i\mathbf{q}_i$ say, column by column, $\mathbf{A}\mathbf{Q} = \mathbf{Q}\boldsymbol\Lambda$. Multiplying on the right by $\mathbf{Q}^{-1}$ gives the eigendecomposition

The second identity follows because the inner factors cancel in pairs, $\mathbf{A}^2 = \mathbf{Q}\boldsymbol\Lambda\mathbf{Q}^{-1}\mathbf{Q}\boldsymbol\Lambda\mathbf{Q}^{-1} = \mathbf{Q}\boldsymbol\Lambda^2\mathbf{Q}^{-1}$, and so on. In plain words, change to the eigenvector coordinates, where the matrix is a separate scalar on each axis, raise each scalar to the $k$-th power, and change back. Powers of a matrix are powers of its eigenvalues.

So the long-run behavior of $\mathbf{A}^k$ is set by the spectral radius $\rho(\mathbf{A}) = \max_i |\lambda_i|$. If $\rho(\mathbf{A}) < 1$, every component shrinks and $\mathbf{A}^k \to \mathbf{0}$. If $\rho(\mathbf{A}) > 1$, the component along the largest eigenvalue grows without bound. You have already used this. Gradient descent on a quadratic in §0.6 multiplies the error by $\mathbf{I} - \eta\mathbf{A}$ every step, whose eigenvalues are $1 - \eta\lambda_i$, and the stability condition $\eta < 2/\lambda_{\max}$ is exactly $\rho(\mathbf{I} - \eta\mathbf{A}) < 1$.

import numpy as np

def eig_sym_2x2(a, b, c):
    raise NotImplementedError

Symmetric matrices and the spectral theorem

The rotation, the shear, and the triangular matrix all misbehave in some way. Symmetric matrices, with $\mathbf{A} = \mathbf{A}^\top$, never do, and they are everywhere in machine learning. Covariance matrices, Hessians, Gram matrices $\mathbf{X}^\top\mathbf{X}$, and kernel matrices are all symmetric. Their good behavior is a theorem.

The spectral theorem. Every real symmetric matrix $\mathbf{A} \in \R^{n \times n}$ can be written as

with real eigenvalues $\lambda_i$ and orthonormal eigenvectors $\mathbf{q}_i$ (the columns of $\mathbf{Q}$). In plain words, a symmetric matrix is a rotation into a special set of perpendicular axes, a separate stretch along each axis, and a rotation back. Since $\mathbf{Q}$ is orthogonal, its inverse is its transpose, and the second form writes $\mathbf{A}$ as a weighted sum of $n$ rank-one matrices $\mathbf{q}_i\mathbf{q}_i^\top$, each of which projects onto one axis.

The Rayleigh quotient

How much does a symmetric matrix stretch a given direction? The Rayleigh quotient answers it,

Expand $\mathbf{x} = \sum_i c_i\mathbf{q}_i$ in the eigenvector basis. Orthonormality makes the cross terms vanish, so $\mathbf{x}^\top\mathbf{A}\mathbf{x} = \sum_i \lambda_i c_i^2$ and $\mathbf{x}^\top\mathbf{x} = \sum_i c_i^2$, and

In plain words, the Rayleigh quotient is a weighted average of the eigenvalues, with weights $c_i^2$ saying how much of $\mathbf{x}$ lies along each eigenvector. So it can never leave the range of the eigenvalues, and it reaches $\lambda_{\max}$ exactly at the top eigenvector. This gives the most useful characterization of the extreme eigenvalues, $\lambda_{\max} = \max_{\mathbf{x} \neq \mathbf{0}} R(\mathbf{x})$ and $\lambda_{\min} = \min_{\mathbf{x} \neq \mathbf{0}} R(\mathbf{x})$.

When $\mathbf{A}$ is a Hessian, $R(\mathbf{x})$ is the curvature of the loss along direction $\mathbf{x}$, and $\lambda_{\max}$ is the curvature of the steepest direction, often called the sharpness. §0.6 showed that gradient descent is stable only if $\eta < 2/\lambda_{\max}$. Cohen et al. (2021) found empirically that in full-batch training of neural networks, the sharpness rises until it reaches about $2/\eta$ and then hovers there, a regime they called the edge of stability.

Positive definite matrices

A symmetric matrix is positive definite (PD) if $\mathbf{x}^\top\mathbf{A}\mathbf{x} > 0$ for every $\mathbf{x} \neq \mathbf{0}$, and positive semidefinite (PSD) if $\mathbf{x}^\top\mathbf{A}\mathbf{x} \ge 0$. Geometrically, the quadratic $\mathbf{x}^\top\mathbf{A}\mathbf{x}$ is a bowl that curves upward in every direction, and its level set $\mathbf{x}^\top\mathbf{A}\mathbf{x} = 1$ is an ellipse with axes along the eigenvectors and half-lengths $1/\sqrt{\lambda_i}$.

Several conditions are equivalent, and each is useful in a different place.

• All eigenvalues are positive. By the Rayleigh quotient, $\mathbf{x}^\top\mathbf{A}\mathbf{x} \ge \lambda_{\min}\lVert\mathbf{x}\rVert^2$, which is positive for every nonzero $\mathbf{x}$ exactly when $\lambda_{\min} > 0$.
• A Cholesky factor exists, $\mathbf{A} = \mathbf{L}\mathbf{L}^\top$ with $\mathbf{L}$ lower triangular and a positive diagonal. You used it in §0.4 to sample from a multivariate Gaussian. Trying a Cholesky factorization is also the standard numerical test for positive definiteness, since it fails partway through if the matrix is not PD and costs a small fraction of an eigendecomposition.
• It is a Gram matrix of independent vectors, $\mathbf{A} = \mathbf{B}^\top\mathbf{B}$ with $\mathbf{B}$ of full column rank.

That last form shows where PSD matrices come from. For any $\mathbf{X} \in \R^{N \times d}$, $\mathbf{v}^\top\mathbf{X}^\top\mathbf{X}\mathbf{v} = \lVert\mathbf{X}\mathbf{v}\rVert^2 \ge 0$, so every Gram matrix is PSD, and it is PD exactly when no nonzero $\mathbf{v}$ has $\mathbf{X}\mathbf{v} = \mathbf{0}$. That fails whenever there are more features than examples ($d > N$), which is why the least-squares matrix $\mathbf{X}^\top\mathbf{X}$ can be singular. A covariance matrix is PSD for the same reason, since $\mathbf{v}^\top\boldsymbol\Sigma\mathbf{v} = \operatorname{Var}(\mathbf{v}^\top\mathbf{x}) \ge 0$, a variance. And adding $\lambda\mathbf{I}$ shifts every eigenvalue up by $\lambda$ without changing the eigenvectors, because $(\mathbf{A} + \lambda\mathbf{I})\mathbf{q}_i = (\lambda_i + \lambda)\mathbf{q}_i$. That is the one-line reason ridge regression in §0.4 always has a unique solution. Finally, at a strict local minimum the Hessian is PSD, and at a saddle point it has eigenvalues of both signs (§0.6).

The singular value decomposition

Eigenvectors need a square matrix, and even then they can be complex, missing, or skewed. The singular value decomposition (SVD) works for every matrix, of any shape, with no exceptions. It says that every matrix $\mathbf{M} \in \R^{m \times n}$ has an orthonormal basis $\mathbf{v}_1, \dots, \mathbf{v}_n$ of its input space and an orthonormal basis $\mathbf{u}_1, \dots, \mathbf{u}_m$ of its output space such that

where $k = \min(m, n)$, $\mathbf{U} \in \R^{m \times m}$ and $\mathbf{V} \in \R^{n \times n}$ are orthogonal, and $\boldsymbol\Sigma \in \R^{m \times n}$ holds the singular values $\sigma_i$ on its diagonal and zeros elsewhere. In plain words, every linear map is a rotation of the input ($\mathbf{V}^\top$), a stretch of each axis by $\sigma_i$ (which also drops or pads dimensions to change the shape), and a rotation into the output ($\mathbf{U}$). The $\mathbf{v}_i$ are the right singular vectors, the $\mathbf{u}_i$ the left singular vectors.

In two dimensions this is the ellipse in the widget. The right singular vectors are the two perpendicular directions on the unit circle that land on the ellipse's axes, the left singular vectors point along those axes, and the singular values are the semi-axis lengths. The number of nonzero singular values is the rank $r$. The first $r$ left singular vectors span the column space (everything $\mathbf{M}$ can output), and the last $n - r$ right singular vectors span the null space (everything it sends to zero). In code you almost always want the thin SVD, which keeps only $k$ columns, $\mathbf{U} \in \R^{m \times k}$, $\boldsymbol\Sigma \in \R^{k \times k}$, $\mathbf{V} \in \R^{n \times k}$.

A worked example you can check by hand. Take the triangular matrix $\mathbf{M} = \begin{bmatrix}3 & 0\\ 4 & 5\end{bmatrix}$ again. Then $\mathbf{M}^\top\mathbf{M} = \begin{bmatrix}25 & 20\\ 20 & 25\end{bmatrix}$, with eigenvalues $45$ and $5$ and eigenvectors $\mathbf{v}_1 = (1, 1)/\sqrt2$ and $\mathbf{v}_2 = (1, -1)/\sqrt2$. So $\sigma_1 = \sqrt{45} = 3\sqrt5 \approx 6.708$ and $\sigma_2 = \sqrt5 \approx 2.236$. The left singular vectors are $\mathbf{u}_1 = \mathbf{M}\mathbf{v}_1/\sigma_1 = (1, 3)/\sqrt{10}$ and $\mathbf{u}_2 = (3, -1)/\sqrt{10}$. As a check, $\sigma_1\sigma_2 = 15 = |\det\mathbf{M}|$, the area scaling.

Compare with the eigenvalues, 3 and 5. The matrix stretches some direction by 6.7, far more than its largest eigenvalue suggests, because its eigenvectors are skewed. For a symmetric matrix the two decompositions agree (the singular values are $|\lambda_i|$ and the singular vectors are eigenvectors, up to sign), but in general they answer different questions. Eigenvalues govern what happens after many repeated applications, $\mathbf{M}^k$. Singular values govern what one application can do. One inequality always holds, $\rho(\mathbf{M}) \le \sigma_1$, since $|\lambda|\,\lVert\mathbf{q}\rVert = \lVert\mathbf{M}\mathbf{q}\rVert \le \sigma_1\lVert\mathbf{q}\rVert$ for any eigenpair.

Matrix norms

How big is a matrix? There are several honest answers, and the three that matter most are all functions of the singular values.

• The spectral norm, $\lVert\mathbf{M}\rVert_2 = \max_{\mathbf{x} \neq \mathbf{0}} \lVert\mathbf{M}\mathbf{x}\rVert/\lVert\mathbf{x}\rVert = \sigma_1$, the largest stretch of any direction (also called the operator norm). To see that it equals $\sigma_1$, write $\lVert\mathbf{M}\mathbf{x}\rVert^2/\lVert\mathbf{x}\rVert^2 = \mathbf{x}^\top\mathbf{M}^\top\mathbf{M}\mathbf{x}/\mathbf{x}^\top\mathbf{x}$, a Rayleigh quotient of $\mathbf{M}^\top\mathbf{M}$, whose maximum is its top eigenvalue $\sigma_1^2$.
• The Frobenius norm, $\lVert\mathbf{M}\rVert_F = \sqrt{\sum_{ij} M_{ij}^2} = \sqrt{\sum_i \sigma_i^2}$, the Euclidean length of all entries. The second equality holds because $\sum_{ij} M_{ij}^2 = \operatorname{tr}(\mathbf{M}^\top\mathbf{M})$, and the trace of a symmetric matrix is the sum of its eigenvalues, here the $\sigma_i^2$.
• The nuclear norm, $\lVert\mathbf{M}\rVert_* = \sum_i \sigma_i$, the sum of all stretches.

In plain words, the spectral norm is the worst case over directions, the Frobenius norm is a root-mean-square over directions times a constant, and the nuclear norm is the total. All three are unchanged by rotations on either side, $\lVert\mathbf{Q}_1\mathbf{M}\mathbf{Q}_2\rVert = \lVert\mathbf{M}\rVert$ for orthogonal $\mathbf{Q}_1, \mathbf{Q}_2$, because rotations do not change singular values.

For a matrix of rank $r$ they are tied together by $\lVert\mathbf{M}\rVert_2 \le \lVert\mathbf{M}\rVert_F \le \sqrt r\,\lVert\mathbf{M}\rVert_2$ and $\lVert\mathbf{M}\rVert_F \le \lVert\mathbf{M}\rVert_* \le \sqrt r\,\lVert\mathbf{M}\rVert_F$ (the upper bounds are Cauchy-Schwarz on the $r$ nonzero singular values). The ratio

is the stable rank, a count of "how many directions really matter" that, unlike the rank, does not jump when tiny singular values appear or vanish. The spectral norm is also submultiplicative, $\lVert\mathbf{A}\mathbf{B}\rVert_2 \le \lVert\mathbf{A}\rVert_2\lVert\mathbf{B}\rVert_2$, since stretching by at most $\sigma_1(\mathbf{B})$ and then at most $\sigma_1(\mathbf{A})$ stretches by at most their product.

The spectral norm deserves one more sentence, because it bounds how much a whole network can amplify a perturbation. With activations whose slopes are at most 1 (ReLU and tanh, for example), a stack of linear layers satisfies $\lVert f(\mathbf{x}) - f(\mathbf{y})\rVert \le \prod_l \lVert\mathbf{W}_l\rVert_2\,\lVert\mathbf{x} - \mathbf{y}\rVert$. Spectral normalization (Miyato et al., 2018) divides each weight matrix by an estimate of its spectral norm, so every layer is at most 1-Lipschitz, which stabilized GAN training. The σReparam fix for attention entropy collapse in §5.7 uses the same estimate.

Low-rank approximation and the Eckart-Young theorem

The SVD writes a matrix as a sum of rank-one pieces, $\mathbf{M} = \sum_i \sigma_i\mathbf{u}_i\mathbf{v}_i^\top$, ordered from the largest stretch to the smallest. The obvious way to compress it is to keep the first $r$ pieces,

the truncated SVD. It needs $r(m + n + 1)$ numbers instead of $mn$. The Eckart-Young theorem says that nothing of rank $r$ does better.

The Eckart-Young theorem (Eckart & Young, 1936, for the Frobenius norm, and Mirsky, 1960, for every rotation-invariant norm). For every matrix $\mathbf{B}$ of rank at most $r$,

and the truncated SVD $\mathbf{M}_r$ achieves both bounds. In plain words, the best rank-$r$ approximation keeps the $r$ largest stretches, and its error is exactly what was thrown away, the next singular value in the spectral norm and the root of the sum of the remaining squared singular values in the Frobenius norm.

The equality part needs no work. $\mathbf{M} - \mathbf{M}_r = \sum_{i > r}\sigma_i\mathbf{u}_i\mathbf{v}_i^\top$ is itself written as an SVD, with singular values $\sigma_{r+1}, \sigma_{r+2}, \dots$, so its spectral norm is $\sigma_{r+1}$ and its Frobenius norm is $\sqrt{\sum_{i>r}\sigma_i^2}$. The lower bounds are the real content.

The fraction of the squared Frobenius norm that the rank-$r$ approximation keeps, $\sum_{i \le r}\sigma_i^2 / \sum_i\sigma_i^2$, is often called the energy it captures, and one minus it is the squared relative error. So everything about how compressible a matrix is lives in how fast its singular values fall. The cell builds a $60 \times 40$ matrix whose singular values decay like $0.8^i$, checks that the truncated SVD's errors are exactly the predicted tails, and then tilts the kept subspace slightly in 2,000 random ways. Every tilt is worse.

import numpy as np

def low_rank(M, r):
    raise NotImplementedError

def low_rank_errors(s, r):
    raise NotImplementedError

How low rank is a real model's weight matrix?

Random matrices are not compressible. The singular values of an $m \times n$ Gaussian matrix spread fairly evenly over a band (the Marchenko-Pastur law describes the shape), so a small rank keeps only a small share of the energy. Trained weights are different, and how different depends on the matrix. The widget computes the full spectrum of every weight matrix of a real trained Transformer, the two-layer character-level GPT that Chapter 6 trains on Shakespeare's sonnets, and compares each with a random Gaussian matrix of the same shape and Frobenius norm.

The differences are large. The learned position embeddings (a $128 \times 64$ matrix) have a stable rank of 1.4 and 99% of their energy in eight directions. Adjacent positions have an average cosine similarity of 0.98, so the 128 vectors trace a smooth curve through a few dimensions, and position 0, whose vector is about twice as long as the rest, supplies much of the top direction. The query and key matrices of the first layer keep 71% and 83% of their energy at rank 8, against 38% for a random matrix, which suggests that this model's attention compares characters along only a few directions. The value and output matrices sit much closer to random. At real scale the same questions matter in two ways. Trained weight matrices of large models have heavy-tailed spectra rather than a few dominant directions (Martin & Mahoney, 2021), so they compress poorly by plain truncation. But the change that fine-tuning makes is often very low rank (Aghajanyan et al., 2021, could fine-tune RoBERTa to 90% of full performance on one task with only 200 parameters), which is the bet LoRA makes.

Power iteration

The opening experiment was an algorithm. Power iteration finds the top eigenvector of a matrix by multiplying a random vector by it over and over, normalizing each time so nothing overflows,

Why does it work, and how fast? Take $\mathbf{A}$ symmetric with eigenvalues ordered $|\lambda_1| > |\lambda_2| \ge \dots$ and write the start as $\mathbf{x}_0 = \sum_i c_i\mathbf{q}_i$. Normalization only rescales, so $\mathbf{x}_k$ points along

Every term but the first shrinks by at least $|\lambda_2/\lambda_1|$ per step, so the tangent of the angle between $\mathbf{x}_k$ and $\mathbf{q}_1$ is at most $|\lambda_2/\lambda_1|^k$ times its starting value. In plain words, each multiplication amplifies the top direction more than any other, and the error shrinks geometrically at a rate set by the gap between the two largest eigenvalues. The eigenvalue estimate converges even faster. The Rayleigh quotient is a weighted average of eigenvalues with weights proportional to $c_i^2\lambda_i^{2k}$, so its error is of order $(\lambda_2/\lambda_1)^{2k}$, the square of the vector's error. This is a general feature of Rayleigh quotients for symmetric matrices, a first-order error in the vector gives only a second-order error in the eigenvalue.

Three ways power iteration fails are worth knowing. If $|\lambda_1| = |\lambda_2|$ (eigenvalues $3$ and $-3$, say) the direction never settles. If the start has $c_1 = 0$ exactly, it converges to the wrong eigenvector, although a random start makes this a probability-zero event and rounding error rescues it anyway. And if the gap is small, it is slow. The number of steps to gain a factor of $e$ is $1/\ln(\lambda_1/\lambda_2)$, about $\lambda_1/(\lambda_1 - \lambda_2)$, so a gap of 1% needs about a hundred steps per factor of $e$.

The same loop gives the spectral norm of any matrix $\mathbf{M}$, by running it on $\mathbf{M}^\top\mathbf{M}$ without ever forming that product, alternating $\mathbf{u} \leftarrow \mathbf{M}\mathbf{v}/\lVert\mathbf{M}\mathbf{v}\rVert$ and $\mathbf{v} \leftarrow \mathbf{M}^\top\mathbf{u}/\lVert\mathbf{M}^\top\mathbf{u}\rVert$, with $\sigma_1 \approx \mathbf{u}^\top\mathbf{M}\mathbf{v}$. Spectral normalization runs a single step of this per training step and keeps $\mathbf{u}$ between steps, since the weights change little, and the σReparam fix of §5.7 is built on the same estimate. And since a Hessian-vector product costs about as much as a gradient (§0.3), power iteration with Hessian-vector products finds the sharpness of a network with billions of parameters without ever forming its Hessian.

To get the top $k$ eigenvectors at once, multiply a block of $k$ vectors and re-orthonormalize them with a QR factorization each step, which keeps them from all collapsing onto $\mathbf{q}_1$. This is orthogonal iteration (or subspace iteration), and it converges at the rate $|\lambda_{k+1}/\lambda_k|$. A final Rayleigh-Ritz step, an eigendecomposition of the small $k \times k$ matrix $\mathbf{X}^\top\mathbf{A}\mathbf{X}$, rotates the block into the best individual eigenvector estimates. Production solvers refine the same ideas. Lanczos methods (scipy.sparse.linalg.eigsh) reuse every vector the iteration produces, and randomized SVD (Halko, Martinsson & Tropp, 2011, implemented as torch.svd_lowrank) uses one random block multiplication plus a few power steps to get the top singular vectors of a huge matrix.

import numpy as np

def power_iteration(A, iters, rng):
    raise NotImplementedError

def top_k_eigh(A, k, iters, rng):
    raise NotImplementedError

Condition numbers

The condition number of an invertible matrix is the ratio of its largest to its smallest stretch,

and it is infinite for a singular matrix. For a symmetric positive definite matrix it is $\lambda_{\max}/\lambda_{\min}$, the same $\kappa$ that set the speed of gradient descent in §0.6. It has a second meaning, just as important, which is how much a matrix amplifies errors.

Solve $\mathbf{A}\mathbf{x} = \mathbf{b}$, but with the right-hand side slightly wrong, $\mathbf{b} + \delta\mathbf{b}$ (from measurement noise, or just rounding). The solution moves by $\delta\mathbf{x} = \mathbf{A}^{-1}\delta\mathbf{b}$, so $\lVert\delta\mathbf{x}\rVert \le \lVert\delta\mathbf{b}\rVert/\sigma_{\min}$. Meanwhile $\lVert\mathbf{b}\rVert = \lVert\mathbf{A}\mathbf{x}\rVert \le \sigma_{\max}\lVert\mathbf{x}\rVert$, so $1/\lVert\mathbf{x}\rVert \le \sigma_{\max}/\lVert\mathbf{b}\rVert$. Multiplying the two inequalities,

In plain words, a relative error in the data can come back $\kappa$ times larger in the answer. The bound is tight, reached when $\mathbf{b}$ lies along the top left singular vector and $\delta\mathbf{b}$ along the bottom one. A practical rule follows. You lose about $\log_{10}\kappa$ significant digits. Float64 carries about 16 digits, so a matrix with $\kappa = 10^{12}$ leaves you four. Float32 carries about 7 and bfloat16 under 3 (§0.9), so a solve with $\kappa = 1000$ in bfloat16 returns noise.

This is why you should never solve least squares by forming $\mathbf{X}^\top\mathbf{X}$ (the "normal equations") when accuracy matters. The singular values of $\mathbf{X}^\top\mathbf{X}$ are $\sigma_i^2$, so $\kappa(\mathbf{X}^\top\mathbf{X}) = \kappa(\mathbf{X})^2$, and squaring the condition number doubles the digits lost. np.linalg.lstsq works with $\mathbf{X}$ directly through the SVD. Ridge regularization helps for the same reason, $\kappa(\mathbf{X}^\top\mathbf{X} + \lambda\mathbf{I}) = (\sigma_1^2 + \lambda)/(\sigma_n^2 + \lambda)$, which a modest $\lambda$ can bring down by orders of magnitude.

The error grows with $\kappa$ and stays below the bound $\kappa\,\varepsilon_{\text{mach}}$, where $\varepsilon_{\text{mach}} \approx 2.2 \times 10^{-16}$ is the float64 rounding error. At $n = 12$ the Hilbert matrix has $\kappa \approx 10^{16}$, every digit is gone, and the "solution" is 10% wrong even though nothing went wrong in the code. In the polynomial fit, the normal equations lose seven orders of magnitude of accuracy that lstsq keeps.

Conditioning shows up throughout training. Badly scaled input features make a linear model's Hessian $\mathbf{X}^\top\mathbf{X}/N$ ill-conditioned, which is why standardizing features speeds up gradient descent. Adam's per-parameter scaling is a crude diagonal attempt to reduce the condition number. And Jordan, who introduced Muon, observed that the updates SGD with momentum and Adam produce for Transformer weight matrices have very high condition numbers, dominated by a few directions, which Muon fixes by setting every singular value of the update to 1 (§8.6).

Projections and least squares

Projecting a vector onto a line is the most basic linear-algebra operation in machine learning. For a unit vector $\mathbf{u}$, the closest point to $\mathbf{x}$ on the line through $\mathbf{u}$ is $(\mathbf{u}^\top\mathbf{x})\,\mathbf{u} = \mathbf{u}\mathbf{u}^\top\mathbf{x}$, and what remains, $(\mathbf{I} - \mathbf{u}\mathbf{u}^\top)\mathbf{x}$, is perpendicular to $\mathbf{u}$. For a subspace with an orthonormal basis given by the columns of $\mathbf{Q} \in \R^{n \times k}$, add up the pieces, $\mathbf{P} = \mathbf{Q}\mathbf{Q}^\top$.

If the subspace is spanned by the columns of a general matrix $\mathbf{A} \in \R^{n \times k}$ of full column rank, find the closest point $\mathbf{A}\mathbf{w}$ to $\mathbf{x}$ by requiring the residual to be perpendicular to every column, $\mathbf{A}^\top(\mathbf{x} - \mathbf{A}\mathbf{w}) = \mathbf{0}$. That gives $\mathbf{w} = (\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\mathbf{x}$ and

In plain words, projecting twice is the same as projecting once, and the projection is symmetric. Its eigenvalues are therefore only 0 and 1 (from $\mathbf{P}^2 = \mathbf{P}$, any eigenvalue satisfies $\lambda^2 = \lambda$), with 1 on the subspace and 0 on its orthogonal complement, and its trace is the dimension $k$. Look again at the equation for $\mathbf{w}$. It is the least-squares normal equation. So least squares is a projection, of the targets $\mathbf{y}$ onto the column space of the features, and the fitted values are $\mathbf{P}\mathbf{y}$.

The SVD gives the general version. The pseudoinverse $\mathbf{M}^+ = \mathbf{V}\boldsymbol\Sigma^+\mathbf{U}^\top$, where $\boldsymbol\Sigma^+$ inverts the nonzero singular values and leaves the zeros, returns the minimum-norm least-squares solution $\mathbf{w} = \mathbf{M}^+\mathbf{y}$ even when $\mathbf{M}$ is rank deficient, and $\mathbf{M}\mathbf{M}^+$ is the projection onto its column space.

The last two lines remove a single direction from a vector, which is a real technique. Removing a "gender direction" from word embeddings was an early debiasing method (Bolukbasi et al., 2016), and projecting one learned direction out of a model's internal activations is how researchers delete a concept the network represents, the idea behind activation steering. Every PCA plot is a projection onto two directions, and the output of an attention head always lies in the span of its value vectors.

Whitening

Whitening transforms data so that its features are uncorrelated with unit variance, its covariance becoming $\mathbf{I}$. Center the data $\mathbf{X} \in \R^{N \times d}$ (rows are examples) to $\mathbf{X}_c$, and let $\mathbf{C} = \mathbf{X}_c^\top\mathbf{X}_c/N = \mathbf{Q}\boldsymbol\Lambda\mathbf{Q}^\top$ be its sample covariance. In this subsection $\mathbf{W} \in \R^{d \times d}$ is a whitening matrix, applied to each example as $\mathbf{z} = \mathbf{W}\mathbf{x}_c$, so the whitened data are $\mathbf{Z} = \mathbf{X}_c\mathbf{W}^\top$ with covariance $\mathbf{W}\mathbf{C}\mathbf{W}^\top$. Any $\mathbf{W}$ with $\mathbf{W}\mathbf{C}\mathbf{W}^\top = \mathbf{I}$ whitens. Two choices are standard,

Both work, since for instance $\mathbf{W}_{\text{PCA}}\mathbf{C}\mathbf{W}_{\text{PCA}}^\top = \boldsymbol\Lambda^{-1/2}\mathbf{Q}^\top\mathbf{Q}\boldsymbol\Lambda\mathbf{Q}^\top\mathbf{Q}\boldsymbol\Lambda^{-1/2} = \mathbf{I}$. In plain words, PCA whitening rotates the data onto its principal axes and rescales each axis to unit variance. ZCA whitening does the same and then rotates back, so the features keep their original orientation.

In fact every whitening matrix is one of these times a rotation. If $\mathbf{W}_1$ and $\mathbf{W}_2$ both whiten, then $\mathbf{R} = \mathbf{W}_2\mathbf{W}_1^{-1}$ satisfies $\mathbf{R}\mathbf{R}^\top = \mathbf{W}_2\mathbf{W}_1^{-1}\mathbf{W}_1^{-\top}\mathbf{W}_2^\top = \mathbf{W}_2\mathbf{C}\mathbf{W}_2^\top = \mathbf{I}$, using $\mathbf{W}_1^{-1}\mathbf{W}_1^{-\top} = (\mathbf{W}_1^\top\mathbf{W}_1)^{-1} = \mathbf{C}$, which is the whitening condition for $\mathbf{W}_1$ rearranged. So the choice among whitenings is a choice of rotation, and ZCA is special. It is the whitening that moves the data the least.

Two practical notes. Divide by $\sqrt{\lambda_i + \epsilon}$ rather than $\sqrt{\lambda_i}$, because directions with nearly zero variance are mostly noise, and whitening would blow that noise up to unit size. And whitening is everywhere once you look for it. BatchNorm and LayerNorm (§0.7) are whitening with a diagonal $\mathbf{C}$, rescaling each feature without decorrelating them. The Mahalanobis distance of §0.4 is the ordinary distance after ZCA whitening. Whitening sentence embeddings before comparing them with cosine similarity is a cheap fix for their anisotropy (Su et al., 2021), as §1.9 shows. And Muon's update is a whitened gradient. If the gradient is $\mathbf{G} = \mathbf{U}\boldsymbol\Sigma\mathbf{V}^\top$ (thin SVD, full rank), then $(\mathbf{G}\mathbf{G}^\top)^{-1/2}\mathbf{G} = \mathbf{U}\boldsymbol\Sigma^{-1}\mathbf{U}^\top\mathbf{U}\boldsymbol\Sigma\mathbf{V}^\top = \mathbf{U}\mathbf{V}^\top$.

import numpy as np

def whiten(X, method="zca", eps=0.0):
    raise NotImplementedError

PCA is the same decomposition again

Principal component analysis, which §1.6 uses to draw embedding spaces, is the spectral theorem applied to $\mathbf{C}$ and the SVD applied to $\mathbf{X}_c$, and the two agree. If $\mathbf{X}_c = \mathbf{U}\boldsymbol\Sigma\mathbf{V}^\top$, then $\mathbf{C} = \mathbf{X}_c^\top\mathbf{X}_c/N = \mathbf{V}(\boldsymbol\Sigma^2/N)\mathbf{V}^\top$. So the principal directions are the right singular vectors of the centered data, and the variance along direction $i$ is $\sigma_i^2/N$. By the Rayleigh quotient, the top direction is the unit vector $\mathbf{v}$ that maximizes the variance $\mathbf{v}^\top\mathbf{C}\mathbf{v}$ of the projected data. By Eckart-Young, projecting onto the top $k$ directions gives the best rank-$k$ approximation of the data. Maximum variance and minimum reconstruction error are the same choice. The latent semantic analysis of §1.4 is the truncated SVD of a word co-occurrence matrix.

Where this shows up in modern models

Gradients through a recurrence

Backpropagating through $T - t$ steps of a recurrent network multiplies the gradient by $T - t$ Jacobians, each $\operatorname{diag}(\phi'(\cdot))\,\mathbf{W}_{hh}$ (§0.8). If the activation's slope is at most $\gamma$, submultiplicativity bounds the whole product by $(\gamma\,\lVert\mathbf{W}_{hh}\rVert_2)^{T-t}$, so gradients must vanish when $\gamma\sigma_1 < 1$. In the long run, $\lVert\mathbf{W}^k\rVert_2^{1/k} \to \rho(\mathbf{W})$ (Gelfand's formula), so the spectral radius decides whether a linear recurrence eventually explodes. But for matrices with skewed eigenvectors the short run can look very different. The matrix below has both eigenvalues equal to 0.5, yet it more than triples a vector in one step, and the norms of its powers stay above 1 for the first four steps before they decay.

The last line is the idea behind orthogonal initialization (Saxe et al., 2014). An orthogonal matrix has every singular value equal to 1, so its powers neither shrink nor grow any vector, and gradients pass through many linear steps unchanged. §3.3 carries out the full analysis for trained RNNs.

Attention

The attention scores of one head are $\mathbf{X}\mathbf{W}_Q\mathbf{W}_K^\top\mathbf{X}^\top/\sqrt{d_k}$, where $\mathbf{W}_Q, \mathbf{W}_K \in \R^{d \times d_k}$. The $d \times d$ matrix $\mathbf{W}_Q\mathbf{W}_K^\top$ that defines how tokens are compared has rank at most $d_k$, so each head compares tokens along only $d_k$ directions, and multi-head attention (§4.5) is a sum of $h$ low-rank bilinear forms instead of one full-rank one. The attention matrix itself is row-stochastic (every row is a softmax, so it is nonnegative and sums to 1). Its rows sum to 1, so $\mathbf{A}\mathbf{1} = \mathbf{1}$ and 1 is an eigenvalue, and no eigenvalue is larger in absolute value. Repeated averaging pulls every token's vector toward the same mixture, and Dong, Cordonnier & Loukas (2021) proved that a stack of pure attention layers, without residual connections or MLPs, converges to a rank-one output doubly exponentially fast in depth. The residual stream of §5.1 is part of what prevents it.

One more rank fact limits every language model. The logits for $N$ contexts over a vocabulary of $V$ tokens are $\mathbf{H}\mathbf{E}^\top$, with hidden states $\mathbf{H} \in \R^{N \times d}$ and output embeddings $\mathbf{E} \in \R^{V \times d}$, so the logit matrix has rank at most $d$ even though it is $N \times V$. Yang et al. (2018) called this the softmax bottleneck and argued that natural language's true log-probability matrix has much higher rank than $d$.

LoRA

LoRA (Hu et al., 2021) freezes a pretrained weight $\mathbf{W}_0 \in \R^{d_{\text{out}} \times d_{\text{in}}}$ and learns an update $\Delta\mathbf{W} = \mathbf{B}\mathbf{A}$ with $\mathbf{B} \in \R^{d_{\text{out}} \times r}$ and $\mathbf{A} \in \R^{r \times d_{\text{in}}}$, so the update has rank at most $r$ and costs $r(d_{\text{in}} + d_{\text{out}})$ parameters. For a $4096 \times 4096$ projection with $r = 8$ that is 65,536 instead of 16.8 million. Eckart-Young tells you exactly what you give up. If full fine-tuning would have made the update $\Delta\mathbf{W}^\star$, the best LoRA could hope for leaves an error of $\sqrt{\sum_{i > r}\sigma_i(\Delta\mathbf{W}^\star)^2}$. In the LoRA paper, ranks from 1 to 8 performed about equally well when adapting GPT-3 175B's query and value projections on the two tasks where they compared ranks, which suggests those updates have a few dominant directions. §14.5 measures the spectra of real fine-tuning updates.