SVD and PCA

SVD: turn, stretch, turn

How does a computer find the longest directions? It uses the singular value decomposition, or SVD. It works on any table of numbers, of any shape.

You've seen that a matrix moves points around. SVD says every matrix makes that move in three simple stages. It turns the space, stretches it along a few directions at right angles, then turns it again. (A "turn" here may also include a mirror flip.)

Here $X$ is your data table, with one row per customer. Read the right-hand side from right to left. $V^\top$ ("V transposed") is the first turn. Its rows are the directions. For centered data, they are the principal components.

Next comes $\Sigma$, a capital sigma. It is the stretch. It holds the singular values, one per direction, sorted from biggest to smallest. We write them $\sigma_1, \sigma_2$ and so on, with a small sigma. Last, $U$ is the final turn.

In plain words: any table of numbers can be rebuilt from a turn, a set of stretches and another turn. For centered data, the biggest stretch points along the direction where your data spreads most.

Which share of the spread?

Each direction's share of the spread comes from its singular value. This works on centered data, as in step 1. Number the directions with $i$: 1 for the first, 2 for the second. Square every singular value. A direction's share is its square divided by the total of all the squares:

In plain words: square the stretches, then ask what fraction of the total each square is. This fraction is called the explained variance ratio. All the shares add up to 1.

You'll also see the bottom written as $\sum \sigma^2$. The big $\sum$ is another sigma, used to mean "add them all up".

Here is a small example to practice the formula. The matrix below stretches the page 3 times across and 1 time up. It isn't a centered data table, but the arithmetic is the same. No turning is needed, so the singular values are 3 and 1. Square them: 9 and 1, which total 10. The shares are 9 ÷ 10 = 0.9 and 1 ÷ 10 = 0.1. The first direction gets 90% of the total.

import numpy as np

X = np.array([[3.0, 0.0],
              [0.0, 1.0]])
U, S, Vt = np.linalg.svd(X)
print("singular values:", S)
print("shares:", (S ** 2 / np.sum(S ** 2)).round(2))

It prints [3. 1.] and [0.9 0.1], the numbers you worked out by hand.

PCA on 120 customers

Now the real thing. The code below invents 120 customers with the four features from earlier. About half are bargain hunters and half are premium buyers. It standardizes the table, runs SVD and prints each component's share of the spread. It also prints the first component's weights.

import numpy as np

rng = np.random.default_rng(7)
n = 120
segment = rng.integers(0, 2, size=n)                 # 0 = bargain hunter, 1 = premium
orders = np.where(segment == 0, 14, 5) + rng.normal(0, 2, n)
value = np.where(segment == 0, 40, 120) + rng.normal(0, 12, n)
discount = np.where(segment == 0, 0.6, 0.1) + rng.normal(0, 0.08, n)
clicks = 0.8 * orders + rng.normal(0, 2, n)
X = np.column_stack([orders, value, discount, clicks])

Z = (X - X.mean(axis=0)) / X.std(axis=0)             # standardize
U, S, Vt = np.linalg.svd(Z, full_matrices=False)
explained = S ** 2 / np.sum(S ** 2)
print("explained variance ratio:", explained.round(3))
print("first component (weights on the 4 features):", Vt[0].round(2))

The first component alone explains about 90% of the spread. Its weights show what it measures. Orders, discounts and clicks pull one way. Order value pulls the other. That's the "bargain hunter or premium" story. PCA found it without ever being told which customer was which. (Your signs may all come out flipped. An arrow and its reverse lie on the same line, so the meaning is the same.)

import numpy as np

orders = np.array([2.0, 4.0, 6.0, 8.0])
clicks = np.array([3.0, 4.0, 8.0, 9.0])
X = np.column_stack([orders, clicks])

centered = X - X.mean(axis=0)
print("centered:")
print(centered)

U, S, Vt = np.linalg.svd(centered, full_matrices=False)
direction = Vt[0] * np.sign(Vt[0, 0])      # make the arrow point the positive way
print("direction:", direction.round(2))
print("share of spread:", (S ** 2 / np.sum(S ** 2)).round(2))
print("one number per customer:", (centered @ direction).round(1))

Seeing it: from 4 numbers to 2

Multiply the centered data by the top two components. You get two numbers per customer. Two numbers make a point on a flat chart, so you can plot it:

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(7)
n = 120
segment = rng.integers(0, 2, size=n)
orders = np.where(segment == 0, 14, 5) + rng.normal(0, 2, n)
value = np.where(segment == 0, 40, 120) + rng.normal(0, 12, n)
discount = np.where(segment == 0, 0.6, 0.1) + rng.normal(0, 0.08, n)
clicks = 0.8 * orders + rng.normal(0, 2, n)
Z = np.column_stack([orders, value, discount, clicks])
Z = (Z - Z.mean(axis=0)) / Z.std(axis=0)

U, S, Vt = np.linalg.svd(Z, full_matrices=False)
P = Z @ Vt[:2].T                                     # (120, 2): coordinates on PC1 and PC2
for seg, label in [(0, "bargain hunters"), (1, "premium")]:
    plt.scatter(P[segment == seg, 0], P[segment == seg, 1], s=18, label=label)
plt.xlabel("first principal component")
plt.ylabel("second principal component")
plt.title("120 customers, 4 features, shown in 2D")
plt.legend()

Two clear groups appear. PCA never saw the segment labels. The colors are added afterwards, only so you can check the result.

People plot embeddings the same way. They take 1,536 numbers down to 2, which is often enough to see topics cluster. Two directions can't hold everything, though. For pictures like this, nonlinear methods such as UMAP often separate clusters more clearly. PCA is the straight-line, exact and fast starting point.

Compression: keep the big stretches

Saving a photo at lower quality keeps the big shapes and drops the fine detail. SVD can do the same to any table. Because the singular values are sorted, you can keep only the first $k$ of them, where $k$ is a number you choose. Throw the rest away. The result is the best possible rank-$k$ approximation: the closest copy of the table you can build from $k$ directions. (The rank of a table is how many truly separate directions it holds.)

In plain words: keep the top $k$ stretches and their directions, drop the rest, and multiply back. The sign $\approx$ means "roughly equal". The small $k$ means "only the first k pieces".

The savings can be huge. A 1,000 × 1,000 table holds a million numbers. A rank-10 copy needs 1,000 × 10 = 10,000 numbers for $U_k$, 10 singular values, and 10 × 1,000 = 10,000 numbers for $V_k^\top$. That's about 20,000 numbers in all, roughly 2% of the original.

The code below builds a 200 × 50 table that secretly has only 3 directions. Then it adds a little noise. It rebuilds the table from 1, 2, 3 and 5 directions. The relative error is how far each copy is from the original, as a share of the original's size. 0 would be a perfect copy.

import numpy as np

rng = np.random.default_rng(0)
low_rank = rng.normal(size=(200, 3)) @ rng.normal(size=(3, 50))   # secretly rank 3
X = low_rank + 0.05 * rng.normal(size=(200, 50))                  # plus a little noise
U, S, Vt = np.linalg.svd(X, full_matrices=False)
for k in [1, 2, 3, 5]:
    approx = U[:, :k] * S[:k] @ Vt[:k]
    error = np.linalg.norm(X - approx) / np.linalg.norm(X)
    print(f"rank {k}: relative error {error:.3f}")

The error falls from 0.750 at rank 1 to almost nothing at rank 3. After that it barely moves. The data really only had three directions. What's left over is the noise.

For later: LoRA, the most common way to fine-tune large models cheaply, assumes the change to a huge weight matrix is low-rank. It learns that change as two thin matrices instead of the full one. Same idea as above, applied to training.

sizes = [8.0, 1.5, 0.5]
kept = sum(s * s for s in sizes[:1]) / sum(s * s for s in sizes)
print(round(kept, 2))