Linear algebra

Nearly everything in machine learning is held as a list of numbers or a grid of them, and nearly every operation is one of those multiplied by another. This section covers what those objects represent, why the dot product gets used as a measure of similarity, and the rules about shapes that decide whether a multiplication is allowed at all. It is also the practical cure for the error you will hit most often, which is two things whose dimensions do not line up.

A vector is a list of numbers with a meaning#

A vector is an ordered list of numbers, written [2, 1] or [0.31, -0.08, 1.24, ...]. Ordered matters: the first number always means the same thing as the first number of every other vector you compare it with. That is the whole definition. Everything else said about vectors is a claim about what those numbers are for.

In machine learning a vector is usually a representation: an embedding of a word, the pixels of an image flattened into one row, the features of one customer. A row from a table with columns age, income and months_active becomes [34, 52000, 18]. An embedding from a model is 768 or 1536 numbers, and none of them individually means anything a person can read. Relationships between vectors do: which ones point the same way, which are far apart.

A matrix is a grid of numbers, rows and columns, written as a list of rows. Its shape is two numbers, rows first: a 2 by 3 matrix has two rows of three numbers each. Two very different things are stored this way and it is worth keeping them apart in your head. A batch of data is a matrix with one row per example and one column per feature, so 32 by 768 is thirty-two items, each described by 768 numbers. A weight matrix is a transformation, and its shape says what it turns into what: 768 by 512 takes something of width 768 and produces something of width 512.

scalar        7                        a single number, shape ()
vector        [2, 1]                   shape (2,)      one thing
matrix        [[1, 2, 3], [4, 5, 6]]   shape (2, 3)    2 rows, 3 columns
tensor        anything with more axes  shape (32, 128, 768)   batch, tokens, width

Read a shape right to left: the last number is the width of one item,
and everything before it counts how many items there are.

The dot product#

Multiply corresponding elements and add them up. That single operation is how similarity is measured everywhere in this field. A large dot product means the vectors point in a similar direction.

Do one by hand. Take a = [2, 1] and b = [4, 2]. The dot product is 2 times 4 plus 1 times 2, which is 10. Now take c = [-1, 2]: 2 times -1 plus 1 times 2 is 0. A dot product of zero means the two vectors are at right angles, and in the geometric picture that is exactly what it means: a and c share no direction at all. If the dot product is negative the vectors point away from one another.

The problem with the raw dot product is that it grows with length. The vector [20, 10] points in exactly the same direction as a = [2, 1], and its dot product with a is 50 rather than 10, purely because it is ten times longer. Cosine similarity divides that out: divide the dot product by both lengths and you get 1.0 for both b and [20, 10], because both point the same way as a. Cosine similarity is not a different measure from the dot product. It is the dot product of the two vectors after each has been shrunk to length 1.

import numpy as np

a = np.array([2.0, 1.0])
b = np.array([4.0, 2.0])        # same direction, twice as long
c = np.array([-1.0, 2.0])       # at right angles
long = np.array([20.0, 10.0])   # same direction, ten times as long

print(a @ b, a @ c, a @ long)                     # 10.0 0.0 50.0
print(np.linalg.norm(a))                          # 2.2360679... = sqrt(5)

def cosine(u, v):
    return (u @ v) / (np.linalg.norm(u) * np.linalg.norm(v))

print(cosine(a, b), cosine(a, c), cosine(a, long))
# 0.9999999999999998 0.0 0.9999999999999999

import numpy as np

a = np.array([1, 2, 3])
b = np.array([2, 3, 4])

np.dot(a, b)          # 1*2 + 2*3 + 3*4 = 20

# Cosine similarity is the dot product with length divided out,
# which is why it compares direction rather than magnitude.
np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b))

Matrix multiplication, done by hand#

A matrix multiply is many dot products at once. Each output entry is one row of the first matrix dotted with one column of the second. That is why the inner dimensions must match: you cannot dot vectors of different lengths.

Nothing replaces doing one. Take A as 2 by 3 and B as 3 by 2. The result is 2 by 2, and each of its four entries is a separate dot product between a row of A, which has three numbers, and a column of B, which also has three.

A = [[1, 2, 3],        B = [[ 7,  8],
     [4, 5, 6]]             [ 9, 10],
                            [11, 12]]
     2 x 3                   3 x 2        ->  result is 2 x 2

entry [0,0] = row 0 of A . column 0 of B = 1*7  + 2*9  + 3*11 = 7 + 18 + 33 = 58
entry [0,1] = row 0 of A . column 1 of B = 1*8  + 2*10 + 3*12 = 8 + 20 + 36 = 64
entry [1,0] = row 1 of A . column 0 of B = 4*7  + 5*9  + 6*11 = 28 + 45 + 66 = 139
entry [1,1] = row 1 of A . column 1 of B = 4*8  + 5*10 + 6*12 = 32 + 50 + 72 = 154

A @ B = [[ 58,  64],
         [139, 154]]

import numpy as np

A = np.array([[1, 2, 3], [4, 5, 6]])          # 2 x 3
B = np.array([[7, 8], [9, 10], [11, 12]])     # 3 x 2

print(A @ B)
# [[ 58  64]
#  [139 154]]

print(B @ A)          # legal too, and a different size and meaning
# [[ 39  54  69]
#  [ 49  68  87]
#  [ 59  82 105]]

Two consequences follow from that hand calculation and they cover most of what you need. The inner dimensions must agree, because a dot product needs two lists of the same length. And the cost is one multiply-add per entry per inner element: an m by n times n by p multiply does m times n times p multiply-adds. A 4096 by 4096 weight matrix against a batch of 32 rows is about 537 million of them, which is why this operation, and not anything else in a network, is what GPUs are bought for.

(m x n) @ (n x p)  ->  (m x p)

The two n values must be equal. The result takes the outer pair.

(32 x 768) @ (768 x 512)  ->  (32 x 512)      valid
(32 x 768) @ (512 x 768)  ->  error           inner dims differ

Shape errors are the most common bug in machine learning code, and they are common for a reason that has nothing to do with the maths being hard. Every array in the program is a grid of anonymous numbers, so nothing except the shape says whether you are holding a batch of examples or a weight matrix, and nothing stops you passing one where the other belongs. The error, when you get one, names the two shapes and nothing else. Read it as a question: which of these two did I mean to transpose?

A = np.array([[1, 2, 3], [4, 5, 6]])     # 2 x 3
try:
    A @ A
except ValueError as e:
    print(e)
# matmul: Input operand 1 has a mismatch in its core dimension 0,
# with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 3)

print((A @ A.T).shape)      # (2, 2)   row against row: 2 examples compared
print((A.T @ A).shape)      # (3, 3)   column against column: 3 features compared

Why a neural network is matrix multiplies#

A layer takes a batch of inputs, multiplies by a weight matrix, adds a bias and applies a non-linearity. The whole network is that chain repeated. Once you see it this way, the architecture diagrams become readable.

# batch of 32 items, each 768 numbers
x = np.random.randn(32, 768)
W = np.random.randn(768, 512)     # weights: in_features x out_features
b = np.random.randn(512)

h = np.maximum(0, x @ W + b)      # (32, 512), ReLU applied

Three things are worth noticing in those four lines. The batch dimension, 32, is untouched by the multiply and comes out the other side unchanged: the same weights are applied to every item independently, which is why a batch is just a way of keeping the hardware busy. The weight matrix is written in_features by out_features, so reading W as 768 by 512 tells you the layer narrows its input. And the bias has shape (512,) while the result has shape (32, 512), so adding them is not a matrix operation at all. It is broadcasting, which is the next thing to understand.

Broadcasting#

NumPy stretches smaller arrays to match larger ones when the shapes are compatible. It is what lets you add a bias of shape (512,) to a result of shape (32, 512). It is also a common source of silent bugs, because an unintended broadcast produces a valid array of the wrong shape rather than an error.

a = np.ones((32, 512))
b = np.ones(512)
(a + b).shape         # (32, 512)  intended

c = np.ones((32, 1))
(a * c).shape         # (32, 512)  probably intended

d = np.ones((1, 32))
(a.T * d).shape       # (512, 32)  almost certainly a bug

The rule NumPy applies is mechanical: line the shapes up from the right, and for each pair of axes either the two must be equal or one of them must be 1, in which case that one is repeated. So (32, 512) and (512,) line up as 512 against 512, fine. (32, 512) and (32,) line up as 512 against 32, which is neither equal nor 1, so that raises an error — a bias vector applied along the wrong axis crashes, which is the good case.

The bad case is the one below, and it is worth more attention than anything else on this page because it produces no error at all. A vector of predictions with shape (3,) minus a column of targets with shape (3, 1) lines up as 3 against 1, so the 1 is stretched: you get a 3 by 3 matrix of every prediction minus every target. The mean of its squares is a real number, computed without complaint, and it is wrong.

pred = np.array([1.0, 2.0, 3.0])          # shape (3,)
true = np.array([[1.0], [2.0], [3.0]])    # shape (3, 1) -- a column, from a dataframe

err = pred - true
print(err.shape)                # (3, 3)   not (3,)
print((err ** 2).mean())        # 1.3333333333333333

print(((pred - true.ravel()) ** 2).mean())   # 0.0 -- the predictions are perfect

The habit that prevents this costs nothing: after any operation that combines two arrays, print the shape of the result and check it against the shape you expected. A loss should be a scalar, a prediction vector should have one number per example. If a loss has a shape at all, you have this bug.

Where each of these shows up in real code#

operation             where you meet it
dot product           an embedding search scoring a query against a document
cosine similarity     the same search when document lengths differ
matrix multiply       every layer of every network: x @ W
the transpose         Q @ K.T in attention, scoring every token against every token
broadcasting          adding a bias, subtracting a mean, applying a mask
shapes                the error message you will read more often than any other

How far to go#

For ML engineering: shapes, transposes, broadcasting, dot products and matrix multiplication. That is enough to read papers and debug training. Eigenvalues, decompositions and change of basis matter for research and rarely for engineering.

What is worth memorising is short: the inner dimensions must match and the outer two survive; the dot product is elementwise multiply then sum; cosine is that with both lengths divided out; broadcasting stretches an axis of length 1. Everything else — determinants, matrix inverses computed by hand, Cramer's rule, row reduction — is what a first-year linear algebra course spends its time on and what you can skip entirely. If you ever need an inverse, and you rarely should, a solver does it better than the formula.

This is the part of linear algebra that turns up in papers rather than in code you write yourself. The section covers four ideas in the terms you need in order to read them: measuring the size of a vector, projecting one vector onto another, the directions a matrix leaves pointing the same way, and writing any matrix as a rotation, a stretch and another rotation. The last one earns its place on its own, because keeping only the largest few stretch factors is the reason LoRA works at all.

Norms#

A norm measures the size of a vector. The L2 norm is ordinary length. The L1 norm is the sum of absolute values. The distinction matters because regularisation penalises one or the other, and they behave differently: L2 shrinks every weight a little, L1 drives many weights to exactly zero.

import numpy as np

v = np.array([3, 4])

np.linalg.norm(v)          # 5.0   L2, the length
np.linalg.norm(v, ord=1)   # 7.0   L1, sum of absolute values
np.linalg.norm(v, ord=np.inf)   # 4.0   the largest single entry
v / np.linalg.norm(v)      # a unit vector in the same direction

Take the difference between L1 and L2 seriously, because it explains a result people otherwise memorise. Suppose two weights must together supply a fixed amount of signal. Under L2 the penalty is the sum of squares, and 0.5 squared plus 0.5 squared is 0.5, while 1.0 squared plus 0.0 squared is 1.0 — so L2 prefers the two weights shared evenly. Under L1 the penalty is the sum of absolute values, and both arrangements cost exactly 1.0, so L1 is indifferent between spreading and concentrating. That indifference at equal cost, combined with the corner in the shape of the L1 penalty, is why L1 pushes weights to exactly zero and L2 only shrinks them.

Gradient norms are where you meet this daily. Gradient clipping compares the L2 norm of all gradients stacked into one long vector against a threshold and scales them down if it is exceeded, so the direction is kept and only the step size is capped.

Projection#

Projecting a vector onto another finds how much of it points in that direction. This is what a dot product with a unit vector gives you, and it underlies PCA, least squares and the geometry of embeddings.

The arithmetic is short enough to do in your head. Project a = [3, 4] onto the first axis, direction d = [1, 0], which already has length 1. The dot product is 3, so the projection is 3 times [1, 0], which is [3, 0]. What is left over, [0, 4], is the part of a that the direction d cannot explain, and it is at right angles to d. Every least squares fit is that split: the part of the target your features can reach, and a residual perpendicular to all of them.

a = np.array([3.0, 4.0])
d = np.array([2.0, 2.0])
u = d / np.linalg.norm(d)          # unit vector along d

projection = (a @ u) * u
print(a @ u)                        # 4.9497474683058325  how far along d
print(projection)                   # [3.5 3.5]
print(a - projection)               # [-0.5  0.5]  the residual
print((a - projection) @ d)         # 3.55e-15  perpendicular, as it must be

Eigenvectors#

Most vectors change direction when a matrix is applied to them. An eigenvector does not: the matrix only stretches it, by a factor called the eigenvalue. For a covariance matrix, the eigenvectors are the directions in which the data varies, and the eigenvalues say by how much. That is exactly what PCA computes.

Here is the smallest case where you can see it without a computer. Take the matrix [[2, 1], [1, 2]] and apply it to [1, 1]: the first row gives 2 times 1 plus 1 times 1, which is 3, and the second row gives the same. So the result is [3, 3], which is 3 times the vector you started with. The direction did not move; only the length did. Now apply it to [1, -1]: the first row gives 2 minus 1, which is 1, the second gives 1 minus 2, which is -1, so the result is [1, -1] itself, unchanged. Two eigenvectors, eigenvalues 3 and 1. Every other vector is a mixture of these two, and the matrix stretches the mixture by 3 in one direction and leaves it alone in the other.

B = np.array([[2.0, 1.0], [1.0, 2.0]])

print(B @ np.array([1.0, 1.0]))     # [3. 3.]  = 3 x the input
print(B @ np.array([1.0, -1.0]))    # [ 1. -1.]  = 1 x the input

values, vectors = np.linalg.eigh(B)
print(values)                        # [1. 3.]  smallest first
print(vectors)                       # columns are the directions, normalised to length 1
# [[-0.70710678  0.70710678]
#  [ 0.70710678  0.70710678]]

Singular value decomposition#

Any matrix can be written as a rotation, a scaling along axes, and another rotation. The scaling factors are the singular values. Keeping only the largest few gives the best possible low-rank approximation of the matrix, which is one of the most useful facts in applied linear algebra.

Do the smallest one. The same matrix [[2, 1], [1, 2]] has singular values 3 and 1 — for a symmetric matrix with positive eigenvalues they coincide with the eigenvalues. Reconstruct using only the first singular value and you get a matrix of four 1.5s, and the error, measured as the square root of the summed squared differences, is exactly 1.0, which is the singular value you threw away. That is the rule in general: the error of the best rank-k approximation is governed by the singular values you dropped. The singular values are a list of how much structure each direction carries, in order.

U, S, Vt = np.linalg.svd(B)
print(S)                                      # [3. 1.]

rank1 = S[0] * np.outer(U[:, 0], Vt[0])
print(rank1)                                  # [[1.5 1.5]
                                              #  [1.5 1.5]]
print(np.linalg.norm(B - rank1))              # 1.0  exactly the discarded singular value

The reason this matters for large matrices is that real weight matrices are not random. Build a 256 by 256 matrix as the product of a 256 by 8 and an 8 by 256 — by construction it has rank 8 — then add small noise, and look at how much of it the first few singular values account for.

rng = np.random.default_rng(0)

noise = rng.standard_normal((256, 256))                      # no structure at all
low_rank = rng.standard_normal((256, 8)) @ rng.standard_normal((8, 256))
W = low_rank + 0.01 * rng.standard_normal((256, 256))        # structure plus a little noise

for name, M in (("random", noise), ("structured", W)):
    S = np.linalg.svd(M, compute_uv=False)
    total = (S ** 2).sum()
    shares = [round(float((S[:k] ** 2).sum() / total), 3) for k in (4, 8, 32, 128)]
    print(name, shares)

# random     [0.058, 0.113, 0.378, 0.895]
# structured [0.614, 1.0, 1.0, 1.0]

U, S, Vt = np.linalg.svd(W, full_matrices=False)

k = 16
W_approx = U[:, :k] @ np.diag(S[:k]) @ Vt[:k, :]   # rank 16

# If W is 1024 x 1024, this stores 2 x 1024 x 16 + 16 numbers
# instead of a million, and may lose surprisingly little.

Count it for the 256 by 256 case above. The full matrix is 65,536 numbers. Storing it as a rank-8 product is 256 times 8 twice, which is 4,096 numbers, or 6.25 per cent — and since the true rank was 8, nothing is lost beyond the noise you added. Whether a real weight matrix behaves like the random one or the structured one is an empirical question, and the answer that made parameter-efficient fine-tuning possible is that the update a fine-tune applies behaves much more like the structured one than anybody expected.

What to skip#

You do not need to compute an SVD or an eigendecomposition by hand, ever, and no engineering interview will ask. Characteristic polynomials, Gram-Schmidt by hand, Jordan normal form: skip all of it. What you need is the reading: singular values are a sorted list of how much each direction carries, truncating the list is the best approximation you can make at that size, and eigenvectors of a covariance matrix are the directions the data spreads along. Those three sentences are the whole return on this section.

Linear algebra is arithmetic on whole lists of numbers at once. Since a model represents everything — a word, an image, a user — as a list of numbers, it is the arithmetic of the entire field. This section covers the handful of operations that actually recur: combining two lists into a single number that says how alike they are, using a grid of numbers to move points around, and finding the line that best fits a cloud of them. Each is done by hand on small numbers first and then in NumPy, so the meaning and the code arrive together.

A vector, and the dot product#

import numpy as np

a = np.array([3.0, 4.0])
b = np.array([4.0, 3.0])

print(a @ b)                  # 3*4 + 4*3 = 24.0   the dot product
print(np.linalg.norm(a))      # sqrt(3^2 + 4^2) = 5.0   the length
cos = (a @ b) / (np.linalg.norm(a) * np.linalg.norm(b))
print(cos)                    # 0.96: pointing in nearly the same direction

The dot product multiplies matching entries and adds them up. Divided by both lengths, it gives the cosine of the angle between the vectors: 1 for the same direction, 0 for unrelated, -1 for opposite. This is exactly how embedding similarity works.

Check the arithmetic against the output before going further. a is [3, 4] and b is [4, 3], so the dot product is 3 times 4 plus 4 times 3, which is 24. The length of a is the square root of 9 plus 16, which is the square root of 25, which is 5, and b has the same length by symmetry. So the cosine is 24 divided by 25, which is 0.96 — close to 1, because the two vectors are nearly the same direction but not quite.

The reason to write this as a @ b rather than a Python loop is not elegance. NumPy hands the work to compiled code that operates on the whole array at once, and the difference is not subtle.

import time

rng = np.random.default_rng(0)
u = rng.standard_normal(1_000_000)
v = rng.standard_normal(1_000_000)

start = time.perf_counter()
total = 0.0
for i in range(len(u)):
    total += u[i] * v[i]
loop_seconds = time.perf_counter() - start

start = time.perf_counter()
fast = u @ v
numpy_seconds = time.perf_counter() - start

print(round(loop_seconds, 3), round(numpy_seconds, 6))   # 0.219 0.000359
print(round(loop_seconds / numpy_seconds))               # 611, and 464 on a second run
print(total, fast)   # -204.34915522872583 -204.34915522874348

A matrix moves points#

theta = np.pi / 2
rotate = np.array([[np.cos(theta), -np.sin(theta)],
                   [np.sin(theta),  np.cos(theta)]])
stretch = np.array([[2.0, 0.0],
                    [0.0, 0.5]])

p = np.array([1.0, 0.0])
print(rotate @ p)             # [0, 1]: turned a quarter circle
print(stretch @ p)            # [2, 0]: doubled along x
print(stretch @ rotate @ p)   # rotate first, then stretch: [0, 0.5]

Multiplying a matrix by a vector transforms it: rotates, stretches, projects. Multiplying two matrices makes one transformation that does both, and order matters. A neural network layer is a learned transformation of this kind followed by a non-linearity.

Two details in that output are worth stopping on. The first entry of the rotation result prints as 6.123234e-17 rather than 0, because the cosine of pi over 2 is not exactly zero in floating point; treat anything of that size as zero and never test a float with an equality. And read the last two lines as composition applied right to left: stretch @ rotate @ p rotates first and then stretches, giving [0, 0.5], whereas rotate @ stretch @ p stretches first and gives [0, 2]. Same two operations, results four times apart. This is the same fact as A @ B not equalling B @ A, seen geometrically.

print(rotate @ stretch @ p)      # [1.2246468e-16 2.0000000e+00]   stretch first
print(stretch @ rotate @ p)      # [1.2246468e-16 5.0000000e-01]   rotate first

One more thing about that stretch matrix is worth noticing, because it is the bridge to everything in the next section. It leaves the direction [1, 0] pointing exactly where it was and only doubles its length, and it does the same to [0, 1], halving it. Those two directions are the eigenvectors of that matrix, and 2 and 0.5 are its eigenvalues. Eigenvectors are not an abstraction added later; they are the directions a transformation does not turn.

Matrix multiplication, one entry at a time#

A = np.array([[1, 2], [3, 4], [5, 6]])    # 3 x 2
B = np.array([[7, 8, 9], [10, 11, 12]])   # 2 x 3

C = np.zeros((3, 3))
for i in range(3):
    for j in range(3):
        C[i, j] = A[i, :] @ B[:, j]       # row i of A dotted with column j of B

print(np.allclose(C, A @ B))              # True

Solving equations#

# 2x + y = 5  and  x - y = 1
A = np.array([[2.0, 1.0], [1.0, -1.0]])
b = np.array([5.0, 1.0])
print(np.linalg.solve(A, b))     # [2, 1]

Substitute to check: 2 times 2 plus 1 is 5, and 2 minus 1 is 1. Both hold, so [2, 1] is the answer. Writing the system as a matrix and a right-hand side is what turns two lines of school algebra into one call that scales to ten thousand unknowns.

It fails in exactly one interesting way, and the failure is worth meeting deliberately. If one equation is a multiple of another, they say the same thing and there is no unique answer.

singular = np.array([[1.0, 2.0], [2.0, 4.0]])    # row 2 is row 1 doubled
try:
    np.linalg.solve(singular, np.array([1.0, 2.0]))
except np.linalg.LinAlgError as e:
    print(e)        # Singular matrix

Never invert a matrix to solve a system. np.linalg.inv(A) @ b gives the same answer in theory and is slower and less accurate in practice, because it computes more than you asked for and loses precision doing it. Reach for solve, or for lstsq when the system is overdetermined.

Least squares: when there is no exact answer#

With more equations than unknowns, such as fitting a line through a hundred noisy points, there is usually no exact solution. Least squares finds the one that minimises the total squared error. This is linear regression.

rng = np.random.default_rng(0)
x = rng.uniform(0, 10, 100)
y = 3 * x + 2 + rng.normal(0, 1, 100)       # true slope 3, intercept 2, plus noise

X = np.c_[x, np.ones_like(x)]               # a column of ones gives the intercept
coef, *_ = np.linalg.lstsq(X, y, rcond=None)
print(coef)                                  # [2.9869275  1.99254002]

The column of ones is worth understanding rather than copying. Without it the model is y equals slope times x, a line forced through the origin. The ones column gives every row the same constant input of 1, so its weight is free to be whatever the intercept needs to be. Every framework does this internally under the name bias, and it is the same trick.

There is a closed-form solution, the normal equations, and it is worth computing once to see that lstsq is not doing anything mysterious. It is also worth not using afterwards, for the reason above: it inverts a matrix.

normal = np.linalg.inv(X.T @ X) @ X.T @ y
print(np.abs(coef - normal).max())      # about 1e-14: the same answer

residual = y - X @ coef
print(X.T @ residual)      # about [4e-13, 4e-14]: zero, to floating point

What to skip#

Determinants beyond knowing that zero means not invertible, computing an inverse by hand, row-reducing by hand, and every identity involving cofactors: none of it will appear in your work or in an engineering interview. What is worth having at your fingertips is this section's five calls — the dot product, the norm, the matrix product, solve and lstsq — and the ability to say what each one means geometrically when asked.

Data with hundreds of columns cannot be plotted and is expensive to store, and most of that width is usually redundant: the points really lie close to a much smaller number of directions. Principal component analysis finds those directions so you can keep a few of them instead of all the original columns. This section builds it out of the pieces from the previous one, then uses it on embeddings — to draw them on a flat plot, and to cut their size without losing much.

The goal#

You have data with many dimensions, such as 768-dimensional embeddings. You want the few directions along which the data varies most, to plot it, compress it, or remove noise. Those directions are the principal components.

A principal component is a direction, written as a vector with one number per original feature. Those numbers are the recipe for mixing the original features into a new one. If your features are height and weight and the first component is [0.7, 0.7], that component is roughly overall size, and each person's score on it is 0.7 times their height plus 0.7 times their weight. The components are ordered so that the first captures as much of the spread in the data as any single direction can, the second captures as much of what is left as any direction perpendicular to the first can, and so on down.

In five lines#

import numpy as np

def pca(X, k):
    Xc = X - X.mean(axis=0)                      # centre each feature
    U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
    components = Vt[:k]                          # top k directions, one per row
    explained = (S ** 2) / (S ** 2).sum()        # share of variance per direction
    return Xc @ components.T, components, explained[:k]

The SVD splits the centred data into directions ordered by how much variance they capture. Projecting onto the first k gives the best k-dimensional summary of the data in the least squares sense.

Five points, by hand#

Run it on something small enough to check. Five points that lie almost on the diagonal line y equals x: the spread is nearly all along one direction, so the first component should be about [0.71, 0.71] and it should account for almost all the variance.

X = np.array([[1.0, 1.0],
              [2.0, 2.1],
              [3.0, 2.9],
              [4.0, 4.2],
              [5.0, 4.8]])

print(X.mean(axis=0))              # [3. 3.]  the centre the data is measured from

Z, comps, explained = pca(X, 2)
print(comps.round(4))
# [[ 0.7162  0.6979]     first component: along the diagonal
#  [ 0.6979 -0.7162]]    second: at right angles to it
print(explained.round(4))          # [0.9976 0.0024]
print(Z[:, 0].round(4))            # [-2.8282 -1.3443 -0.0698  1.5537  2.6886]

Read those last five numbers as the point of the whole method. The middle point scores -0.07, meaning it sits essentially at the centre; the first point scores -2.83, far down the diagonal. Two columns became one, and 0.24 per cent of the variation was thrown away with them. Reconstruct from that single component and the largest error on any coordinate is 0.116, which is the size of the wobble off the line that you deliberately discarded.

Now break it, in the way people break it in practice. Multiply the first column by 1000, as though it were measured in grams rather than kilograms. Nothing about the data's structure changed, and the answer changes completely.

Y = np.c_[X[:, 0] * 1000, X[:, 1]]     # same data, first column in different units

print(pca(Y, 2)[1].round(3))
# [[ 1.     0.001]      the first component is now simply 'the first column'
#  [ 0.001 -1.   ]]
print(pca(Y, 2)[2].round(3))           # [1. 0.]

How many components#

The honest answer comes from the cumulative explained variance, not from a rule. Build something embedding-shaped to see the pattern: 2,000 items described by 300 numbers, but generated from only 20 underlying directions with noise on top.

rng = np.random.default_rng(0)
latent = rng.standard_normal((2000, 20))            # 20 real underlying directions
basis = rng.standard_normal((20, 300))
embeddings = latent @ basis + 0.5 * rng.standard_normal((2000, 300))

Z, comps, explained = pca(embeddings, 300)
cumulative = np.cumsum(explained)
print(cumulative[[1, 9, 19, 49]].round(3))    # [0.147 0.609 0.989 0.991]
print(int(np.searchsorted(cumulative, 0.9)) + 1)   # 18 dimensions reach 90%

Two things follow. Two components, the number you need for a plot, cover only 14.7 per cent here, so a two-dimensional picture of a corpus is a heavy summary and should be read as one. And where the curve flattens tells you roughly how many directions the data actually has, which for real embedding models is far fewer than the 768 or 1536 they are stored in.

Check against a library#

from sklearn.decomposition import PCA

ref = PCA(n_components=2).fit(X)
print(np.abs(ref.components_).round(4))
# [[0.7162 0.6979]
#  [0.6979 0.7162]]
print(np.allclose(np.abs(ref.components_), np.abs(pca(X, 2)[1])))   # True; signs may flip

The absolute values are needed because a component pointing one way and the same component pointing the other way describe the same line, so any library may return either. If you ever compare two PCA runs and half the signs disagree, that is the reason and not a bug. Run the same comparison on the 2,000 by 300 embeddings above and it prints False, with a largest disagreement of about 0.004: scikit-learn switches to a different solver on data that shape, and where two components have nearly equal variance the directions it returns are rotated slightly against yours. Components with close singular values are not individually stable, which is a reason not to interpret any single one of them too hard.

You have the operations; this is where they are actually used. A transformer layer is almost entirely matrix multiplication, and following the shapes through one layer — what size goes in, what comes out at each step — is the quickest way to stop a language model being a black box. It also makes the cost visible, which is what lets you see in a line of arithmetic why LoRA fine-tuning is so much cheaper than updating every weight.

Tokens become vectors#

An embedding table is a matrix with one row per vocabulary token. Looking up a token means picking its row. A sentence of 10 tokens becomes a 10 by d matrix.

It is worth seeing that a lookup and a matrix multiply are the same operation here, because papers write it both ways and beginners assume they are two different mechanisms. Multiplying a one-hot row — all zeros except a single 1 in position i — by the embedding table selects row i and nothing else, since every other row is multiplied by zero. Implementations index instead of multiplying because indexing is enormously cheaper, but the maths written in the paper is the multiply.

import numpy as np

table = np.array([[0.1, 0.2, 0.3],      # token 0
                  [1.0, 0.0, 0.0],      # token 1
                  [0.5, 0.5, 0.5]])     # token 2

one_hot = np.array([0.0, 1.0, 0.0])     # 'this is token 1'
print(one_hot @ table)                   # [1. 0. 0.]
print(table[1])                          # [1. 0. 0.]  the same thing, done cheaply

tokens = [2, 0, 1, 1]                    # a four-token sequence
print(table[tokens].shape)               # (4, 3)  four rows of width d

Attention is three matrix products#

X        [tokens, d]            the input
Q = X Wq [tokens, d_head]       what each token is looking for
K = X Wk [tokens, d_head]       what each token offers
V = X Wv [tokens, d_head]       what each token passes on

scores  = Q K^T / sqrt(d_head)  [tokens, tokens]   every token against every token
weights = softmax(scores)       each row sums to 1
out     = weights V             [tokens, d_head]   a weighted mix of values

The tokens by tokens score matrix is why attention cost grows with the square of the sequence length. Put numbers on it. At 8,192 tokens that matrix has 8,192 squared entries, which is 67 million numbers, or 134 megabytes in 16-bit precision — for one head, in one layer, for one sequence in the batch. Double the context and it quadruples. Every long-context technique you will read about exists to avoid holding that matrix.

Counting the parameters of a layer#

Model sizes stop being folklore once you can count them, and the count is four matrix shapes. Take d equal to 4,096. Attention holds four d by d matrices, for queries, keys, values and the output projection: four times 4,096 squared, which is 67.1 million. The feed-forward block holds matrices that widen to 4d and come back, and in the gated form used by most current models that is three matrices of d by 4d: three times 4,096 times 16,384, which is 201.3 million. So one layer is 268.4 million parameters, and sixteen d squared is a useful thing to remember as the shape of that number.

d, layers, vocab = 4096, 32, 128256

attention = 4 * d * d              # Q, K, V and the output projection
mlp = 3 * d * (4 * d)              # gated feed-forward: two up, one down
per_layer = attention + mlp

print(per_layer / 1e6)             # 268.435456 million per layer
print(layers * per_layer / 1e9)    # 8.589934592 billion in the layers
print(vocab * d / 1e6)             # 525.336576 million in the embedding table

That count is also the memory bill. At two bytes per parameter in 16-bit precision, 8.6 billion parameters is 17.2 gigabytes just to hold the weights, before activations, before the key-value cache, before the optimiser state that training adds on top. When someone says a model will not fit on a 24 gigabyte card, this multiplication is the whole of what they did.

Low rank: why LoRA is cheap#

A full weight update for a 4096 x 4096 matrix:   16,777,216 numbers

LoRA writes the update as B A, with B [4096 x 8] and A [8 x 4096]:
  4096 x 8 + 8 x 4096 = 65,536 numbers   (about 0.4% of the full update)

A product of a tall thin matrix and a short wide one has rank at most 8 here: it can only change the weights along 8 directions. The bet behind LoRA is that fine-tuning mostly needs a few directions, and in practice that often holds.

Scale that to the whole model with the count from above. Applying rank-8 adapters to all four attention matrices in all 32 layers is 32 times 4 times 65,536, which is 8.4 million trainable parameters against 8.6 billion in the model: under a tenth of one per cent. That is the difference between a fine-tune that needs a cluster and one that runs on a single card, and it comes entirely from the shape of two matrices.

rank = 8
lora_per_matrix = 2 * d * rank                    # 65,536
lora_total = layers * 4 * lora_per_matrix         # four attention matrices per layer

print(lora_per_matrix)                             # 65536
print(lora_total / 1e6)                            # 8.388608 million trainable
print(100 * lora_total / (layers * per_layer))     # 0.09765625 per cent

One caution on the arithmetic: fewer trainable parameters is not proportionally less compute. The forward pass still runs every one of the 8.6 billion weights, and the backward pass still has to propagate gradients through the frozen layers to reach the adapters. What LoRA removes is optimiser state and gradient storage for the frozen weights, which with Adam is the largest part of a training run's memory, and it makes each fine-tuned variant a few megabytes to store rather than a few gigabytes.