PCA, and the next step

Back in Lesson 1 we said a search engine turns a sentence into a vector of 768 numbers. Data like that cannot be drawn, and the more columns it has, the heavier the computation gets. But in real data the columns are often tangled together, the way taller people also tend to weigh more. So instead of carrying many numbers around separately, you can pick just a few directions in which the data spreads the widest, keep only the coordinates along them, and lose almost no information.

The method that turns this idea into a recipe is principal component analysis (PCA). It is the first tool people reach for in data visualization (a 768-D → 2-D map), machine learning preprocessing, noise removal, compressing face photos, and more.

PCA brings almost nothing new. Thread together the dot product and shadows from Lesson 3, change of basis from Lesson 10, the eigenvectors of symmetric matrices from Lesson 11, projection from Lesson 13 and the SVD from Lesson 14, and you have PCA. This lesson ties the whole course together that way, and at the end previews the next course, numerical analysis.

15-1. Data is a matrix (5 min)

A data table is a matrix as it stands. There is one convention: one row = one sample (a person, a photo, a sentence), one column = one feature (height, weight, brightness). Let's write down the height (cm) and weight (kg) of five people.

Each row is a 2-D vector, that is, a single point on the plane. The whole data set is a point cloud.

>>> import numpy as np
>>> X = np.array([[170, 65], [160, 52], [180, 80], [175, 70], [165, 60]])
>>> X.shape
(5, 2)

Read (5, 2) as "5 samples, 2 features". From here on, a data matrix always has this shape (number of samples, number of features).

axis: which way to fold

When you take a mean, you have to decide which direction to average along. That is what axis is for.

>>> X.mean(axis=0)
array([170. ,  65.4])
>>> X.mean(axis=1)
array([117.5, 106. , 130. , 122.5, 112.5])
>>> X.mean()
np.float64(117.7)

axis=0 folds away axis 0 (the direction the rows are stacked in, i.e. top to bottom). The five rows fold into one, leaving one mean per column: "mean height 170, mean weight 65.4". axis=1 folds axis 1 (left to right), giving each person "the mean of their height and weight", which is a meaningless number here. Leave out axis and you get the mean of all 10 numbers.

A way to remember it: axis=k removes position k of the shape. On (5, 2), axis=0 leaves (2,) and axis=1 leaves (5,). In PCA it is almost always axis=0 (per feature).

15-2. Subtracting the mean, and covariance (8 min)

Centering. First subtract the mean from every point. This slides the whole point cloud so its center sits at the origin (Lesson 2: an arrow moved elsewhere is the same vector). Only then can you measure "which way it spreads" with vectors reaching out from the origin.

Variance. How far one feature spreads from its mean is measured by the variance, "the average of the squared differences from the mean". When computing it from $n$ samples, the convention is to divide by $n-1$ instead of $n$ (just think of it as a correction for the fact that estimating the whole from a sample comes out a little too small).

Covariance. This is how much two features move together. Multiply their differences from the mean and average them. When one is above its mean and the other is too, the product is positive; when they go opposite ways, negatives pile up. It is the same picture as the sign of the dot product in Lesson 3 (same side / perpendicular / opposite). In fact, the covariance is the dot product of the two centered columns divided by $n-1$.

Gathering all of these in one place gives the covariance matrix. Calling the centered data $X_c$,

The variance of each feature sits on the diagonal, and the covariance sits off the diagonal. What matters is that it is a symmetric matrix, $C^{\mathsf T} = C$. In Lesson 11 we saw that a symmetric matrix has perpendicular eigenvectors.

Some books write the covariance matrix as $\Sigma$. It is the same letter as the singular value matrix Σ from Lesson 14, so tell them apart by context. These notes use $C$.

Let's compute it step by step in NumPy.

>>> Xc = X - X.mean(axis=0)
>>> Xc
array([[  0. ,  -0.4],
       [-10. , -13.4],
       [ 10. ,  14.6],
       [  5. ,   4.6],
       [ -5. ,  -5.4]])
>>> np.allclose(Xc.mean(axis=0), 0)
True

X has shape (5, 2) and X.mean(axis=0) has shape (2,). The subtraction works despite the different shapes thanks to broadcasting from Lesson 2: the length-2 mean is applied to each of the five rows alike. After centering, the mean should be 0, but floating-point crumbs (values like -5.7e-15) can remain, so we checked with np.allclose.

Compare the covariance matrix built from the formula with what np.cov gives.

>>> Xc.T @ Xc / (len(X) - 1)
array([[ 62.5,  82.5],
       [ 82.5, 110.8]])
>>> C = np.cov(X, rowvar=False)
>>> C
array([[ 62.5,  82.5],
       [ 82.5, 110.8]])
>>> np.var(X[:, 0], ddof=1)
np.float64(62.5)

The variance of height is 62.5, the variance of weight is 110.8, and the covariance is 82.5 (positive: taller means heavier). np.cov subtracts the mean and divides by n − 1 for you.

rowvar=False means "rows are not the variables (columns are)". np.cov defaults to rowvar=True, so if you pass the data without saying anything, it reads each row as a feature. In our table the rows are people, so you must pass rowvar=False. The ddof=1 in np.var means divide by n − 1. np.var divides by n by default, so you write this to match np.cov.

15-3. Principal components = eigenvectors of the covariance matrix (9 min)

Now for the PCA question: from which direction does the point cloud look the most spread out?

Pick a direction $\mathbf{w}$ of length 1. The length of the shadow each centered point $\mathbf{x}$ casts onto $\mathbf{w}$ is the dot product $\mathbf{x}\cdot\mathbf{w}$ from Lesson 3. The shadow lengths of all the points at once are $X_c\mathbf{w}$, and computing the variance of these shadows gives a neat formula.

The length-1 $\mathbf{w}$ that makes $\mathbf{w}^{\mathsf T}C\mathbf{w}$ largest is the eigenvector of C for its largest eigenvalue, and the value it reaches is that eigenvalue $\lambda_1$. The lab from Lesson 14 makes this easy to see. Move the unit circle with the symmetric matrix C and you get an ellipse; the direction of the long radius is the direction stretched the most, and for a symmetric matrix that direction is an eigenvector.

• PC1 (the first principal component): the eigenvector for C's largest eigenvalue $\lambda_1$. The direction in which the point cloud spreads the widest.
• PC2: the eigenvector for the next eigenvalue $\lambda_2$. Since C is symmetric, it is perpendicular to PC1.
• $\lambda_i$ = the variance measured along that direction.
• Explained variance ratio: $\dfrac{\lambda_1}{\lambda_1 + \lambda_2}$. The share of the total spread that PC1 alone holds.

The denominator $\lambda_1 + \lambda_2$ equals the sum of C's diagonal, $\operatorname{Var}(\text{x}) + \operatorname{Var}(\text{y})$. Rotate the axes and measure again, and the total spread stays the same; PCA is choosing the direction that piles as much of that spread as possible onto PC1.

The lab below holds a cloud of 60 points. The yellow arrow is PC1 and the gray arrow is PC2. Each arrow's length is twice the square root of the variance in that direction (the standard deviation), that is, 2√λ.

Things to notice in the lab:

• Dragging the black handle (the center) does not change the covariance matrix or PC1. Because the mean is subtracted first, only the shape of the spread counts, wherever the point cloud sits.
• The starting view has a tilt of 30°, yet PC1 is about (0.855, 0.518), an angle of about 31°. The 60 points are a random sample, so they do not match the direction they were built with exactly. In real data, too, what PCA finds is the widest direction "in this sample".
• Make 'Short spread' equal to 'Long spread' and the cloud turns round, with the ratio dropping to around 59%. It is not exactly 50% because the 60 points happen to spread a little more one way, and PC1 then points in that accidental direction regardless of the 'Tilt' setting. In other words, when the data has no particularly wide direction, PC1 is hard to trust. Push 'Short spread' even higher and PC1 swings around by 90°. Whatever the sliders are called, PC1 is always the side that is actually wider.
• Turn on 'Project onto PC1 (2-D → 1-D)' and each point drops straight down onto the PC1 line (Lesson 13's projection). Instead of two numbers per point, only one number remains, its position along the line: a 2-D → 1-D reduction.

Reading the NumPy panel

In the starting view, the "Check it with NumPy" panel below the lab shows this.

>>> X = np.array(points)  # the 60 points from the lab
>>> X.shape
(60, 2)
>>> np.round(X.mean(axis=0), 4)
array([0.5   , 0.5001])
>>> C = np.cov(X, rowvar=False)
>>> np.round(C, 4)
array([[1.6861, 0.8735],
       [0.8735, 0.7732]])
>>> vals, vecs = np.linalg.eigh(C)
>>> np.round(vals, 4)
array([0.244 , 2.2152])
>>> np.round(vecs[:, -1], 4)
array([-0.8553, -0.5181])
>>> np.round(vals[-1] / vals.sum(), 4)
np.float64(0.9008)
>>> z = (X - X.mean(axis=0)) @ vecs[:, -1]
>>> z.shape
(60,)

(The lab's list of points, points, only lives inside the page, so this block is for reading, not code you can type along with.)

Let's read it line by line.

• np.round(…, 4): the rounding function from Lesson 14. It shows only four decimal places so the digits do not trail on. The values themselves are unchanged; only the display is rounded. So a number with a trailing zero dropped, like 0.244, means the fourth decimal is 0.
• X.shape is (60, 2): 60 points (samples), 2 features (x, y).
• X.mean(axis=0): the center of the point cloud, which is where the black handle sits. The points are drawn with their coordinates rounded to three decimals, so it is off by a hair, as in 0.5001.
• np.linalg.eigh(C): the eigenvalue function for symmetric matrices from Lesson 11. It lists the eigenvalues smallest first. So the largest eigenvalue is the last one, vals[-1], and its eigenvector is the last column of vecs, vecs[:, -1]. That is PC1.
• vals[-1] / vals.sum(): the explained variance ratio. 0.9008 matches the lab's "90.1%". vals.sum() adds up all the numbers in the array.
• z = (X - X.mean(axis=0)) @ vecs[:, -1]: the shadow lengths of the centered points projected onto PC1. 60 points have become 60 numbers, so the shape is (60,).

The lab shows PC1 as (0.855, 0.518), and NumPy as (−0.8553, −0.5181). This is the same sign issue you saw in Lesson 14. An eigenvector times −1 is still an eigenvector, pointing along the same line. In PCA what matters is the line the direction lies on, not which end the arrowhead is at.

15-4. Reducing dimensions by projection (6 min)

Back to the height and weight data, and let's carry PCA all the way through.

>>> vals, vecs = np.linalg.eigh(C)
>>> np.round(vals, 3)
array([  0.688, 172.612])
>>> np.round(vecs, 3)
array([[-0.8,  0.6],
       [ 0.6,  0.8]])
>>> v = vecs[:, -1]
>>> round(float(vals[-1] / vals.sum()), 4)
0.996

PC1 is the direction (0.6, 0.8), "the direction in which height and weight grow together", and it explains 99.6% of the total variance on its own. PC2 is the direction (−0.8, 0.6), which measures "heavy or light for one's height". The two directions are perpendicular.

>>> bool(np.isclose(v @ vecs[:, 0], 0))
True

Now find each person's shadow length on PC1. It is the dot product from Lesson 3, as is.

>>> z = Xc @ v
>>> np.round(z, 2)
array([ -0.32, -16.72,  17.68,   6.68,  -7.32])
>>> np.round(np.var(z, ddof=1), 3)
np.float64(172.612)

Each person's two numbers (height, weight) have become one number (their PC1 coordinate). The variance of these numbers is exactly λ₁. Person 3 (17.68) is at the "tallest and heaviest" end, and person 2 (−16.72) at the opposite end.

How well can a single number bring back the original point? Start from the mean and go z along the PC1 direction. That is where the points landed when you turned on 'Project onto PC1 (2-D → 1-D)' in the lab.

>>> X_hat = X.mean(axis=0) + np.outer(z, v)
>>> np.round(X_hat, 1)
array([[169.8,  65.1],
       [160. ,  52. ],
       [180.6,  79.5],
       [174. ,  70.7],
       [165.6,  59.5]])

np.outer(z, v) is the outer product from Lesson 14: a (5, 2) matrix that stacks z[i] * v one row per person. It differs from the original table by only about 1 or 2 in each entry. What was lost is the spread along PC2, that is, λ₂ = 0.688.

Seen through the eyes of Lesson 10, PCA is a change of basis. It takes PC1 and PC2 as a new basis, re-measures the coordinates, and then throws away the coordinates with small variance. vecs is an orthogonal matrix, with perpendicular columns of length 1, so there is no need for an inverse: Xc @ vecs is the new coordinates directly.

15-5. The whole pipeline, and the same answer from SVD (9 min)

Let's build our own cloud of 200 points, similar to the lab's, and run the whole PCA in one go. Take a cloud with a horizontal spread of 2 and a vertical spread of 0.5, rotate it by 30°, and move its center to (3, 1).

>>> rng = np.random.default_rng(0)
>>> base = rng.standard_normal((200, 2)) * [2.0, 0.5]
>>> th = np.radians(30)
>>> R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
>>> X = base @ R.T + [3, 1]
>>> X.shape
(200, 2)

* [2.0, 0.5] uses broadcasting to multiply the first column by 2 and the second by 0.5. R is the rotation matrix from Lesson 6. But notice we wrote base @ R.T, not R @ base. The points are stored as rows, so to multiply each point by R, you multiply by the transpose on the right. (For a row vector p, p Rᵀ = (R p)ᵀ.)

>>> mean = X.mean(axis=0)
>>> Xc = X - mean
>>> C = np.cov(X, rowvar=False)
>>> vals, vecs = np.linalg.eigh(C)
>>> np.round(vals, 3)
array([0.254, 3.856])
>>> pc1 = vecs[:, -1]
>>> np.round(pc1, 3)
array([-0.859, -0.513])
>>> if pc1[0] < 0:
...     pc1 = -pc1
>>> round(float(np.degrees(np.arctan2(pc1[1], pc1[0]))), 1)
30.8
>>> round(float(vals[-1] / vals.sum()), 3)
0.938

The sign came out negative, so we flipped it with an if. Lining it up so that "the x component is positive" is a common habit that makes comparisons easier. np.arctan2(y, x) gives the angle of the vector (x, y) in radians, and np.degrees turns it into degrees (Lesson 3). We got 30.8°, close to the 30° we built in, and the variances (horizontal spread² = 4, vertical spread² = 0.25) roughly match too.

The same answer from SVD (the link to Lesson 14)

In Lesson 14 we saw that the eigenvectors of AᵀA are V and its eigenvalues are σ². The covariance matrix is Xcᵀ Xc divided by n − 1, so PCA falls straight out of the SVD of the centered data.

>>> U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
>>> U.shape, S.shape, Vt.shape
((200, 2), (2,), (2, 2))
>>> np.round(Vt[0], 3)
array([-0.859, -0.513])
>>> np.round(S ** 2 / (len(X) - 1), 3)
array([3.856, 0.254])
>>> np.allclose(S ** 2 / (len(X) - 1), vals[::-1])
True

• Vt[0] (the first right singular vector) is PC1 (it even came out with the same sign, but there is no guarantee it always will).
• S ** 2 / (n - 1) are the eigenvalues. SVD lists them largest first and eigh smallest first, so we reversed the order with vals[::-1] to compare. [::-1] is slicing that takes the items "one by one, backwards from the end".
• With 200 rows we passed full_matrices=False. Without it, U would be 200×200.

PCA libraries used in practice (scikit-learn, for example) do not build the covariance matrix; they compute it with SVD like this. That avoids the floating-point error introduced by multiplying out Xcᵀ Xc, which is exactly the subject of the next course, numerical analysis.

When there are many features

With 3 or more features, the code stays the same. Make a third feature out of "first + second + a little noise", and the third column carries essentially no new information.

>>> t = rng.standard_normal(200)
>>> X3 = np.column_stack([X, X[:, 0] + X[:, 1] + 0.1 * t])
>>> X3.shape
(200, 3)
>>> vals3 = np.linalg.eigvalsh(np.cov(X3, rowvar=False))
>>> np.round(vals3[::-1], 2)
array([11.11,  0.27,  0.  ])
>>> np.round(np.cumsum(vals3[::-1]) / vals3.sum(), 4)
array([0.9764, 0.9997, 1.    ])

np.linalg.eigvalsh is the version of eigh that returns only the eigenvalues. np.cumsum is the cumulative sum, adding up from the front one at a time. One principal component explains 97.6% and two explain 99.97%, so reducing this 3-D data to 2-D loses almost nothing. With real data, a common approach is to keep principal components up to the point where this cumulative ratio passes 90% or 95%.

15-6. Common mistakes (3 min)

Forgetting rowvar=False. Call plain np.cov(X) on data whose rows are samples and you get the covariance between samples. The shape gives it away immediately.

>>> X = np.array([[170, 65], [160, 52], [180, 80], [175, 70], [165, 60]])
>>> np.cov(X).shape
(5, 5)
>>> np.cov(X, rowvar=False).shape
(2, 2)

Forgetting eigh's order. np.linalg.eigh gives the smallest eigenvalue first. vecs[:, 0] is not PC1 but the narrowest direction. PC1 is vecs[:, -1].

Mixing features with different units as they are. Variance changes size with the unit. Write only the heights of the same data in meters and the variance of height shrinks to 1/10000, so PC1 points almost purely at weight.

>>> Xm = X * [0.01, 1]
>>> vals_m, vecs_m = np.linalg.eigh(np.cov(Xm, rowvar=False))
>>> np.round(vecs_m[:, -1], 3)
array([0.007, 1.   ])
>>> Xs = (X - X.mean(axis=0)) / X.std(axis=0, ddof=1)
>>> vals_s, vecs_s = np.linalg.eigh(np.cov(Xs, rowvar=False))
>>> np.round(vecs_s[:, -1], 3)
array([0.707, 0.707])

When you mix features with different units, subtract each column's mean and divide by its standard deviation (standardization) so every feature has "variance 1", then run PCA. X.std(axis=0, ddof=1) is the standard deviation of each column.

Running SVD without subtracting the mean. If you skip centering when doing PCA by SVD, the first singular vector is not "the direction of widest spread" but the direction from the origin toward the center of the point cloud.

>>> np.round(np.linalg.svd(X)[2][0], 3)
array([-0.933, -0.361])
>>> np.round(X.mean(axis=0) / np.linalg.norm(X.mean(axis=0)), 3)
array([0.933, 0.359])

Instead of PC1, (0.6, 0.8), it came out almost the same as the normalized mean (170, 65.4). np.cov subtracts the mean for you, but SVD does not.

15-7. Looking back over 15 lessons (4 min)

We started from a single bundle of numbers and made it all the way here. Here is the heart of each lesson, with one line of NumPy that puts it to work.

Looking back, two pictures run through the whole course. One is a matrix = a transformation that moves space (Lesson 5). The determinant was that transformation's area scale factor, the inverse was undoing it, eigenvectors were the axes whose direction does not change, and singular values were the radii when a circle becomes an ellipse. The other is Av = a linear combination of the columns (Lesson 4). Systems of equations, the column space and least squares were all the same question: "can b be built from A's columns, and if not, what is the closest you can get?"

With these two pictures in your head, if you open the linear algebra book you once closed, you should be able to see the picture behind the symbols.

15-8. Coming up next: numerical analysis (3 min)

There is a scene we ran into again and again in this course: 0.30000000000000004, 1.9999999999999996, np.allclose, "use solve instead of inv" (Lesson 8), "singular values can come out like 1e-17" (Lesson 14). Calculations that are clean on paper drift a little inside a computer. The next course, numerical analysis, takes on exactly this drift head-on.

>>> 0.1 + 0.2
0.30000000000000004
>>> A = np.array([[2, 1], [0.5, 1.5]])
>>> S = np.linalg.svd(A, compute_uv=False)
>>> round(float(S[0] / S[1]), 3)
2.618
>>> round(float(np.linalg.cond(A)), 3)
2.618
>>> round(float(np.linalg.cond(np.array([[1, 1], [1, 1.0001]]))))
40002

The last two lines are the condition number. It is defined with the singular values from Lesson 14 (σ₁ ÷ σ₂, how flat the ellipse is), and when it is large, even a tiny wobble in the input makes the answer wobble a lot. The second matrix, which is almost flat, has a condition number over 40,000.

Here is a preview of what the next course covers.

• Floating point: how a computer writes real numbers, and why errors appear and pile up.
• Condition numbers and stability: telling apart a problem that is sensitive in itself from a method that magnifies the error. The real reason solve beats inv(A) @ b.
• LU decomposition: what np.linalg.solve does inside (Gaussian elimination, written as matrices).
• Iterative methods: how to creep closer to the answer when the matrix is millions × millions and cannot be solved in one shot (Jacobi, Gauss–Seidel, conjugate gradient).
• The power method and the QR algorithm: refining the power method you got a taste of in Lesson 12, all the way to how np.linalg.eig actually finds eigenvalues.

If this course taught you "what it does, in pictures", the next one teaches you "how the computer does it, and how accurately".

Exercises (5 min)

1. In the lab, if you set 'Tilt' to 0°, what does the off-diagonal entry of the covariance matrix (the covariance) get close to? Which way does PC1 point then?
2. For the height and weight data X, find the variance of the PC2 coordinates (Xc @ vecs[:, 0]) and check whether it equals λ₂.
3. For X3 in 15-5, how many principal components do you need for the cumulative ratio to pass 99%?

Solutions

>>> Xc = X - X.mean(axis=0)
>>> vals, vecs = np.linalg.eigh(np.cov(X, rowvar=False))
>>> z2 = Xc @ vecs[:, 0]
>>> bool(np.isclose(np.var(z2, ddof=1), vals[0]))
True
>>> np.round(vals[0], 3)
np.float64(0.688)
>>> int(np.argmax(np.cumsum(vals3[::-1]) / vals3.sum() > 0.99)) + 1
2

Question 1: with a tilt of 0°, the cloud spreads along the x-axis and x and y do not move together, so the covariance gets close to 0, and PC1 points near (1, 0) (or (−1, 0)). The covariance matrix becomes almost diagonal, so the variances on its diagonal are the eigenvalues as they stand.

Question 2: the variance of the PC2 coordinates equals λ₂ = 0.688. The result from 15-4, that the variance of each principal component's coordinates is its eigenvalue, holds for PC2 as well.

Question 3: np.cumsum(…) > 0.99 is an array that is True wherever the cumulative ratio passes 99%, and np.argmax gives the position of the first True (counting from 0). Add 1 to turn it into a count, and you need 2.

That completes all 15 lessons of Linear Algebra You Can See, with NumPy. In the next course, numerical analysis, you will see how the solve, eig and svd you used today are computed inside the computer, and how far you can trust their answers.