PCA from scratch

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.