Norms, projections and decompositions

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.