PCA from Scratch — Eigenvectors as the Axes of Your Data

A cloud of correlated data has a grain — a direction along which it varies most. Principal component analysis is the machinery that finds that grain. This article builds it from one picture: center the data, ask which direction carries the most variance, and discover that the answer is an eigenvector of the covariance matrix. Eigenvalues turn out to be the variance along each axis, the top few axes give you dimensionality reduction, and the explained-variance ratio tells you what you threw away.

On Lines and Planes of Closest Fit to Systems of Points in Space

Too many numbers, most of them redundant

A dataset with many features is a room full of columns, and most of those columns are quietly telling you the same thing twice. Height and weight move together; a customer's number of purchases and total spend move together; two nearby pixels in an image are almost always the same color. When features are correlated, the data does not really fill up all the dimensions you gave it — it huddles near a lower-dimensional shape inside that space, and most of the coordinates are wasted describing the same variation over and over.

Principal component analysis (PCA) is the tool that finds the shape the data actually lives on. It rotates your coordinate axes so they line up with the directions the data genuinely varies in, and then lets you keep only the few axes that matter. This article builds PCA in the usual order — the picture, the direction of maximum variance, the covariance matrix and its eigenvectors, projection for dimensionality reduction — and then the explained-variance ratio, code, and the two ways it bites you in practice.

The picture: a cloud has a grain

Picture a scatter of points in two dimensions that lie in a tilted, cigar-shaped cloud. The cloud is not aligned with the $x$ and $y$ axes you happened to measure in — it runs diagonally. There is clearly one direction along which the points are spread out a lot (down the long axis of the cigar) and a perpendicular direction along which they are barely spread at all (across its narrow width).

Those two directions are the principal components. The first is the direction of greatest variance; the second is the direction of greatest remaining variance perpendicular to the first; and so on in higher dimensions. PCA is nothing more than the procedure for finding these axes, in order, and the whole method falls out of taking that single sentence — "the direction of greatest variance" — completely literally.

First, center the data

Before you can talk about the direction a cloud spreads in, you have to move it to the origin. Compute the mean of each feature and subtract it, so every column now has mean zero:

$$x_i \leftarrow x_i - \bar{x}, \qquad \bar{x} = \frac{1}{n}\sum_{i=1}^{n} x_i$$

In words: find the average position $\bar{x}$ of all the points, then subtract it from each one. Here $n$ is the number of samples and $x_i$ is the $i$-th point; after the subtraction the cloud's center of mass sits exactly at the origin.

This is not a cosmetic step. Variance is spread about the mean, so if you skip the centering, the largest "direction of variation" the math finds will be dominated by how far the whole cloud sits from the origin, not by how the points scatter among themselves. Centering makes the data matrix $X$ (with $n$ rows, one per sample, and $d$ columns, one per feature) express pure deviations, which is exactly what the next step measures.

The direction of maximum variance

Now make the goal precise. We want the unit vector $u$ — a direction, so $\|u\| = 1$ — along which the projected data has the largest variance. Projecting a centered point $x_i$ onto $u$ gives the single number $u^\top x_i$, and the variance of those projections across the dataset is

$$\operatorname{Var}(u) = \frac{1}{n}\sum_{i=1}^{n} (u^\top x_i)^2 = u^\top \left(\frac{1}{n} X^\top X\right) u.$$

In words: drop each point's shadow onto the direction $u$, giving the length $u^\top x_i$, square it, and average — that is how spread out the data is along $u$. (Because the data is centered, the mean of those shadows is zero, so the mean square is the variance.) The algebra just moved $u$ outside the sum, and the thing left in the middle, $\frac{1}{n} X^\top X$, is a single $d \times d$ matrix that no longer depends on $u$. That matrix is the star of the whole show.

The covariance matrix and its eigenvectors

That middle matrix is the covariance matrix:

$$C = \frac{1}{n} X^\top X.$$

In words, it is a table of how strongly every pair of features moves together. Its diagonal entries are the variances of each feature, and its off-diagonal entries are the covariances between pairs of features — the numbers that were large precisely because features were correlated. $C$ is symmetric and positive semi-definite, and that is what makes everything work.

So the question "which direction has the most variance?" is now the clean optimization problem: maximize $u^\top C u$ subject to $\|u\| = 1$. Set it up with a Lagrange multiplier $\lambda$ for the constraint and differentiate, and the condition for the maximum is startlingly simple:

In words: find the special directions $u$ that $C$ does not rotate at all — it only stretches them by a factor $\lambda$. That is the definition of an eigenvector. The directions of maximum variance are exactly the eigenvectors of the covariance matrix. And the eigenvalue is not a leftover — substitute $Cu = \lambda u$ back into the variance we were maximizing and $u^\top C u = u^\top (\lambda u) = \lambda \|u\|^2 = \lambda$. The eigenvector points along a principal axis; its eigenvalue is the variance of the data along that axis. Rank the eigenvectors by eigenvalue, largest first, and you have the principal components in order of importance. Because $C$ is symmetric, these axes are guaranteed to be mutually perpendicular — the rotated coordinate system PCA promised.

Dimensionality reduction: project onto the top axes

Finding the axes is only half the point. The payoff is throwing most of them away. Keep the top $k$ eigenvectors — the ones with the largest eigenvalues — stack them as the columns of a matrix $W$ (size $d \times k$), and project every data point onto them:

$$Z = X W.$$

In words: measure each point's shadow along each of the top $k$ axes, and use those lengths as its new coordinates. $Z$ is an $n \times k$ matrix, one row of new coordinates per sample. Each original point, which had $d$ coordinates, is now described by just $k$ numbers — its coordinates along the $k$ most important axes. If your 100-feature dataset really lives near a 5-dimensional shape, you can set $k = 5$ and keep almost everything that mattered while discarding 95 columns of near-redundant noise. That is dimensionality reduction: fewer numbers per sample, faster downstream models, and plots you can actually look at when $k = 2$ or $3$.

What you lose is the variance along the axes you dropped — and because you dropped the smallest eigenvalues, that loss is as small as it can possibly be for any $k$-dimensional linear projection. That optimality is the reason PCA, first written down by Pearson in 1901, is still the default first move.

The explained-variance ratio: how much did you keep?

How do you choose $k$, and how do you know the reduction was safe? The eigenvalues answer both, because each one is literally the variance captured by its axis. The fraction of the total variance carried by the $i$-th component is its explained-variance ratio:

$$\text{explained variance}_i = \frac{\lambda_i}{\lambda_1 + \lambda_2 + \dots + \lambda_d}.$$

In words: the share of the total variance that this one axis is responsible for. The denominator, the sum of all eigenvalues, is the total variance of the data; the numerator is the variance along this single axis. In two dimensions this is just $\lambda_1 / (\lambda_1 + \lambda_2)$ for the first component. Add these up over the components you keep and you get the total fraction of variance retained — say, "the top 5 components explain 92% of the variance." That single percentage is how you defend a choice of $k$: plot the cumulative explained variance against the number of components (the scree plot), and pick the point where the curve flattens and each new axis buys you almost nothing.

Seeing it in code

PCA is a handful of lines on top of one linear-algebra call. Here it is on a small correlated cloud, done by hand so every step above is visible:

import numpy as np

# n samples, 2 correlated features
X = np.array([[2.5, 2.4], [0.5, 0.7], [2.2, 2.9], [1.9, 2.2],
              [3.1, 3.0], [2.3, 2.7], [2.0, 1.6], [1.0, 1.1]])

X = X - X.mean(axis=0)              # 1. center: every column now has mean 0
C = (X.T @ X) / len(X)             # 2. covariance matrix  C = (1/n) XᵀX

vals, vecs = np.linalg.eigh(C)     # 3. eigenvalues/vectors (eigh: C is symmetric)
order = np.argsort(vals)[::-1]     #    sort largest eigenvalue first
vals, vecs = vals[order], vecs[:, order]

print("explained variance:", vals / vals.sum())   # e.g. [0.96, 0.04]

W = vecs[:, :1]                    # 4. keep the top k=1 axis
Z = X @ W                          # 5. project: 2-D points → 1 number each

Run it and the first component explains the overwhelming majority of the variance — the two features were nearly a straight line, so one axis is almost enough to describe the whole cloud. Note the call is eigh, not eig: because the covariance matrix is symmetric, the specialized routine is faster and returns real, orthogonal eigenvectors. (Production code usually reaches for the singular value decomposition (SVD) of $X$ instead, which is numerically steadier and gives the same axes without ever forming $C$ — but the eigen-picture above is what it is computing.)

The pitfalls: scale and linearity

PCA is at the mercy of your units. The covariance matrix is built from raw numbers, so a feature measured in grams will have a variance a million times larger than the same feature in kilograms — and PCA will hand its first component almost entirely to whichever feature happens to carry the biggest numbers, regardless of whether it is the most meaningful. When features are in different units (age in years, income in dollars, height in metres), you must standardize first — subtract the mean and divide by the standard deviation of each column — which is equivalent to running PCA on the correlation matrix instead of the covariance matrix. Forgetting this is the single most common way PCA gives a useless answer.

PCA can only find straight lines. Every axis it returns is a linear combination of the original features, so it can only capture linear structure. If your data lies on a curved shape — points wrapped around a spiral, or on the surface of a sphere — no straight axis describes them well, and PCA will report a muddle. The variance along a curled-up direction looks small even when that direction carries all the interesting structure. When the structure is nonlinear, you reach for kernel PCA or manifold-learning methods (t-SNE, UMAP — covered in the clustering and dimensionality reduction article) that are built to unroll curved shapes; PCA is the linear baseline you try first and keep only if it already works.

How this shows up on the job

PCA is the standard preprocessing and compression step. Before feeding high-dimensional features into a model, reducing to the top components cuts training time and can improve results by stripping out low-variance noise directions. In scikit-learn this is PCA(n_components=0.95) — keep however many axes are needed to retain 95% of the variance — and it drops straight into a pipeline before your estimator.

It is the fastest way to see high-dimensional data. Projecting to the top two or three components and plotting is the first thing most people do with a new embedding or dataset — clusters, outliers, and gross structure jump out of a 2-D scatter that were invisible in 100 dimensions. Just remember the axes are abstract combinations of features, not any single measurement.

Watch the explained-variance curve before you trust the reduction. If the top few components explain almost everything, the data really was low-dimensional and PCA is safe. If the explained variance is spread thinly across dozens of components — a flat scree plot — the data genuinely needed those dimensions, and aggressive reduction will quietly throw away signal. The percentage is the honesty check on the whole exercise.

Summary

• Correlated features mean data huddles near a lower-dimensional shape; PCA finds that shape by rotating the axes to line up with the directions of greatest variance
• Center the data first, then the direction of maximum variance is the top eigenvector of the covariance matrix $C = \frac{1}{n} X^\top X$
• Each eigenvalue is the variance along its axis; rank eigenvectors by eigenvalue to get the principal components in order
• Project onto the top $k$ eigenvectors ($Z = XW$) to reduce dimensions with the smallest possible loss; the explained-variance ratio $\lambda_i / \sum_j \lambda_j$ tells you what you kept
• Two traps: PCA is scale-dependent (standardize features in different units) and strictly linear (curved structure needs kernel PCA or manifold methods)

Check your understanding

1. Why must you center the data before computing the covariance matrix? What would the "first principal component" tend to capture if you skipped it?
2. The first eigenvalue of a covariance matrix comes out as $\lambda_1 = 8.0$ and the second as $\lambda_2 = 2.0$. What fraction of the variance does the first principal component explain, and what does keeping only that component throw away?
3. You run PCA on a dataset where one feature is a distance in millimetres and another is a temperature in degrees. Without standardizing, which feature will dominate the first component, and why is that a problem?
4. A colleague's data lies on a tight spiral. They run PCA, keep the top component, and complain it "explains only 40% of the variance." What is going wrong, and what class of method should they reach for instead?