03 — Linear Algebra

This chapter builds the linear algebra that machine learning runs on, starting from geometry. The destination is a single claim you can verify by hand: attention is two matrix multiplications with a softmax between them.

Each idea is given its geometric meaning before its formula: vector, dot product, norm, matrix product, rank, basis, eigenvector, and singular value decomposition (SVD).

These ideas decide concrete questions:

• Why cosine similarity ignores length (Vectors norms and what cosine similarity actually measures).
• Why principal component analysis (PCA) is an eigenvalue problem (Pca derived from the eigendecomposition).
• Why a rank-50 approximation of a 1,000-column matrix keeps 98.9% of the information (Svd low rank approximation and a worked compression).
• Why backpropagation (backprop) runs backwards (Matrix calculus derived not memorized).
• Why an output token costs more than an input token (Matrix multiplication cost and why batching is a matmul).

Each result is stated with the assumption it rests on, in a paragraph labelled Assumptions, marking where it stops being true.

Probability decides what your model should optimize. Linear algebra decides what it costs to run and whether the optimizer can get there. Every embedding is a vector, every layer is a matrix, every attention score is a dot product, and every training step is a walk down a surface whose shape is an eigenvalue spectrum.

What goes in and what comes out

Linear algebra has exactly one input type and one output type, and everything else in this chapter is a special case of that one sentence.

The input is a list of numbers. A vector is an ordered list of numbers — (3, -4, 0, 12) is a vector of four numbers. Geometrically it is an arrow from the origin to a point in space: four numbers means an arrow in four-dimensional space, where the numbers are how far you travel along each of the four perpendicular directions. R^n, read aloud as “R to the n,” is the name for the set of all such lists of n real numbers — R^2 is the flat plane, R^3 is the space you live in, and R^768 is where a typical text embedding lives (an embedding is a vector of numbers a model produces to stand in for a word, sentence, or image, arranged so that similar things get nearby vectors).

The output is another list of numbers. A matrix is a rectangular grid of numbers. Feeding a vector into a matrix produces another vector: an m × n matrix (read “m by n” — m rows, n columns) takes a vector of n numbers in and returns a vector of m numbers out. That is the whole contract. A neural network layer is a matrix, so “a vector of 768 numbers goes in, a vector of 3,072 numbers comes out” describes both a layer and a matrix, and they are the same object.

The chapter’s destination makes this concrete. A transformer is the architecture behind every current large language model, and its defining layer is attention.

Here is attention stated purely as shapes. A token is the chunk of text the model actually reads, roughly a word or word-piece, and d is the length of each token’s vector. In one attention layer, n token vectors go in stacked as one n × d matrix, and n vectors of the same shape come out — each output vector a blend of the input vectors, weighted by relevance.

Same shape in, same shape out. What changes is the content: each row leaves carrying information from the other rows. Section 4 shows the blending is nothing but two matrix products, derives its cost, and checks the arithmetic on three tokens.

Reading the symbols

Every symbol below is plain ASCII, and each is also defined again at the point it is first used. Refer back to this table when a symbol is unfamiliar.

1. Vectors, norms, and what cosine similarity actually measures

There are only two measurements you can make on a vector: how long it is, and which way it points. Cosine similarity is the second measurement with the first one deliberately thrown away.

Norms

A norm is a measurement of length: it takes a vector and returns a single non-negative number saying how big that vector is. Geometrically, the ordinary norm is the length of the arrow from the origin to the point. Algebraically there is a whole family of them, and which one you pick changes what “big” means.

Every genuine norm has to satisfy three properties:

1. It is zero only for the all-zeros vector.
2. Scaling the vector scales the length: ||c·x|| = |c|·||x||, so doubling a vector doubles its norm.
3. Going directly is never longer than going via a detour: ||x + y|| <= ||x|| + ||y||, the triangle inequality.

Those three properties are what make the comparisons below meaningful, and one entry in the list below fails them.

Here are the four norms that come up, all computed on the same vector x = (3, -4, 0, 12). Read each line as: name, notation, formula, the answer for this x, nickname. In the notation, ||x|| is read “the norm of x,” the subscript names which norm, and sum_i means “add up one term for each coordinate i.”

L0    ||x||_0   = count of nonzero entries    = 3     not a real norm; the target of sparsity
L1    ||x||_1   = sum_i |x_i|                 = 19    "taxicab"
L2    ||x||_2   = sqrt( sum_i x_i^2 )         = 13    Euclidean, sqrt(9+16+0+144)
Linf  ||x||_inf = max_i |x_i|                 = 12    worst coordinate

The L2 norm is the one from the Pythagorean theorem — the straight-line distance from the origin. sqrt(9 + 16 + 0 + 144) = sqrt(169) = 13. When someone says “the norm” with no qualifier, this is it.

The L1 norm is the distance a taxi drives on a grid of streets: you cannot cut the diagonal, so you add up each coordinate’s absolute value, 3 + 4 + 0 + 12 = 19.

Linf, read “L infinity,” reports only the single worst coordinate — here 12 — and ignores the rest entirely.

L0 counts how many entries are nonzero. That count is what “sparsity” means. It is not a norm at all: doubling x to (6, -8, 0, 24) leaves the count at 3, so property 2 above fails. That failure is exactly why sparsity is hard to optimize directly, and why L1 is used as its tractable stand-in (Norms and regularization geometry).

The ordering, and why the last step needs a trick

The numbers came out 12 <= 13 <= 19, and that ordering is not a coincidence of this example. It holds for every vector:

||x||_inf  <=  ||x||_2  <=  ||x||_1  <=  sqrt(p)·||x||_2
   12          13          19           2 × 13 = 26

The first two steps are easy to see. ||x||_inf <= ||x||_2 because the largest coordinate squared is one of the terms in the sum of squares, so the sum is at least that large. ||x||_2 <= ||x||_1 because squaring sum_i |x_i| produces all the x_i^2 terms plus cross terms that are all non-negative, so ||x||_1^2 >= ||x||_2^2.

The last step is the one that needs work, and it rests on the Cauchy-Schwarz inequality: a dot product can never exceed the product of the two lengths, |a · b| <= ||a||·||b||. (The dot product is defined in the next subsection; for now treat it as the sum of products of matching coordinates. Cauchy-Schwarz is that subsection’s cos(theta) <= 1 in disguise.)

The obstacle is that Cauchy-Schwarz talks about dot products while ||x||_1 is a sum of absolute values. The key step is to write the sum as a dot product with the all-ones vector:

||x||_1 = sum_i |x_i|  =  |x| · 1          |x| = the vector of absolute values
                                           1   = the all-ones vector, p of them

        <=  || |x| ||_2 · ||1||_2          Cauchy-Schwarz, now that it is a dot product
        =   ||x||_2 · sqrt(p)              || |x| ||_2 = ||x||_2, and ||1||_2 = sqrt(p)

here:      19  <=  13 × 2  =  26

p is the number of coordinates — four in this example, so sqrt(p) = 2 and the bound reads 19 <= 26.

Two things to take from that derivation. Without the first line the appeal to Cauchy-Schwarz has nothing to apply to; that rewrite is the entire content of the proof. And the sqrt(p) is not mysterious — it is the length of the all-ones vector (1, 1, 1, 1) and nothing else, which is why the bound grows with the dimension.

The gap between L1 and L2 grows with how spread out the mass is, which is exactly why L1 penalizes many-small-coefficients more heavily than L2 does (Norms and regularization geometry).

Assumptions. The chain of inequalities holds for every real vector with no conditions attached. The two ends are tight rather than loose: ||x||_inf = ||x||_2 = ||x||_1 exactly when the vector has at most one nonzero entry, and ||x||_1 = sqrt(p)·||x||_2 exactly when every coordinate has the same absolute value. Those two extremes are the concentrated and the spread-out cases, which is the sentence above stated precisely.

Dot products and the geometry

The dot product of two vectors is a single number that measures how much the two arrows point the same way. Take them tip to tail: if they point in the same direction the dot product is positive, if they are perpendicular it is exactly zero, and if they point in opposite directions it is negative. Two vectors whose dot product is zero are called orthogonal, which is the general-dimension word for “at right angles.”

Algebraically it is the sum of the products of matching coordinates. The block below states that, then states the geometric reading as a second equality, then divides through to get cosine similarity.

a · b = sum_i a_i·b_i = ||a||·||b||·cos(theta)

cosine_similarity(a, b) = (a · b) / ( ||a||·||b|| )   in [-1, 1]

Read the first line as “a dot b equals the sum over i of a_i times b_i, which also equals the length of a times the length of b times the cosine of the angle between them.” theta is that angle. Since the cosine of an angle runs from -1 at 180 degrees through 0 at 90 degrees to 1 at 0 degrees, the dot product carries the alignment and the two lengths multiplied together.

The second form is the first one solved for cos(theta): divide out both lengths and what remains is pure direction, a number between -1 and 1 called the cosine similarity.

Here are both computed on a = (1, 2, 3) and b = (2, 4, 6). The last line shows the two metrics disagreeing completely about these vectors.

a · b   = 1·2 + 2·4 + 3·6 = 2 + 8 + 18 = 28
||a||   = sqrt(1+4+9)   = sqrt(14) = 3.7417
||b||   = sqrt(4+16+36) = sqrt(56) = 7.4833
product = 3.7417 × 7.4833 = 28.000   cos = 28 / 28 = 1.000

a - b = (-1, -2, -3)
Euclidean distance ||a - b|| = sqrt(1 + 4 + 9) = 3.7417

Cosine says these vectors are identical; Euclidean distance says they are 3.74 apart. b is 2a — same direction, twice the magnitude. That is the entire content of the metric: cosine discards magnitude and keeps only direction. For text embeddings that is usually the invariance you want — a 200-word and a 600-word article on the same topic should not rank apart because one is longer, and pooled embeddings do grow in norm with length. It is also information permanently discarded: confidence, intensity, and specificity partly live in the magnitude.

This measurement recurs throughout the chapter. An attention score is a dot product between a query vector and a key vector, which is why relevance in a transformer is an angle rather than an exact match; Matrix multiplication cost and why batching is a matmul shows that computing every such score at once is a matrix product. The same quantity reappears as projected variance in Pca derived from the eigendecomposition and as a gradient contraction in Matrix calculus derived not memorized.

Assumptions. Cosine similarity divides by both lengths, so it is undefined for the zero vector — a vector of all zeros has no direction, and a pooled embedding of an empty string can produce one. The identity a · b = ||a||·||b||·cos(theta) is a definition of angle in dimensions above three rather than a theorem about them, and it is well posed precisely because Cauchy-Schwarz guarantees the ratio never leaves [-1, 1].

Why vector search normalizes first

Once vectors are scaled to length one, three different ways of ranking search results collapse into the same ranking — so a vector database may as well use the cheapest.

A vector is unit-normalized when it has been divided by its own length, so ||a|| = 1 and only its direction survives.

Now expand the squared distance between two unit vectors:

||a - b||^2 = (a - b)·(a - b)
            = a·a - 2·(a·b) + b·b        expanding the square
            = 1   - 2·(a·b) + 1          since a·a = ||a||^2 = 1, same for b
            = 2 - 2·(a·b)

The whole point is the last line. Squared distance is 2 minus twice the dot product, and on unit vectors the dot product is the cosine, since the denominator ||a||·||b|| equals 1. So distance goes down exactly when cosine goes up, with no wiggle room.

On unit-normalized vectors, ranking by cosine, by inner product, and by Euclidean distance produce the identical ordering. (“Inner product” is another name for the dot product.)

That is why every vector database normalizes at write time and uses raw inner product at query time. It is the cheapest of the three — no square roots, no per-query norms — and provably equivalent to the other two.

Assumptions. The equivalence needs both vectors to be unit length. If documents are normalized but queries are not, cosine and inner product still agree — the query’s length is a constant that scales every candidate identically — but if the documents are not normalized, raw inner product starts preferring long vectors over well-aligned ones, and the three rankings genuinely diverge. That failure mode has a name, maximum inner product search, and it needs different index structures than cosine search does.

Pooling dilution: what one token in 200 is worth

Dense retrieval turns a whole passage into one vector and searches by cosine, so a rare and decisive word gets averaged in with 199 ordinary ones — the failure mode described in Embeddings and why dense search misses err_4021. Here is the arithmetic behind it.

Set up a deliberately simple model of a passage.

• The passage has 200 tokens, each with its own vector.
• Those 200 vectors are near-orthonormal: each has length one, and each is roughly perpendicular to all the others. That is what high-dimensional random directions approximately are.
• The pooled passage embedding is the sum S of the 200 token vectors. (Real pooling takes the mean, which is S/200; cosine ignores scale, so the sum gives the same answer with less notation.)

One consequence to get straight before the swap. Perpendicular lengths add in quadrature, not linearly — the Pythagorean theorem in 200 dimensions — so ||S||^2 = 200 and ||S|| ≈ sqrt(200) ≈ 14.14, not 200.

Now swap exactly one token: remove token vector u, add a fresh token vector v.

S' = S - u + v                     ||S'|| ≈ sqrt(200)  (still 200 unit vectors)

S · S' = S·S  -  S·u  +  S·v
       = 200  -   1   +   0
       = 199

cos(S, S') = 199 / (sqrt(200) · sqrt(200)) = 199 / 200 = 0.995

Each of the three terms is read straight off the model:

• S·S = ||S||^2 = 200, from the quadrature line above.
• S·u = 1, because u is one of the 200 unit vectors making up S, so it contributes its own squared length and nothing from the others.
• S·v = 0, because a fresh direction is perpendicular to everything already there.

Replacing one token in 200 moves the cosine by 0.005 — and the margin between the passage you want and a merely topical one is often smaller than that.

Now apply it to a real case. The rare literal ERR_4021 fragments into four low-information tokens. That count is Tokens’s, where the BPE split ["ERR", "_", "40", "21"] is worked out and the consequence is stated qualitatively.

The 2% is this chapter’s own arithmetic on top of it, not a figure quoted from there: on the near-orthonormal model above, four tokens out of 200 own 4/200 = 2% of the pooled vector’s direction.

So the embedding does not “miss” the string. Geometry gave the string a vote proportional to its token count, and that vote is tiny. BM25 — “Best Match 25,” the standard keyword-scoring function behind lexical search engines — weights the same term by how rare it is, the exact opposite weighting. That is why hybrid retrieval is a correction rather than a hedge.

Assumptions. The 0.995 is the output of a model, not a measurement: real token embeddings are not orthonormal, they are anisotropic, meaning they cluster into a cone rather than spreading evenly over all directions, which makes typical pairwise dot products positive instead of zero. That makes the true cosine higher than 0.995, so the conclusion — one token barely moves a pooled vector — is if anything understated. What would break the argument is short passages: at 5 tokens rather than 200, swapping one moves the cosine by roughly 0.2, and dense retrieval over short fields behaves completely differently.

2. Matrices as linear maps

So far vectors have only been measured. Matrices are what act on them, and a matrix is better read as a geometric transformation than as a grid of numbers — which is also what makes it obvious why multiplying matrices in a different order gives a different answer.

An m × n matrix A is a linear map from R^n to R^m — a function that takes an n-number vector and returns an m-number vector, with the property that it maps straight lines to straight lines and leaves the origin fixed.

That is what linear means, stated precisely: A·(x + y) = A·x + A·y and A·(c·x) = c·(A·x). Grids stay grids; nothing bends.

How to read a matrix: look at where it sends the axes

A key fact about a matrix is that its columns are the images of the basis vectors.

A basis is a minimal set of directions from which every other vector can be built by scaling and adding. In R^2 the standard basis is e_1 = (1,0) and e_2 = (0,1), read “e one” and “e two” — the unit steps along each axis.

Multiply A by e_1 and everything except column 1 gets multiplied by zero, so A·e_1 is column 1. Everything else follows by linearity: any vector is a combination of basis vectors, and the map preserves combinations.

So to read what a matrix does, look at where it sends (1,0) and (0,1). Three examples, each written as a list of rows:

scale x by 2      [[2, 0],
                   [0, 1]]      (1,0) -> (2,0)     (0,1) -> (0,1)

rotate 90 deg     [[0, -1],
                   [1,  0]]     (1,0) -> (0,1)     (0,1) -> (-1,0)

shear             [[1, 1],
                   [0, 1]]      (1,0) -> (1,0)     (0,1) -> (1,1)

Careful with the reading direction: A·e_1 picks out the first column, so the image of (1,0) is read down the left-hand column, not across the top row.

The first matrix stretches the plane horizontally and leaves it alone vertically. The rotation is a quarter turn counterclockwise. The shear leaves (1,0) alone and slides (0,1) over to (1,1), tilting vertical lines while keeping horizontal ones.

Order matters, because composition does

Composition of maps is matrix multiplication. A·B means “do B first, then A,” and its columns are the images of the basis vectors under the combined operation. (“Matmul” is the standard abbreviation for matrix multiplication.)

Matmul is non-commutative because function composition is. Rotating then scaling is not scaling then rotating, and the gap is not small:

R = [[0,-1],[1,0]]   (rotate 90 deg)     S = [[2,0],[0,1]]   (scale x by 2)

S·R = [[0,-2],[1,0]]        applied to (1,0)  ->  (0, 1)
R·S = [[0,-1],[2,0]]        applied to (1,0)  ->  (0, 2)

Follow the point (1,0) through each order rather than trusting the products.

In S·R the rotation runs first, sending (1,0) to (0,1). The horizontal stretch that follows does nothing to a vector with no horizontal component, so it stays at (0,1).

In R·S the scaling runs first, sending (1,0) to (2,0). The rotation then carries that to (0,2).

Same two operations, different answers, and the difference is a factor of two rather than a rounding detail.

Two facts later sections lean on

Transposing a product reverses it: (AB)^T = B^T·A^T. A^T is read “A transpose” and means the matrix flipped so that row i becomes column i. That reversal is why backpropagation gradients come out transposed (Matrix calculus derived not memorized) — it is not a convention someone chose.

An orthogonal matrix preserves every length and every angle. A matrix Q is orthogonal when it is square and Q^T Q = I, where I is the identity matrix that leaves every vector unchanged. Length preservation is one line:

||Q·x||^2 = (Q·x)^T (Q·x) = x^T · Q^T Q · x = x^T · I · x = x^T x = ||x||^2

Read it left to right: a squared length is a vector dotted with itself; (Qx)^T = x^T Q^T by the transpose rule just above; Q^T Q collapses to the identity by definition; and what is left is the original squared length.

Rotations and reflections are exactly the orthogonal matrices. They are the only maps that change nothing measurable about your data, which is why they appear as the length-preserving factors in the SVD (Svd low rank approximation and a worked compression).

Assumptions. “Matrix” and “linear map” are the same thing only after you fix a basis for the input space and one for the output space; the same map has different matrices in different bases, which is precisely the freedom Svd low rank approximation and a worked compression exploits. The length-preservation argument needs Q square: a tall matrix with orthonormal columns also satisfies Q^T Q = I and also preserves lengths, but Q Q^T is not the identity and the map cannot be undone, because it embeds a small space inside a bigger one instead of rearranging one space.

3. Rank, span, null space — and what low rank MEANS for data

Every matrix can be audited: count the directions it can actually reach, and the directions it destroys. Cashed out on a table of data, that count is the single empirical fact behind PCA, embeddings, recommender systems, and LoRA.

Three definitions

Take each one geometrically first, then in notation.

The span of a set of vectors is everything you can build from them by scaling and adding. The span of one vector is a line. The span of two independent vectors is a plane.

The rank of a matrix is how many genuinely independent directions its output can reach. A 3 × 3 matrix of rank 2 squashes all of space onto a plane.

The null space is the set of inputs the matrix flattens to zero — the directions it destroys. Anything sent there cannot be recovered from the output.

Written down:

span(v_1..v_k)  all linear combinations of those vectors
rank(A)         dimension of the column space = number of independent directions A can reach
null(A)         { x : Ax = 0 }, the directions A destroys
rank-nullity    rank(A) + dim(null(A)) = n     (n = number of columns)

Two pieces of notation there. { x : Ax = 0 } is read “the set of all x such that A times x is zero,” and dim means the dimension of that set — the number of independent directions inside it.

The column space is the span of the matrix’s columns. By Matrices as linear maps the columns are the images of the basis vectors, so their span is exactly the set of outputs the matrix can produce.

The rank-nullity theorem says the two counts always add up to the number of input dimensions. Every input direction is either carried through to the output or crushed to zero, and none is both.

Rank and null space on a concrete matrix

Take A = [[1,2,3], [2,4,6], [1,1,1]] — three rows, three columns.

Row 2 is 2 × row 1, so it contributes no direction row 1 did not already contribute. The three rows span only two independent directions, so the matrix has rank 2.

To find what it destroys, solve Ax = 0. Row 2 is redundant, so only two equations are worth writing:

x1 + 2·x2 + 3·x3 = 0
x1 +   x2 +   x3 = 0      subtract:  x2 + 2·x3 = 0  ->  x2 = -2·x3
                          then       x1 = -x2 - x3 = 2·x3 - x3 = x3

null space = span{ (1, -2, 1) }        dim = 1

check:  A·(1,-2,1) = (1-4+3, 2-8+6, 1-2+1) = (0, 0, 0)
rank-nullity:  2 + 1 = 3 = number of columns

Walk the steps. Subtracting the second equation from the first kills x1 and leaves x2 + 2·x3 = 0, so x2 = -2·x3. Substituting back into the second equation gives x1 = -x2 - x3 = 2·x3 - x3 = x3.

Two equations, three unknowns, so one unknown is free. Set x3 = 1 and the other two follow: x2 = -2 and x1 = 1. That single direction (1, -2, 1) is the whole null space, and the check confirms A sends it to the origin.

Rank 2 plus nullity 1 equals the 3 columns, as the theorem requires.

The diagram below shows that split: the input space forks into the part that survives and the part that is destroyed.

flowchart LR
    IN["Input space<br/>R^n"] --> NS["Null space<br/>collapsed to 0<br/>dim = n - r"]
    IN --> RS["Row space<br/>dim = r"]
    RS -->|"A maps 1-to-1"| CS["Column space<br/>everything A can reach<br/>dim = r"]
    NS -->|"A maps to 0"| Z(("0"))
    CS --> OUT["Output space<br/>R^m"]

    style CS fill:#2d6a4f,color:#fff
    style NS fill:#9d0208,color:#fff

Reading it: the input space R^n splits in two.

One part is the null space, of dimension n - r, which the matrix collapses to 0 — that is the red box, and nothing that lands there comes back.

The rest is the row space, of dimension r. On that piece the matrix is one-to-one: distinct inputs give distinct outputs and nothing is lost. Its image is the column space, everything A can reach, which also has dimension r and sits inside the output space R^m — the green box.

On the example above, n = 3, r = 2, and n - r = 1: a plane’s worth of input survives and one line’s worth is destroyed.

Assumptions. Rank-nullity holds for every matrix over every field with no conditions. What does not carry over cleanly is numerical rank: in floating-point arithmetic a matrix is almost never exactly rank-deficient, so rank in practice means “how many singular values (Svd low rank approximation and a worked compression) are above a tolerance,” and the answer depends on the tolerance you choose. Row reduction as done above is also numerically fragile; a production rank check uses the singular value decomposition instead.

What low rank means for data

Rank stops being an algebra exercise the moment the matrix is a table of data.

An n × p data matrix of rank r means every row is a linear combination of just r independent patterns. You logged p columns, but the data only ever varies in r directions. The other p - r are algebraically determined and carry zero additional information.

Exact rank deficiency is the pathological case. It is what “perfect multicollinearity” means — one predictor being an exact linear combination of others. It makes X^T X singular, so it cannot be inverted, and the regression has no unique solution.

The standard way to cause it by accident is a full one-hot encoding plus an intercept. A one-hot encoding turns a category into one column per level, and those columns always sum to the all-ones intercept column, so the intercept is an exact combination of them. Dropping one level is the fix.

The interesting case is approximate low rank. Here the singular values (Svd low rank approximation and a worked compression) decay quickly rather than hitting zero, so p measured features carry r << p effective dimensions — << is read “much less than.”

That single empirical fact is what a surprising number of techniques are betting on. The table below lists five of them and, for each, the specific low-rank claim that has to be true for it to work.

Two words in the first row are used before they are defined. Take them provisionally: an eigenvalue is the factor by which a matrix stretches one of the special directions it does not turn at all, and a covariance matrix is the square table recording how each pair of features varies together. Both get their full treatment in Eigenvalues and eigenvectors. That row only needs “a few big numbers and a long tail of small ones.”

The rest of the vocabulary in the table:

• PCA is principal component analysis, derived in Pca derived from the eigendecomposition.
• A manifold is a low-dimensional surface curving through a high-dimensional space, the way a sheet of paper is a 2-D thing sitting in 3-D.
• Recsys is recommender systems, where the matrix has one row per user and one column per item.
• LoRA is Low-Rank Adaptation, a fine-tuning method derived in Svd low rank approximation and a worked compression.
• Distillation trains a small model to reproduce a large one’s outputs.
• Activations are the intermediate vectors a layer produces, as opposed to its weights.

Low rank is also what rescues vector search

The curse of dimensionality is the observation that in high dimensions all pairwise distances concentrate toward the same value. Taken at face value it says nearest-neighbor search in 1,024 dimensions should be meaningless — every neighbour is about as far away as every other.

Vector search at 1,024 dimensions plainly works, so something in that argument does not apply.

The resolution is low rank. The concentration arithmetic depends on the intrinsic dimension, not the ambient one. The ambient dimension is the length of the stored vector — 1,024. The intrinsic dimension is how many directions the data actually explores, and for real embeddings that is in the tens. The curse applies to the intrinsic count, which is small enough to be harmless.

Assumptions. Every row of that table is an empirical claim about particular data, not a theorem, and each one is falsifiable. Truly high-rank data exists — white noise has a flat spectrum by construction, and a low-rank approximation of it keeps nothing. The right question to ask of any of these methods is not “is the data low rank” but “how fast do the singular values decay,” which is a plot you can produce in one line.

4. Matrix multiplication cost, and why batching is a matmul

One cost formula is worth memorizing, and it pays for itself three ways: it derives that attention is a matrix multiplication, it explains why reordering parentheses can make the same computation five times cheaper, and it settles why an output token costs more than an input token.

The cost formula

One formula covers the cost of every matrix multiplication, and it comes straight from counting how many dot products the output requires.

For A of shape (m × k) and B of shape (k × n), AB costs m·n·k multiply-accumulates, or 2·m·n·k FLOPs. A multiply-accumulate (MAC) is one multiplication plus one addition, and a FLOP is one floating-point operation, hence the factor of two. The shape of the count follows from the definition of the product: the output has m·n entries, and each entry is a dot product of length k, costing k multiplies and k adds. That formula answers more interview questions than any other in this chapter.

Assumptions. The count is for the straightforward algorithm on dense matrices. Asymptotically faster algorithms exist — Strassen’s method and its descendants beat m·n·k for large square matrices — but they are numerically touchier and rarely used inside a neural network. The count also assumes no structure to exploit: sparsity, block-diagonality, or a low-rank factorization each change it, and the last of those is the point of Svd low rank approximation and a worked compression.

Why attention is a matmul

Attention is the mechanism that lets each token in a sequence pull information from the other tokens. The entire operation is two matrix multiplications with one row-wise normalization between them, and the reason is the dot product of Vectors norms and what cosine similarity actually measures: a matrix product is exactly a table of all pairwise dot products, so “score every token against every token” and “multiply two matrices” are the same instruction.

The chain has five links, and each one is checkable.

Link 1: the input is a matrix. A sequence of n tokens, each represented by a vector of d_model numbers, stacks into one matrix X of shape (n × d_model) — one row per token. Stacking is not a modeling choice, it is just writing the vectors down in a grid.

Link 2: the three projections are three matmuls. Each token is turned into three different vectors by three learned matrices: a query (what this token is looking for), a key (what this token advertises about itself), and a value (what this token hands over if selected). Written per token these are three matrix-vector products; written for the whole sequence at once they are three matrix-matrix products, Q = X·W_Q, K = X·W_K, V = X·W_V, of shape (n × d_model)(d_model × d_k) giving (n × d_k). That collapse from n separate products into one is the same batching argument spelled out at the end of this section. The role vocabulary is developed further in Attention derived as content based lookup.

Link 3: one score is a dot product. How much should token i attend to token j? Take q_i · k_j. By Vectors norms and what cosine similarity actually measures that number is ||q_i||·||k_j||·cos(theta), so it is large when the query and the key point the same way and negative when they oppose. Relevance has become an angle.

Link 4: all the scores at once are a matrix product. There are n^2 such pairs, and the definition of matrix multiplication says entry (i, j) of Q·K^T is row i of Q dotted with row j of K. That is precisely q_i · k_j. So the full n × n table of scores is the single expression Q·K^T — no loop, no special-casing, just the product of a matrix with another matrix’s transpose. This is the load-bearing step: the score table is a matmul because a matmul is a table of dot products by definition.

Link 5: the weighted sum of values is another matrix product. The raw scores are divided by sqrt(d_k) and passed through softmax, a function that turns a row of arbitrary numbers into positive weights that sum to 1 by exponentiating each and dividing by the total. The division by sqrt(d_k) is there because a dot product of two d_k-dimensional random vectors has a standard deviation that grows like sqrt(d_k), and without the correction the softmax would saturate into a hard maximum — the derivation is in ml/04 — why sqrt(d_k) exists. Call the resulting n × n matrix of weights A. Each output row is then the weighted average of the value rows using that row of weights, and a weighted sum of rows is again a matrix product: A·V.

Put together, that is softmax(Q·K^T / sqrt(d_k))·V — matmul, normalize each row, matmul.

The whole chain on three tokens

Here is the chain worked end to end on n = 3 tokens, with d_k = 4 so that sqrt(d_k) = 2 and the scaled scores stay whole numbers.

Skip the projection step and take Q, K, V as already given. Q and K are 3 × 4 (three tokens, d_k = 4); V is 3 × 2 (three tokens, value width 2, kept narrow so the output is easy to read).

Q = [[2,0,0,0],      K = [[2,0,0,0],      V = [[1,0],
     [1,1,0,0],           [0,2,0,0],           [0,1],
     [1,1,1,0]]           [0,0,2,0]]           [1,1]]

The three keys point along three different axes, so they are mutually orthogonal — the near-orthonormal picture from Vectors norms and what cosine similarity actually measures made exact. That is what keeps the dot products clean.

Now every score is one dot product, and the whole table is one product. Four stages follow: raw scores, scaled scores, softmax weights, output.

Q·K^T = [[4, 0, 0],        row i, column j  =  q_i · k_j
         [2, 2, 0],
         [2, 2, 2]]

/ sqrt(4) = [[2, 0, 0],
             [1, 1, 0],
             [1, 1, 1]]

softmax by row:  [[0.786986, 0.106507, 0.106507],     e^2 / (e^2 + 2)   = 0.786986
                  [0.422319, 0.422319, 0.155362],     e   / (2e + 1)    = 0.422319
                  [0.333333, 0.333333, 0.333333]]     all equal -> 1/3

A·V = [[0.893493, 0.213014],
       [0.577681, 0.577681],
       [0.666667, 0.666667]]

Check the first row by hand, one stage at a time.

Scores. q_1 = (2,0,0,0) and k_1 = (2,0,0,0) give 2·2 = 4. Then q_1 · k_2 = 0 and q_1 · k_3 = 0, because those keys lie along axes where q_1 has no component. Row 1 of Q·K^T is (4, 0, 0).

Scale. Divide by sqrt(d_k) = 2, giving (2, 0, 0).

Softmax. Exponentiate each: e^2 = 7.389056, e^0 = 1, e^0 = 1. The total is 7.389056 + 2 = 9.389056. Divide each by the total:

7.389056 / 9.389056 = 0.786986
1.000000 / 9.389056 = 0.106507      (twice)

check: 0.786986 + 0.106507 + 0.106507 = 1.000000

Output. Weight the three rows of V and add:

0.786986·(1,0) + 0.106507·(0,1) + 0.106507·(1,1)
  first  coordinate: 0.786986 + 0        + 0.106507 = 0.893493
  second coordinate: 0        + 0.106507 + 0.106507 = 0.213014

Token 1 got mostly value row 1, which is what a strongly matched query should produce.

Row 3 is the opposite extreme. Its scaled scores are all equal at 1, so softmax returns a flat 1/3 each, and the output (0.666667, 0.666667) is the plain average of the three value rows — ((1+0+1)/3, (0+1+1)/3).

Uniform attention is an average, and peaked attention is a lookup; everything in between is a matmul away.

The code below reproduces all four stages and asserts every number printed above. Note how short attention is — two @ operators and a normalization.

import numpy as np

Q_attn = np.array([[2.0, 0, 0, 0], [1, 1, 0, 0], [1, 1, 1, 0]])
K_attn = np.array([[2.0, 0, 0, 0], [0, 2, 0, 0], [0, 0, 2, 0]])
V_attn = np.array([[1.0, 0], [0, 1], [1, 1]])
d_k = 4


def attention(Q, K, V, d_k):
    """Two matmuls with a row-wise softmax between them. That is all attention is."""
    scores = Q @ K.T / np.sqrt(d_k)                    # matmul 1: every pairwise dot product
    e = np.exp(scores - scores.max(axis=1, keepdims=True))
    weights = e / e.sum(axis=1, keepdims=True)         # each row now sums to 1
    return weights @ V, weights                        # matmul 2: weighted sum of value rows


out, W_attn = attention(Q_attn, K_attn, V_attn, d_k)

assert np.allclose(Q_attn @ K_attn.T, [[4, 0, 0], [2, 2, 0], [2, 2, 2]])
assert np.allclose(W_attn.sum(axis=1), 1.0)                       # softmax rows are weights
assert np.allclose(W_attn[0], [0.786986, 0.106507, 0.106507], atol=1e-6)
assert np.allclose(W_attn[2], [1 / 3, 1 / 3, 1 / 3])              # equal scores -> plain average
assert np.allclose(out[0], [0.893493, 0.213014], atol=1e-6)
assert np.allclose(out[2], V_attn.mean(axis=0))                   # row 3 IS the average of V

Assumptions. The derivation assumes the score is a plain dot product, which is what “scaled dot-product attention” means and is not the only choice — additive attention scores with a small network instead, and it is not a matmul. It assumes every token may see every other token; a causal language model masks the upper triangle before the softmax, which changes the weights but not the shape of the computation. And it assumes a single head: real models run several attention heads in parallel over slices of the vector and concatenate the results, which is more matmuls of the same form rather than a different operation.

Associativity: same value, different cost

Matmul is associative in value but not in cost. Associativity means (A·B)·C and A·(B·C) produce the identical matrix; the cost formula above says they need not cost the same, because the intermediate results have different sizes. Take A (100×1000), B (1000×100), C (100×10):

(A·B)·C:   100·100·1000 = 10,000,000   then  100·10·100 =   100,000   total 10.1M
A·(B·C):   1000·10·100  =  1,000,000   then  100·10·1000 = 1,000,000  total  2.0M

The difference is entirely the size of the intermediate.

Grouping to the left builds a 100 × 100 intermediate — the product A·B — at a cost of ten million multiply-accumulates. Grouping to the right builds a 1000 × 10 intermediate for one million, and the second product costs another million.

Five times cheaper from reordering parentheses, with an identical result. That is the whole idea behind low-rank factorization (Svd low rank approximation and a worked compression) and behind linear attention.

Applying it to attention

softmax(QK^T)V has to be computed left to right, because the softmax sits between the two factors and there is nothing to reassociate around. Cost is 4·n^2·d: two matmuls, QK^T and then (scores)·V, at 2·n^2·d FLOPs each.

Remove the softmax and the expression becomes (QK^T)V, a plain triple product, so you may regroup it as Q(K^T V). Still two matmuls, but now 4·n·d^2.

Those two costs are not close. Put realistic numbers in:

n = 8,192 tokens, d = 128 per head

with softmax:     4·n^2·d = 4 × 67,108,864 × 128 = 3.44e10 FLOPs   quadratic in n
reassociated:     4·n·d^2 = 4 × 8,192 × 16,384   = 5.37e8  FLOPs   linear in n

ratio = n/d = 8192/128 = 64x

3.44e10 is scientific notation for 34.4 billion.

The first line is quadratic in the sequence length because the n × n score table has to be built. The second never forms it: K^T V is only d × d — 128 by 128 regardless of how long the sequence is. Every “linear attention” variant is that one algebraic move, and the softmax is precisely what you give up to make it.

Where the memory figure comes from

Two claims get stacked in discussions of attention cost, and they have different sources.

The quadratic scaling is Attention and why context costs what it does’s result. That section draws the n × n score table, counts the n(n+1)/2 filled cells under the causal mask, and concludes that doubling the context roughly quadruples the work. It puts no byte or megabyte figure on the score matrix at all — its only memory sentence is about the KV cache.

Converting the scaling into bytes is this chapter’s arithmetic. At n = 8,192 in bf16 (2 bytes per number), one head’s score matrix is:

8,192 × 8,192 × 2 bytes = 134,217,728 bytes = 134 MB per head
134 MB × 32 heads       = 4.3 GB if materialized

That 4.3 GB is the number FlashAttention exists to avoid. FlashAttention computes the same result in tiles that stay in fast on-chip memory, so the full score matrix is never written out to main memory at all.

Assumptions. Reassociation is legal only because matrix multiplication is associative, and a nonlinearity between two factors destroys that: softmax(QK^T)·V cannot be regrouped, full stop. Linear-attention variants therefore compute a different function, not a cheaper route to the same one, and the accuracy they give up is the price. The n/d speedup also assumes n > d; at short sequence lengths the reassociated form is the more expensive one.

Why batching is a matmul

Feeding one input through a layer is a matrix-vector product. Feeding many at once is a matrix-matrix product, and the reason that is worth doing is not that it saves arithmetic — it does not — but that it reads the weights from memory once instead of once per input.

One row through a d × d layer is a matrix-vector product. B rows is a (B × d)(d × d) matrix-matrix product.

Here is the asymmetry that matters. FLOPs scale linearly with B — twice the rows, twice the arithmetic. But the weight matrix is read from memory exactly once either way, because all B rows use the same weights.

FLOPs       = 2·B·d^2
bytes moved = 2·d^2       (the weights in bf16; activations are small by comparison)

arithmetic intensity = 2·B·d^2 / (2·d^2) = B FLOPs per byte

Three terms in that block.

bf16 is bfloat16, a 16-bit floating-point format. Each weight occupies 2 bytes, so a d × d weight matrix occupies 2·d^2 bytes.

Arithmetic intensity is how many FLOPs the hardware performs per byte it reads from memory. Dividing the two lines cancels the 2·d^2 completely and leaves exactly B. The batch size is the arithmetic intensity.

Machine balance is the same ratio measured for the hardware itself, and it gives you something to compare B against:

machine balance of a modern accelerator = 1e15 FLOP/s / 3.3e12 byte/s ≈ 300 FLOPs/byte

B <  300  ->  memory-bound: the GPU idles waiting for weights
B >  300  ->  compute-bound: the matmul units are actually busy

A GPU (graphics processing unit) that can do 1e15 FLOPs per second while reading 3.3e12 bytes per second needs roughly 300 FLOPs of work per byte to keep its arithmetic units fed. Whichever side of the ratio is smaller is the bottleneck.

Those constants are not universal. They are rounded from the published specification of one current datacenter part, an NVIDIA H100 SXM — roughly 1e15 dense bf16 tensor-core FLOP/s against 3.35e12 bytes/s of HBM3 bandwidth, whose ratio is 295.

Why an output token costs more than an input token

That single inequality derives the entire prefill/decode split of The kv cache the most important mechanism in this chapter.

Language models generate in two phases. Prefill reads the whole prompt at once. Decode then emits one token at a time, reusing cached keys and values from all previous positions.

Prefill processes thousands of tokens in one matmul. B is in the thousands, comfortably above 300, so it is compute-bound and fast.

Decode emits one token per sequence per step. B equals the number of concurrent sequences, usually well under 300, so it is memory-bound and slow — the GPU spends most of its time waiting for weights it will use for a single row.

Output tokens cost more than input tokens because decode runs on the wrong side of that inequality. Continuous batching exists to push B back across it, by packing unrelated requests into one matmul.

Assumptions. The bytes-moved line assumes the weights dominate the traffic and are read exactly once, which requires the working set to fit in on-chip memory and no re-reading across tiles; for very large batches the activations stop being negligible and the intensity flattens out. The 300 is a property of one machine and moves with every hardware generation, so the number to carry is the method, not the constant. And the whole argument assumes a dense matmul: sparsity, or quantization — storing weights in fewer bits than 16 — changes both sides of the ratio at once.

5. Eigenvalues and eigenvectors

Costs are set by matrix shapes; long-run behavior is set by the eigenvalues. The key question is which directions a matrix merely stretches without turning, and raising a matrix to a power is governed entirely by the answer.

Most vectors come out of a matrix pointing somewhere new. A few come out pointing exactly where they went in, only longer or shorter. Those are the eigenvectors, and the stretch factor for each is its eigenvalue:

A·v = lambda·v          v != 0

Read that as “A applied to v equals lambda times v,” where lambda is a plain number. The condition v != 0 excludes the trivial zero vector, which every matrix sends to zero and which would otherwise satisfy the equation for every lambda.

Geometrically: the eigenvectors are the directions the map only stretches, never turns.

Everything else about eigenvalues follows from applying that repeatedly. A·v = lambda·v, so A·(A·v) = A·(lambda·v) = lambda·(A·v) = lambda^2·v, and in general A^k v = lambda^k v. The largest |lambda| therefore dominates any repeated application, because raising the biggest number to a high power leaves the others behind.

Computing them by hand

Take A = [[2, 1], [1, 2]].

The starting move is to turn the definition into something solvable. A·v = lambda·v rearranges to (A - lambda·I)·v = 0, which says the matrix A - lambda·I sends the nonzero vector v to the origin. A matrix can only do that if it flattens space — and “flattens space” is exactly a determinant of zero. The determinant is the factor by which a matrix scales area or volume, so zero means the output has no area at all.

So the eigenvalues are the lambda making det(A - lambda·I) = 0.

A - lambda·I = [[2-lambda, 1], [1, 2-lambda]]

det = (2-lambda)(2-lambda) - 1·1 = (2-lambda)^2 - 1 = 0
     ->  (2-lambda)^2 = 1  ->  2-lambda = ±1  ->  lambda = 1, 3

With the eigenvalues in hand, put each one back and solve for the direction. Writing v = (x, y), the first row of (A - lambda·I)·v = 0 reads (2-lambda)x + y = 0:

lambda = 3:  (2-3)x + y = 0  ->  -x + y = 0  ->  y =  x  ->  v = (1, 1)/sqrt(2)
lambda = 1:  (2-1)x + y = 0  ->   x + y = 0  ->  y = -x  ->  v = (1,-1)/sqrt(2)

check:  A·(1,1) = (2+1, 1+2) = (3,3) = 3·(1,1)         ok
        A·(1,-1) = (2-1, 1-2) = (1,-1) = 1·(1,-1)      ok
trace = 2+2 = 4 = 1 + 3                   det = 4-1 = 3 = 1 × 3

Each eigenvector is divided by sqrt(2) only to give it length one — any scalar multiple of an eigenvector is also an eigenvector, so the scale is yours to choose.

Geometrically this matrix stretches the diagonal (1,1) direction by 3 and leaves the anti-diagonal (1,-1) direction completely alone.

The last two checks in the block are free and worth doing every time. The trace is the sum of the entries down the main diagonal. Trace is the sum of eigenvalues; determinant is their product. Both are one-line computations that catch an arithmetic slip immediately, and both work at any size.

The spectral theorem

A matrix is symmetric when it equals its own transpose, so entry (i,j) matches entry (j,i). [[2,1],[1,2]] above is symmetric.

The spectral theorem says a real symmetric matrix has real eigenvalues and an orthonormal eigenbasis — its eigenvectors are mutually perpendicular and can be scaled to length one. Check it on the example: (1,1)·(1,-1) = 1 - 1 = 0, perpendicular as promised.

That lets you write any real symmetric A as

A = Q·D·Q^T        Q orthogonal, its columns the eigenvectors
                   D diagonal, holding the eigenvalues

Read it right to left as three steps: Q^T rotates into the eigenvector frame, D stretches each axis by its eigenvalue, Q rotates back.

This theorem is why the next three sections work at all, because the matrices they care about are all symmetric:

• A covariance matrix records how each pair of features varies together (Pca derived from the eigendecomposition).
• A Hessian is the matrix of second derivatives of a loss (Positive definiteness the hessian and the condition number).
• A Gram matrix holds all pairwise dot products of a set of vectors.

Repeated application is governed by the spectrum

The spectrum of a matrix is just the collection of its eigenvalues.

A^k scales the top eigendirection by lambda_max^k. Every other direction is scaled by a smaller power, so after enough applications only the top direction survives in relative terms. Three things are the same fact:

• Power iteration finds the top eigenvector by multiplying a random vector by A over and over and renormalizing.
• PageRank is that algorithm run on the web’s link matrix.
• Vanishing and exploding gradients are that algorithm run on the recurrent Jacobian — the matrix of partial derivatives of one time step’s output with respect to the previous step’s. Eigenvalues below 1 shrink gradients to nothing over many steps; eigenvalues above 1 blow them up (Recurrence and the vanishing gradient derived).

Same mechanism, different quantity, and it is worth being exact about which. ml/04 §5 derives the threshold of 1 from a matrix norm, bounding the product by s = max|tanh'| · ||W||. The statement here is about the eigenvalues of the recurrent Jacobian.

The two agree for a symmetric matrix, and in general ||A|| >= max|lambda|. So the norm version is the safer bound and the eigenvalue version the sharper description: a non-symmetric Jacobian can have every eigenvalue below 1 and still grow a gradient for many steps before the asymptotic decay takes over.

Assumptions. Every square matrix has n eigenvalues over the complex numbers, but a real matrix need not have a single real eigenvector — the 90-degree rotation of Matrices as linear maps turns every direction, so no real direction is merely stretched. Writing A = Q·D·Q^T requires more than existence of eigenvalues: it requires n independent eigenvectors, which the spectral theorem guarantees for real symmetric matrices and which fails for defective matrices such as [[1,1],[0,1]], whose only eigenvector direction is (1,0) even though the eigenvalue 1 appears twice. Uniqueness has its own condition: an eigenvector is only ever determined up to scale, and when two eigenvalues are equal there is a whole plane of eigenvectors and no preferred direction inside it. The “top eigenvalue dominates” argument additionally needs |lambda_1| strictly greater than |lambda_2| and a starting vector with a nonzero component along v_1; a tie makes power iteration wander in a subspace instead of converging.

6. PCA, derived from the eigendecomposition

Principal component analysis (PCA) — the standard method for compressing many correlated columns into a few uncorrelated ones — is an eigenvalue problem, not a recipe. Most candidates can say “PCA finds directions of maximum variance.” Far fewer can show why those directions are eigenvectors. It is four lines, reproducible on a whiteboard.

Setting up the problem

Start geometrically. You have a cloud of points. You want the single straight line through the origin such that the shadows of the points on that line are as spread out as possible.

Casting a shadow onto a direction is called projection. For a unit-length direction w, the projection of a point x is just the dot product x · w — one number, saying how far along w the shadow falls.

Spread is measured by variance. So “find the direction of maximum variance” means “find the unit vector w that maximizes the variance of the projections.”

Now the notation. X is the data matrix, n × p: one row per data point, one column per feature. Center it first, meaning subtract each column’s mean so the cloud sits around the origin. Let C = (1/n)·X^T X be the covariance matrix, which is p × p.

Project every row onto w, giving a vector of scores z = Xw, one score per data point. Because X was centered, z has mean zero, so its variance is just the mean of its squares:

Var(z) = (1/n)·z^T z              variance = mean of squares, since mean(z) = 0
       = (1/n)·(Xw)^T (Xw)        substitute z = Xw
       = (1/n)·w^T X^T X w        using (Xw)^T = w^T X^T from §2
       = w^T·C·w                  collecting X^T X / n back into C

So the variance of the shadows is w^T·C·w, read “w transpose C w.” Check the shapes: (1 × p)(p × p)(p × 1) gives 1 × 1, a single number. An expression of that form is called a quadratic form.

Maximizing it

The problem is now maximize w^T C w subject to w^T w = 1.

The constraint w^T w = 1 says w has length one. Without it the problem has no answer: doubling w quadruples w^T C w, so you could inflate the variance forever by lengthening w rather than by pointing it anywhere interesting.

Use a Lagrange multiplier, the standard device for maximizing under a constraint. Build a new function that adds the constraint back in, multiplied by an unknown scalar, then set its derivative to zero.

L(w, lambda) = w^T·C·w - lambda·(w^T·w - 1)

dL/dw = 2·C·w - 2·lambda·w = 0     ->     C·w = lambda·w

Two derivatives are used there. The derivative of w^T·C·w is 2·C·w because C is symmetric — that is result 2 in Matrix calculus derived not memorized, derived from scratch. The derivative of lambda·w^T·w is 2·lambda·w, the same result with C replaced by the identity.

Setting the difference to zero and dividing by 2 gives C·w = lambda·w. That is the eigenvector equation of Eigenvalues and eigenvectors, arrived at with no extra steps and nothing assumed.

The stationary points of “maximize projected variance” are exactly the eigenvectors of the covariance matrix.

Now ask what variance each of those directions actually achieves. Substitute the equation you just derived back into the objective:

w^T·C·w = w^T·(lambda·w)        substitute C·w = lambda·w
        = lambda·(w^T·w)        lambda is a scalar, pull it out front
        = lambda                since w^T·w = 1

The variance captured by a component is its eigenvalue. Not proportional to it — equal to it.

The rest follows. The maximum is the eigenvector with the largest eigenvalue. The second component solves the same problem restricted to the orthogonal complement — the set of directions perpendicular to the first — which gives the second eigenvector, and so on down. That restriction is where the requirement that components be uncorrelated comes from.

Explained variance ratio is lambda_k / sum_j lambda_j. Since eigenvalues sum to the trace (Eigenvalues and eigenvectors), the denominator is just the total variance of the original features.

The same derivation, worked by hand

Here is the whole derivation on four points, chosen small enough that the eigenvalues come out as integers. The points are already centered — their coordinates sum to zero in each column:

(-3,-1), (-1,-3), (1,3), (3,1)          so X is 4 × 2:  n = 4 points, p = 2 features

Step 1: build X^T X. It is 2 × 2. The top-left entry is the sum of squared first coordinates. The off-diagonal is the sum of the products of the two coordinates. The bottom-right is the sum of squared second coordinates:

top-left:      (-3)^2 + (-1)^2 + 1^2 + 3^2         =  9 + 1 + 1 + 9  = 20
off-diagonal:  (-3)(-1) + (-1)(-3) + (1)(3) + (3)(1) = 3 + 3 + 3 + 3 = 12
bottom-right:  (-1)^2 + (-3)^2 + 3^2 + 1^2         =  1 + 9 + 9 + 1  = 20

X^T X = [[20, 12],
         [12, 20]]

Step 2: divide by n to get the covariance.

C = X^T X / 4 = [[5, 3], [3, 5]]              (population divisor; see below)

Step 3: read off the eigenvalues using the two free checks. Trace is 5 + 5 = 10 and determinant is 25 - 9 = 16. By Eigenvalues and eigenvectors the eigenvalues sum to the trace and multiply to the determinant, which is exactly the characteristic quadratic:

lambda^2 - 10·lambda + 16 = 0
lambda = (10 ± sqrt(100 - 64)) / 2 = (10 ± 6) / 2   ->   lambda = 8, 2

lambda = 8:  w1 = (1, 1)/sqrt(2)         lambda = 2:  w2 = (1, -1)/sqrt(2)

The eigenvectors are the same (1,1) and (1,-1) as in Eigenvalues and eigenvectors, because [[5,3],[3,5]] has the same symmetric shape as [[2,1],[1,2]].

Checking it against the actual spread

Do not take the algebra on trust. Project the four points onto each direction and measure the variance directly.

The projections are just dot products. For the first point onto w1: (-3,-1) · (1,1)/sqrt(2) = (-3 + -1)/sqrt(2) = -4/1.4142 = -2.828.

onto w1 = (1,1)/sqrt(2):   (-4, -4, 4, 4)/sqrt(2) = (-2.828, -2.828, 2.828, 2.828)
                            variance = (8 + 8 + 8 + 8)/4 = 8.0      = lambda_1  ok

onto w2 = (1,-1)/sqrt(2):  (-2, 2, -2, 2)/sqrt(2)  = (-1.414, 1.414, -1.414, 1.414)
                            variance = (2 + 2 + 2 + 2)/4 = 2.0      = lambda_2  ok

total variance = 8 + 2 = 10 = trace(C)      PC1 explains 8/10 = 80%

The projected scores have mean zero, so their variance is the mean of their squares — (-2.828)^2 = 8.0 for each of the four, giving 8.0. That is the eigenvalue, exactly as derived. “PC1” is the first principal component, the direction with the largest eigenvalue.

Two practical notes

Both are about assumptions rather than arithmetic, and both are common interview follow-ups.

Centering is mandatory, and scaling usually is. Without centering, the first component points at the mean of the cloud rather than at any structure inside it. Without standardizing, a feature measured in dollars dominates one measured in years for purely dimensional reasons (Scaling which models care and exactly why). Standardizing means dividing each column by its standard deviation so all columns become unitless.

The 1/n versus 1/(n-1) choice does not matter here. It is exactly the n-1 question from Maximum likelihood worked twice — that is math/02, the statistics chapter in this track, not the ml/02 this chapter cites elsewhere. Either divisor rescales every eigenvalue by the same constant, so the components and the explained-variance ratios are unaffected.

Assumptions. The derivation needs C symmetric, which it always is, and the spectral theorem of Eigenvalues and eigenvectors then guarantees a full orthonormal set of components exists. Uniqueness is a separate matter: each component is determined only up to a sign flip, so implementations disagree about which end of the axis is positive, and when two eigenvalues are equal the components inside that subspace are arbitrary — a perfectly circular cloud has no first principal component at all. The larger assumption is that variance is what you care about. PCA finds only linear structure, so data lying on a curved manifold gets flattened wrongly, and a low-variance direction can still be the one that predicts your label, which is why PCA is not feature selection.

7. SVD, low-rank approximation, and a worked compression

The singular value decomposition (SVD) says every linear map is a rotation, a stretch, and another rotation — and that one factorization compresses matrices, explains how PCA is actually implemented, and names what LoRA is betting on.

Every matrix — rectangular, singular, anything — factors as

X = U·S·V^T          U: n × r orthonormal (left singular vectors)
                     S: r × r diagonal, s_1 >= s_2 >= ... >= 0
                     V: p × r orthonormal (right singular vectors)

Read the pieces one at a time.

U and V have orthonormal columns: mutually perpendicular, each of length one. By Matrices as linear maps that makes them rotations, possibly with a reflection thrown in. They change no length and no angle.

S is diagonal: zero everywhere except the main diagonal, where it holds the singular values s_1 >= s_2 >= ..., each one non-negative and sorted largest first. A diagonal matrix does nothing but scale each axis independently.

r is the rank (Rank span null space and what low rank means for data), so the three factors are exactly as wide as the map has directions to work with.

The diagram below runs a single vector through the three factors in order. The coloured box is the only step that changes anything measurable.

flowchart LR
    A["any vector x"] --> B["V^T<br/>ROTATE<br/>into the right basis"]
    B --> C["S<br/>SCALE<br/>axis lengths by s_i"]
    C --> D["U<br/>ROTATE<br/>into the output basis"]
    D --> E["Xx"]

    style C fill:#2d6a4f,color:#fff

Follow any vector x through the diagram left to right. V^T rotates it into the right basis — the input frame in which the map happens to be axis-aligned. S then changes the axis lengths, multiplying coordinate i by s_i and nothing more. Finally U rotates the result into the output basis, and what emerges is Xx, the same answer the original matrix would have given. S is the only green box because it is the only step that changes anything measurable: by Matrices as linear maps the two rotations preserve every length and every angle, so all of the map’s content — its rank, its conditioning, what it amplifies and what it destroys — lives in that one diagonal matrix. Everything the rest of this section does is read off the green box.

Every linear map is a rotation, an axis-aligned scaling, and another rotation. There is nothing else a matrix can do.

SVD and PCA are the same computation

The relationship is algebraic, not analogical. Expand X^T X with X written in factored form:

X^T·X = (U S V^T)^T (U S V^T)
      = V·S^T·U^T · U·S·V^T        by (AB)^T = B^T A^T, applied twice
      = V·S·U^T·U·S·V^T            S is diagonal, so S^T = S
      = V·S·S·V^T                  since U^T·U = I (orthonormal columns)
      = V·S^2·V^T

so    C = X^T X / n = V·(S^2/n)·V^T

Compare the last line to the spectral theorem’s A = Q·D·Q^T from Eigenvalues and eigenvectors. It is the same shape: an orthogonal matrix, a diagonal matrix, the transpose of the orthogonal matrix. Matching them up term by term gives the whole result.

V holds the principal component directions, and the PCA eigenvalues are s_i^2 / n.

Check it against Pca derived from the eigendecomposition’s four points, where X^T X = [[20,12],[12,20]]:

eigenvalues of X^T X:  trace 40, det 400 - 144 = 256   ->  32 and 8
singular values:       s = sqrt(32) = 5.657,  sqrt(8) = 2.828
PCA eigenvalues:       s^2 / n  =  32/4 = 8   and   8/4 = 2

Those are exactly the 8 and 2 computed by hand in §6.

Which is why every real implementation runs SVD on X and never forms the covariance. Building X^T X squares the condition number (Positive definiteness the hessian and the condition number) — you can see it in the derivation above, where S becomes S^2, so a ratio of 1e8 between the largest and smallest singular value becomes 1e16.

float64 is double-precision floating point, which carries about 16 significant decimal digits. A condition number of 1e16 consumes all of them, and the answer has no correct digits left. It is the identical argument as ml/02 — why nobody computes the inverse.

Assumptions. The SVD is the rare decomposition with no preconditions: every real or complex matrix has one, square or not, singular or not. The singular values are unique and their order is forced. The singular vectors are not: u_i and v_i can both be sign-flipped together, and when a singular value is repeated, the corresponding vectors are only determined up to an arbitrary rotation inside that subspace. A zero singular value leaves its u_i unconstrained entirely, which is exactly the null space of Rank span null space and what low rank means for data showing up in the factorization.

Low-rank approximation

Keeping only the largest few singular values gives the best possible simplification of a matrix, and “best” here is a theorem, not an approximation.

The SVD lets you rewrite X as a sum instead of a product:

X = sum_{i=1..r} s_i · u_i · v_i^T

Each term u_i·v_i^T is an outer product — a column vector times a row vector, (n × 1)(1 × p), producing a full n × p matrix of rank one. So the sum builds X out of r rank-one layers, and because s_1 >= s_2 >= ..., the layers arrive heaviest first.

Truncate after k terms and you get X_k, a rank-k matrix. The Eckart-Young theorem says this is not merely a decent approximation but the best possible rank-k approximation, with error

||X - X_k||_F^2 = sum_{i>k} s_i^2        the singular values you threw away, squared

Two norms are named in that theorem. The Frobenius norm of a matrix is the square root of the sum of all its squared entries — the L2 norm of Vectors norms and what cosine similarity actually measures applied to the matrix flattened into one long vector. The spectral norm is the largest singular value, the most any vector’s length can be stretched. Eckart-Young holds in both.

Verifying the error formula

Use the same four points from Pca derived from the eigendecomposition, and take k = 1. The pieces are s_1 = sqrt(32) = 5.657, u_1 = (-0.5,-0.5,0.5,0.5), and v_1 = (1,1)/sqrt(2).

The first row of X_1 is s_1 · u_1[0] · v_1^T = 5.657 × (-0.5) × (0.7071, 0.7071) = (-2, -2). Doing that for all four rows:

X       rows:  (-3,-1), (-1,-3), ( 1,3), (3,1)
X_1     rows:  (-2,-2), (-2,-2), ( 2,2), (2,2)
X - X_1     :  (-1, 1), ( 1,-1), (-1,1), (1,-1)

each residual has squared length 1 + 1 = 2, and there are 4 of them
||X - X_1||_F^2 = 4 × 2 = 8 = s_2^2                    ok
||X||_F^2 = 9+1 + 1+9 + 1+9 + 9+1 = 40 = 32 + 8        ok

Geometrically, the rank-one approximation replaces each of the four points by its shadow on the first principal direction, so (-3,-1) becomes (-2,-2).

The total squared error came out 8, which is exactly s_2^2 — the one singular value that was thrown away. The theorem predicted that number before any of the arithmetic was done.

Assumptions. Eckart-Young holds for the Frobenius and spectral norms, and more generally for any norm unchanged by rotations, but it is not true for arbitrary norms — the best rank-k fit under an L1 criterion is a different and much harder problem. It also assumes you want the best approximation of the matrix as a whole; if some entries are missing, as in a recommender system, truncated SVD is no longer the answer and matrix completion methods — which fit a low-rank factorization to the observed entries only — take over.

A worked compression

Here is what the theorem buys, in floats and in multiply-accumulates, on a dense 1,000 × 1,000 matrix stored in factored form:

dense storage    1,000 × 1,000                   = 1,000,000 floats
rank-k storage   k·(1,000 + 1,000 + 1)           =     2,001·k floats

The 2,001 per retained direction is one column of U (1,000 numbers), one column of V (1,000 numbers), and one singular value.

Put k = 50 in, and find where the two lines cross:

k = 50   ->  2,001 × 50 = 100,050 floats  =  10.0% of the original
break-even at 2,001·k = 1e6  ->  k = 500, so any k below half rank saves storage

Storage is the smaller part of the payoff. The point is that the singular values decay, so a small k keeps almost everything.

The decay assumption, stated honestly

The arithmetic below assumes s_i is proportional to 1/i — the tenth singular value is a tenth of the first, the hundredth a hundredth.

That is a stylized fact, not a derived one, and it is worth naming where it comes from. It is the spectral analogue of Zipf’s law, the same rank-versus-magnitude power law seen in word frequencies, city sizes and citation counts, and it is the shape repeatedly reported for term-document matrices, user-item rating matrices, and natural-image patch matrices.

Nothing forces it. The diagnostic is one plot: put the spectrum on log-log axes, where s_i ∝ 1/i is a straight line of slope -1. The arithmetic below is only as good as that line.

The energy calculation

Energy here means the sum of squared singular values, which by the Eckart-Young formula is exactly the part of the matrix you retain. With s_i ∝ 1/i, s_i^2 ∝ 1/i^2, and the proportionality constant cancels in the ratio:

energy kept = sum_{i<=50} 1/i^2  /  sum_{i<=1000} 1/i^2
            = 1.62513 / 1.64393
            = 98.86%

The sums converge quickly because 1/i^2 falls off fast — 1 + 0.25 + 0.111 + 0.0625 + .... The first 50 terms already account for nearly all of the total, and the remaining 950 contribute under 2% between them.

10% of the storage, 98.9% of the energy.

The same factorization also speeds up inference

Multiplying a vector by the factored form gets the same discount, by the associativity trick of Matrix multiplication cost and why batching is a matmul. Group right to left:

U·(S·(V^T x))     V^T x  is (50 × 1000)(1000)  =   50,000 MACs -> a 50-vector
                  S·(..) is 50 diagonal scalings =      50 MACs -> a 50-vector
                  U·(..) is (1000 × 50)(50)      =  50,000 MACs -> a 1000-vector
                                                   ----------
                                                     100,050 MACs

dense 1000 × 1000 matrix-vector product          = 1,000,000 MACs

Ten times faster, from the same factorization that saved the memory. Right to left keeps every intermediate a 50-number vector. Grouping left to right would rebuild the full 1,000 × 1,000 matrix first and throw the entire saving away.

LoRA is the same arithmetic

LoRA is Low-Rank Adaptation, the standard cheap method for fine-tuning — further training of an already-trained model on a narrower task.

The mechanism: freeze the original weight matrix W, and learn an update dW = B·A where B is (d × r) and A is (r × d). The product B·A has rank at most r by construction, because it factors through an r-dimensional middle.

The parameter count is the whole pitch:

d = 4,096, r = 16

LoRA:  2·d·r = 2 × 4,096 × 16 =    131,072 parameters
full:  d^2   = 4,096 × 4,096  = 16,777,216 parameters

131,072 / 16,777,216 = 0.78%

The bet LoRA makes is precisely a rank claim: that the update needed to adapt a pretrained model lives in a few dozen directions, even though the weights themselves do not. That is why LoRA works well for style and task adaptation and poorly for teaching genuinely new capability.

Assumptions. The compression numbers assume the spectrum decays; s_i ∝ 1/i is an empirical regularity observed in many natural matrices, not a law, and a matrix with a flat spectrum gains nothing from truncation. The 10× speedup assumes the factors are multiplied right to left and that the dense baseline was actually dense. LoRA’s premise is likewise empirical and testable: if the rank-r update cannot express the change you need, training simply plateaus, and raising r is the diagnostic.

8. Positive definiteness, the Hessian, and the condition number

One property decides whether a surface curves upward in every direction, and it matters in three places: it catches corrupt correlation matrices, it explains why deep networks stall at saddles rather than minima, and it predicts how many gradient steps a problem needs.

A symmetric matrix M is positive definite (PD) when x^T·M·x > 0 for all x != 0. Equivalently — and this is the version you will actually use — when all its eigenvalues are strictly positive.

Geometrically, the quadratic form x^T·M·x traces out a bowl. Positive definiteness says the bowl opens upward in every direction, with no flat and no downward-curving direction anywhere.

Positive semi-definite (PSD) relaxes > 0 to >= 0, allowing zero eigenvalues. The bowl may have flat-bottomed valleys, but it still never curves downward.

The property decides something in three separate places. Take them in turn.

Covariance matrices are always PSD

The proof is one substitution. Pca derived from the eigendecomposition showed that w^T·C·w = Var(Xw), the variance of the data projected onto w. A variance is a mean of squares, so it cannot be negative. Therefore w^T·C·w >= 0 for every w, which is the definition of PSD.

That gives you a free correctness check on real data. A “correlation matrix” with a negative eigenvalue is a data bug, not a finding.

The usual cause is pairwise deletion of missing values: each cell of the matrix is computed from whichever rows happened to have both of its two columns present, so every cell comes from a different subsample. The result is not the covariance of any dataset, and nothing requires it to be PSD.

The Hessian classifies critical points

The Hessian is the matrix of all second partial derivatives of a function, so it describes the local curvature of a loss surface. A critical point is a place where the gradient vanishes and the surface is momentarily flat.

At a critical point, the Hessian’s eigenvalue signs say which kind of flat you are standing on:

• All positive (PD) — a local minimum. The surface curves up in every direction.
• All negative — a local maximum.
• Mixed signs — a saddle, curving up along some directions and down along others, like the middle of a horse’s saddle.

Now count how likely each case is in high dimensions. Suppose the d eigenvalue signs behaved like independent coin flips. Then the probability that all d come out positive is 2^-d:

d = 20      2^-20   ≈ 1e-6
d = 1,000   2^-1000 ≈ 1e-301
d = 1e6     zero for any practical purpose

“Stuck in a local minimum” is almost never what is happening in deep learning. Saddles and plateaus are. That is the correct answer to “why does training stall.”

It is also why momentum helps. A saddle has an escape direction — the one with the negative eigenvalue — and momentum accumulates velocity along it. Plain SGD, stochastic gradient descent, which updates weights using the gradient from one small batch at a time, just dithers on the near-flat coordinates instead.

Convexity is a PSD condition; conditioning is an eigenvalue ratio

A function is convex when the straight line between any two points on its graph lies above the graph — equivalently, when it has no downward curvature anywhere.

Stated in this section’s vocabulary: f is convex exactly when its Hessian is PSD everywhere. That is how you prove that log loss is convex and that MSE-over-a-sigmoid is not (Logistic regression); MSE is mean squared error.

Convexity tells you whether the bowl points up. The condition number tells you how lopsided it is:

kappa = lambda_max / lambda_min        (for a symmetric PD matrix)

kappa = 1      a perfect circular bowl
kappa large    a long narrow ravine

Conditioning matters because gradient descent’s error shrinks by a factor of (kappa-1)/(kappa+1) per step. At kappa = 1 that factor is 0; at large kappa it approaches 1, meaning almost no progress per step.

The block below prices three situations, each with the same two features and different preprocessing. Each row gives the covariance, its eigenvalues, kappa, and the resulting step count. The step count is ln(1e-3) / ln((kappa-1)/(kappa+1)) — how many steps to cut the error by 1,000x.

age (sd 15) and income (sd 45,000), uncorrelated:
    C = diag(225, 2.025e9)          kappa = 9,000,000     ~3.1e7 steps
after standardizing:
    C = I                           kappa = 1             1 step
standardized BUT correlated at 0.99:
    C = [[1, 0.99], [0.99, 1]]      eigen 1.99 and 0.01   kappa = 199
    factor 198/200 = 0.99 per step  ->  ln(1e-3)/ln(0.99) = 687 steps

Row by row.

Row 1, raw features. A standard deviation of 15 years gives a variance of 15^2 = 225. A standard deviation of 45,000 dollars gives 45,000^2 = 2.025e9. diag means those two numbers sit on the diagonal with zeros off it, because the features are uncorrelated. For a diagonal matrix the eigenvalues are the diagonal entries, so kappa = 2.025e9 / 225 = 9,000,000. The per-step factor is (9e6 - 1)/(9e6 + 1) = 0.99999978, so it takes roughly 31 million steps to gain three decimal places.

Row 2, standardized. Dividing each column by its own standard deviation makes both variances 1, so C = I, every eigenvalue is 1, and kappa = 1. The per-step factor is 0 and one step lands on the answer.

Row 3, standardized but correlated. Both variances are still 1, but the features now correlate at 0.99. For [[1, c], [c, 1]] the eigenvalues are 1 + c and 1 - c, giving 1.99 and 0.01 and kappa = 199. The bowl is circular no longer — it is a ravine — and 687 steps are needed.

Where those rows come from, and what they add

Rows 1 and 2 are Scaling which models care and exactly why’s example: the same age (sd 15) and income (sd 45,000), the same kappa = 9,000,000, the same kappa = 1 after standardizing. What is new here is that they are derived from the eigenvalues rather than asserted.

Gradient descent and why the surface shape decides everything is the general result those rows instantiate: the largest stable step is set by the sharpest curvature, the step count by the flattest, so the whole problem is their ratio. It runs that argument on a generic kappa = 100 bowl and contains no age/income example.

Row 3 is what both chapters leave implicit, and it is the point of putting the three rows side by side. Standardizing fixes scale-induced conditioning and does nothing at all about correlation-induced conditioning — row 3 is fully standardized and still needs 687 steps.

Whitening fixes both: rotate into the PCA basis and divide each direction by sqrt(lambda_i), so every direction ends up with variance 1 and kappa becomes 1. That is the real reason decorrelating features can speed up training that scaling alone did not.

Testing positive definiteness cheaply

You almost never want to compute eigenvalues just to answer “is this matrix PD.”

The cheap test is to attempt a Cholesky factorization, M = L·L^T, which writes M as a lower-triangular matrix times its own transpose. It exists if and only if M is PD, so the attempt is the test: the factorization simply fails partway through when the matrix is not PD.

It costs n^3/3 flops, half of an LU factorization — the general-purpose “lower-upper” factorization used to solve linear systems — and computes no eigenvalues at all.

Assumptions. Positive definiteness is defined for symmetric matrices; for a non-symmetric matrix the quadratic form only sees the symmetric part, so the test silently answers a different question. The Hessian classification requires the function to be twice continuously differentiable and the critical point to be non-degenerate: a single zero eigenvalue makes the second-derivative test inconclusive, and that is exactly the plateau case. The 2^-d saddle count is a heuristic, not a theorem — real Hessian eigenvalue signs are correlated, and random matrix theory gives a sharper statement — though the qualitative conclusion survives every refinement. The convergence-rate formula is a worst-case bound for gradient descent with a fixed optimal step size on a quadratic (or a strongly convex, smooth function); with momentum the exponent improves to roughly sqrt(kappa), and stochastic gradients change the story again. Cholesky’s guarantee is exact for PD matrices and delicate for merely PSD ones, where the factorization exists but is not unique and numerical noise can push a zero eigenvalue slightly negative.

9. Matrix calculus, derived not memorized

Every gradient formula in backpropagation can be reconstructed from the shapes of the matrices involved, so none of them has to be memorized — including the two-line result behind every classifier you have written.

One convention makes all of this mechanical: a gradient has the same shape as the thing you differentiate with respect to.

A gradient is the vector of partial derivatives of a single output number with respect to every input, and it points in the direction of steepest increase. So dL/dW, read “the gradient of the loss L with respect to the weights W,” is shaped like W.

That is not a mnemonic. It is what lets you reconstruct every backprop formula by shape analysis alone. (Backprop is short for backpropagation, the algorithm that computes these gradients layer by layer.)

Four results with scalar outputs

Each of these is derived rather than quoted. Read the block first, then the walkthrough underneath.

1)  f = a^T·x                       df/dx = a
    componentwise: f = sum_i a_i·x_i, so df/dx_k = a_k.

2)  f = x^T·A·x                     df/dx = (A + A^T)·x   = 2·A·x if A symmetric
    f = sum_i sum_j x_i·A_ij·x_j
    df/dx_k = sum_j A_kj·x_j + sum_i x_i·A_ik = (A·x)_k + (A^T·x)_k

3)  f = ||X·w - y||^2               df/dw = 2·X^T·(X·w - y)
    let r = X·w - y, f = r^T·r, df/dr = 2r, and dr/dw = X.
    Chain rule with the shape constraint: df/dw must be p x 1, and the only
    contraction of X (n x p) with 2r (n x 1) giving that is X^T·(2r).
    Setting it to zero gives the normal equations of ml/02 section 2.

4)  softmax + cross-entropy         dL/dz = p - y
    p_k = e^{z_k} / sum_j e^{z_j},  L = -sum_k y_k·log p_k
    dp_k/dz_i = p_k·(delta_ki - p_i)
    dL/dz_i = -sum_k (y_k / p_k)·p_k·(delta_ki - p_i)
            = -sum_k y_k·(delta_ki - p_i)
            = -y_i + p_i·sum_k y_k  =  p_i - y_i          since sum_k y_k = 1

Result 1 is the linear case. a^T·x is the dot product written as a matrix product, so each input x_k contributes exactly a_k and nothing else. The gradient is a itself.

Result 2 differentiates a quadratic form. The subtlety is that x_k appears twice in the double sum sum_i sum_j x_i·A_ij·x_j — once when i = k and once when j = k. Differentiating picks up both occurrences, which is where the two terms (A·x)_k and (A^T·x)_k come from. When A is symmetric they are equal, giving 2·A·x. That is the result Pca derived from the eigendecomposition used to turn the PCA objective into an eigenvector equation.

Result 3 is least squares. r = X·w - y is the residual vector of prediction errors. df/dw has to be p × 1, and there is only one way to combine X (which is n × p) with 2r (which is n × 1) to get that shape: X^T·(2r). The shape rule picked the formula. The normal equations referred to there are the closed-form solution of linear regression, derived in Linear regression the fit.

Result 4, with the skipped step filled in

Two pieces of vocabulary first. The logits z are the raw, unnormalized scores a classifier produces; softmax turns them into probabilities p. And delta_ki, the Kronecker delta, is 1 when k = i and 0 otherwise. It lets one formula cover both “differentiate a probability with respect to its own logit” and “with respect to a different logit.”

The line dp_k/dz_i = p_k·(delta_ki - p_i) is where people stall, because it is usually quoted. Here it is derived, using the quotient rule on p_k = e^{z_k} / D where D = sum_j e^{z_j}.

Note first that dD/dz_i = e^{z_i}, since only the j = i term in the sum involves z_i. Then split by cases.

Case k = i — differentiating a probability with respect to its own logit. Both numerator and denominator depend on z_i:

dp_i/dz_i = [ e^{z_i}·D  -  e^{z_i}·e^{z_i} ] / D^2        quotient rule
          = (e^{z_i}/D) - (e^{z_i}/D)^2
          = p_i - p_i^2  =  p_i·(1 - p_i)

Case k != i — the numerator e^{z_k} does not involve z_i, so only the denominator contributes:

dp_k/dz_i = [ 0·D  -  e^{z_k}·e^{z_i} ] / D^2
          = -(e^{z_k}/D)·(e^{z_i}/D)
          = -p_k·p_i

The Kronecker delta packs both cases into one line: p_k·(delta_ki - p_i) is p_k·(1 - p_i) when k = i and -p_k·p_i otherwise. That is exactly the two cases above.

The rest of result 4 is then bookkeeping. The p_k from the delta identity cancels the 1/p_k coming from differentiating log p_k, which is the cancellation that makes the whole thing collapse. y is the one-hot label vector, so sum_k y_k = 1 in the last step.

Result 4 is why the last line of every classifier you have written is literally probs - onehot. It is the same fact as Entropy cross entropy kl, where cross-entropy fell out of maximum likelihood. One derivation gives you the loss, another gives you its gradient, and neither needs to be memorized.

Assumptions. Result 2 assumes A does not itself depend on x. Result 4 assumes the labels form a proper distribution summing to 1, which one-hot labels satisfy exactly and softened labels satisfy as long as they still sum to 1; it also assumes the softmax and the cross-entropy are fused into one operation, and computing them separately both loses the cancellation and invites numerical overflow. The layout convention matters too: this chapter uses the convention where dL/dW matches the shape of W, and texts using the opposite convention state every formula transposed.

The linear layer, and the shape rule as a theorem

Here the shape rule stops being a convenience and starts doing the derivation, because for a linear layer there is exactly one way to combine the available matrices that produces the required output shape.

The forward pass is Y = X·W + b. Name every shape before going further:

X   (B × d_in)       the batch of inputs, one row per example
W   (d_in × d_out)   the weight matrix
b   (d_out)          the bias vector, added to every row
Y   (B × d_out)      the output

B is the batch size, d_in the input width, d_out the output width.

Backprop hands the layer one thing: the upstream gradient G = dL/dY, of shape (B × d_out). That is the gradient the layer above has already computed and passed down. From G, X, and W, the layer must produce three gradients.

Here is the trick. The available matrices are G (B × d_out), X (B × d_in), and W (d_in × d_out). Ask what has to come out, and the answer is forced:

dL/dW = X^T·G          (d_in × B)(B × d_out) = (d_in × d_out)     shape of W    ok
dL/dX = G·W^T          (B × d_out)(d_out × d_in) = (B × d_in)     shape of X    ok
dL/db = sum over the batch rows of G                              shape of b    ok

Those are the only contractions with the right shapes, which is why shape-matching is a derivation and not a trick.

Checking it by hand

Use numbers small enough to verify every entry. Take X = [[1, 2]] (batch of one, two inputs), W = [[1, 3], [2, 4]], target t = [4, 10], and loss 0.5·||Y - t||^2.

The leading 0.5 is there so the derivative of the squared error comes out as exactly Y - t, with no stray factor of two to carry around.

Y = X·W = [1·1 + 2·2,  1·3 + 2·4] = [5, 11]
G = dL/dY = Y - t = [5-4, 11-10] = [1, 1]

dL/dW = X^T·G = [[1],[2]] · [[1, 1]] = [[1·1, 1·1],   = [[1, 1],
                                        [2·1, 2·1]]      [2, 2]]
dL/dX = G·W^T = [1,1] · [[1,2],[3,4]] = [1·1 + 1·3, 1·2 + 1·4] = [4, 6]

Now verify two of those entries without using the formulas at all, by tracing what each number touched in the forward pass.

dL/dW_21. The weight W_21 (row 2, column 1) multiplies X_2 = 2 on its way into Y_1. So nudging W_21 by eps changes Y_1 by 2·eps, and the loss by G_1 · 2·eps. Therefore dL/dW_21 = 1 · 2 = 2, matching the matrix answer.

dL/dX_1. The input X_1 reaches Y_1 through W_11 = 1 and reaches Y_2 through W_12 = 3. Both paths contribute, so dL/dX_1 = G_1·1 + G_2·3 = 1 + 3 = 4, again matching.

That is the point of the exercise: every gradient entry can be read off the forward pass by asking which outputs a given number touched and with what weight, and the matrix formulas reproduce those answers exactly.

The code below runs the same computation, then does a finite-difference check: nudge one weight by 1e-6, see how much the loss actually moves, and confirm the ratio matches the gradient that was derived.

import numpy as np


def linear_backward(X, W, G):
    """Gradients for Y = X @ W. Shapes force every line."""
    return X.T @ G, G @ W.T, G.sum(axis=0)     # dW, dX, db


X = np.array([[1.0, 2.0]])
W = np.array([[1.0, 3.0], [2.0, 4.0]])
G = np.array([[1.0, 1.0]])                     # = Y - t for t = [4, 10]

dW, dX, db = linear_backward(X, W, G)
assert np.allclose(dW, [[1.0, 1.0], [2.0, 2.0]])
assert np.allclose(dX, [[4.0, 6.0]])

# finite-difference check on W[1, 0]
eps = 1e-6
t = np.array([[4.0, 10.0]])
loss = lambda Wm: 0.5 * ((X @ Wm - t) ** 2).sum()
Wp = W.copy()
Wp[1, 0] += eps
assert abs((loss(Wp) - loss(W)) / eps - dW[1, 0]) < 1e-4

In that block @ is Python’s matrix-multiplication operator and .T is the transpose, so X.T @ G is the first formula above transcribed directly.

Assumptions. Shape-matching identifies the right formula uniquely only when the shapes are distinguishable. If B, d_in, and d_out all happen to be equal, several contractions have the correct shape and only the componentwise derivation settles which is right — so treat the shape rule as a fast check that is provably sufficient in the general case and needs care in the square case.

Why backprop is reverse-mode

Why are gradients computed from the loss backwards rather than from the inputs forwards? The answer is entirely the associativity argument of Matrix multiplication cost and why batching is a matmul applied to a chain of Jacobians.

The chain rule over L layers is a product of Jacobians. A Jacobian is the matrix of partial derivatives of one layer’s outputs with respect to its inputs, so composing layers means multiplying their Jacobians.

Two facts collide to decide the direction.

Matrix products are associative in value but not in cost (Matrix multiplication cost and why batching is a matmul). You may bracket the chain however you like and get the same answer at different prices.

The loss is a single number. So one end of the product is a 1 × d row vector — there is only one output to differentiate — while every other factor is d × d.

Start from the scalar end and every intermediate stays a row vector, costing d^2 per layer. Start from the other end and there is no row vector to start from, so every step is a full matrix-by-matrix product at d^3.

The diagram contrasts the two. Green is the cheap total, red the expensive one.

flowchart LR
    subgraph REV["Reverse mode -- multiply right to left"]
        R1["scalar dL/dy<br/>1 x d"] --> R2["x J_L<br/>row vector times matrix<br/>d^2 work"] --> R3["x J_L-1<br/>d^2 work"] --> R4["total: L·d^2"]
    end
    subgraph FWD["Forward mode -- multiply left to right"]
        F1["J_1<br/>full d x d Jacobian"] --> F2["x J_2<br/>matrix times matrix<br/>d^3 work"] --> F3["total: L·d^3<br/>in ONE pass"]
        G1["same total, other accounting:<br/>one tangent VECTOR, L·d^2"] --> G2["run it once per input:<br/>d x L·d^2 = L·d^3"]
    end

    style R4 fill:#2d6a4f,color:#fff
    style F3 fill:#9d0208,color:#fff

The upper track is reverse mode, multiplying right to left. It starts from the scalar dL/dy as a 1 × d row. Each subsequent step is a row vector times a matrix, which costs d^2, so across L layers the total is L·d^2.

The lower track is forward mode, multiplying left to right. There is no scalar at that end to start from, so the running product is a full d × d matrix. It begins at the Jacobian J_1, each step is a matrix times a matrix at d^3, and across L layers that is L·d^3 — in a single pass. That is the accounting the subgraph’s top row shows, and it is the one to state in an interview.

The subgraph’s second row is the same total reached a different way, not an extra cost to multiply in. Forward mode is more usually implemented one tangent vector at a time: seed the derivative with respect to a single input, push a d-vector through the chain at d^2 per layer for L·d^2, then repeat once per input — d × L·d^2 = L·d^3.

These are alternatives, and the standard slip is to charge for both. “d^3 per step and one pass per input” would give L·d^4, over-counting by a factor of d. Each accounting is already complete on its own: a full-Jacobian sweep needs only one pass because it carries every input’s derivative at once, and a per-input sweep costs only d^2 per layer precisely because it carries one.

Put realistic numbers on the two totals:

d = 1,024, L = 48 layers

reverse mode:  L·d^2 = 48 × 1,048,576     = 5.0e7
forward mode:  L·d^3 = 48 × 1,073,741,824 = 5.2e10        ~1,000x more

The gap is exactly a factor of d, which at a width of 1,024 is about a thousand. At a real model’s width it is worse.

Backprop is reverse-mode because there is one scalar output and millions of inputs; if it were the other way round, forward mode would win. That is the whole reason, and it is also why Jacobian-vector products (forward mode) are the right tool for sensitivity analysis with few inputs, while vector-Jacobian products (reverse mode) are the right tool for training.

Assumptions. Reverse mode wins on arithmetic and pays for it in memory: it needs the forward-pass activations kept alive until the backward pass reaches them, which is why activation checkpointing — recomputing some activations on the backward pass instead of storing them — exists. The comparison also assumes one scalar output; with m outputs, reverse mode costs m passes and the crossover moves. And it assumes every layer is differentiable, which hard thresholds and discrete sampling steps are not.

10. Norms and regularization geometry

Why does the L1 penalty produce exactly-zero coefficients while the L2 penalty never does? A single geometric property — rotation invariance — settles it, and the same machinery read in the SVD basis shows what ridge regression is really doing.

Regularization means adding a penalty on the size of the coefficients to a fitting problem, to stop it from chasing noise. Regularization why l1 zeroes and l2 does not derives L1 sparsity twice — geometrically from the constraint region’s corners, algebraically from the subgradient, which is the stand-in for a derivative at the sharp corner where the absolute value function has none. Here is the linear-algebra layer underneath both, which makes the result feel inevitable rather than lucky.

L2 is the only p-norm invariant under rotation. For orthogonal Q, ||Qx||_2^2 = x^T Q^T Q x = x^T x = ||x||_2^2 — the same one-line argument as Matrices as linear maps, since Q^T Q = I. No other p has this property.

Watch it fail for L1 on the smallest possible example. Take x = (1, 0) and rotate it 45 degrees, which sends it to (0.7071, 0.7071):

||x||_1 = 1 + 0 = 1        ||Q·x||_1 = 0.7071 + 0.7071 = 1.4142     changed by 41%
||x||_2 = 1                ||Q·x||_2 = sqrt(0.5 + 0.5) = 1          unchanged

The L1 norm of the rotated vector is sqrt(2) = 1.4142, a 41% increase, purely because the vector stopped lying on an axis. The L2 norm did not move, because turning an arrow does not change its length.

That is the deep statement about sparsity. A rotation-invariant penalty cannot possibly prefer the coordinate axes, because it cannot tell the coordinate axes apart from any other orthonormal basis.

Ridge — the name for L2-penalized regression — therefore treats every direction identically and can never single out “coefficient 3 should be zero.” It has no way to express that preference.

L1 is not rotation-invariant, and that asymmetry is exactly the mechanism by which it privileges the original basis. Which is the right thing to want when your columns are age, zip, and n_prior_claims rather than arbitrary linear combinations of them.

It is the same argument as “MLPs are rotation-invariant and trees are not” (Why gbdts still beat neural nets on tabular data) — an MLP is a multi-layer perceptron, the plain stack of dense layers, and a tree splits on one original column at a time.

A method that cannot distinguish the original basis from a rotated one must learn which basis matters from data. A method that is basis-aware gets it for free.

Assumptions. The claim is about p-norms specifically: L2 is the only one of that family unchanged by every rotation. Rotation invariance is also a statement about the penalty, not about the whole fitting problem, and it produces sparsity only when the original coordinates are meaningful — after PCA, where the columns are already arbitrary linear combinations, an L1 penalty sparsifies a basis nobody chose and the interpretability argument evaporates.

Ridge in the SVD basis, which explains what it actually does

Rewriting both ordinary least squares and ridge in terms of the singular values turns “it penalizes large weights” into a precise statement about which directions survive.

Write X = U·S·V^T. Ordinary least squares (OLS) — plain linear regression with no penalty — and ridge become sums over singular directions:

OLS:    w = sum_i  ( u_i^T·y / s_i ) · v_i
ridge:  w = sum_i  [ s_i / (s_i^2 + lambda) ] · (u_i^T·y) · v_i

shrinkage of direction i, relative to OLS  =  s_i^2 / (s_i^2 + lambda)

Both lines build the coefficient vector one singular direction at a time, adding up p contributions.

u_i^T·y is how much of the target lies along direction i. OLS divides that by s_i, so a direction the data barely explores — a tiny s_i — gets divided by a tiny number and its coefficient explodes.

Ridge replaces 1/s_i with s_i/(s_i^2 + lambda). As s_i shrinks toward zero, that expression shrinks toward zero too instead of blowing up. Dividing the ridge factor by the OLS factor gives the shrinkage ratio s_i^2 / (s_i^2 + lambda).

Here lambda is the penalty strength, not an eigenvalue. The symbol is overloaded across the literature and this is the other meaning.

The table below evaluates that shrinkage ratio at lambda = 1 for four different directions, from well-explored to barely explored. Read down the first column as “how much variance the data has along this direction.”

Each row is one division: 100/101 = 0.990, 10/11 = 0.909, 1/2 = 0.500, 0.01/1.01 = 0.0099.

Ridge does not shrink coefficients uniformly. It shrinks each singular direction by how little variance the data has along it. High-variance directions survive untouched; directions the data barely explores are deleted.

That is a materially better answer than “it penalizes large weights,” and it explains three things at once.

Why ridge fixes collinearity. Collinearity is a tiny s_i. OLS’s 1/s_i blows up, while ridge’s s_i/(s_i^2 + lambda) is capped at 1/(2·sqrt(lambda)) — the maximum of that expression, reached at s_i = sqrt(lambda).

Why ridge needs standardized inputs. The s_i carry the units of the columns, so comparing them to a single lambda is only meaningful once the columns are unitless.

What the effective degrees of freedom are. This is a continuous count of how many parameters the model is really spending, and it is just the sum of the shrinkage factors: sum_i s_i^2/(s_i^2 + lambda). At lambda = 0 every factor is 1 and the sum is p; as lambda grows the sum falls toward 0.

Assumptions. The OLS line requires every s_i to be strictly positive, which is to say X must have full column rank; with an exactly zero singular value OLS has no unique solution at all, and that is Rank span null space and what low rank means for data’s perfect multicollinearity reappearing. Ridge has no such requirement — adding lambda to the denominator is exactly what makes the problem solvable for any X — which is the deeper reason it is the standard fix. Both lines assume the penalty is applied to all coefficients uniformly and that the intercept has been left out of it, and the whole reading assumes the columns were standardized first, since the singular values carry the units of the columns.

The code below checks the PCA-through-SVD identity and the shrinkage table against the numbers derived above:

import numpy as np


def pca_via_svd(X, k):
    """PCA without ever forming X^T X. Returns scores, components, eigenvalues."""
    Xc = X - X.mean(axis=0)
    U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
    # 1/n convention. math/02 section 2 covers 1/(n-1) -- math/02, not ml/02,
    # which this chapter also cites and which a bare "ch 02" would not separate.
    eigenvalues = s ** 2 / len(X)
    return U[:, :k] * s[:k], Vt[:k], eigenvalues[:k]


X = np.array([[-3.0, -1.0], [-1.0, -3.0], [1.0, 3.0], [3.0, 1.0]])
scores, comps, eig = pca_via_svd(X, k=2)
assert np.allclose(eig, [8.0, 2.0])                       # matches the hand derivation
assert np.allclose(scores.var(axis=0), [8.0, 2.0])        # variance IS the eigenvalue
assert np.allclose(np.abs(comps[0]), [2 ** -0.5] * 2)     # w1 = (1,1)/sqrt(2)


def ridge_shrinkage(s, lam):
    """Per-singular-direction shrinkage relative to OLS."""
    return s ** 2 / (s ** 2 + lam)


assert abs(ridge_shrinkage(np.sqrt(0.01), 1.0) - 0.00990) < 1e-5
assert abs(ridge_shrinkage(np.sqrt(100.0), 1.0) - 0.99010) < 1e-5

Cheat sheet

Every row below is derived somewhere above. This table is for recall, not first contact — if a row does not make sense, the derivation is in the section it came from.

Which results hold unconditionally, and which do not

This is the table interviewers probe. For each decomposition or theorem: when does the object exist at all, when is it unique, and what makes it fail. “Unique” matters more than it sounds — a non-unique answer means two correct implementations can disagree, which is why PCA components flip sign between libraries.

Next: SQL patterns.