Principal Component Analysis Explained!

Imagine a dataset with a column for height in centimetres and another for height in inches. You have two columns, but only one piece of information. Any model you train has to spend effort discovering that they're the same thing.

Real data is full of softer versions of this. House size and number of rooms. Ad spend and website visits. Pixel 1,000 and pixel 1,001 of a photo. Features move together, so a dataset with 50 columns might really only have 5 or 6 things going on.

PCA finds those few things. It looks for the directions in which the data actually varies, sorts them from most to least important, and lets you describe each row with just the first few. That's useful for compressing data, removing redundant features, plotting high-dimensional data in 2D, and speeding up whatever model comes next.

This post builds PCA from scratch in NumPy, checks it against scikit-learn, and then looks inside the result to see what it found.

The idea: rotate the axes

Picture a scatter plot of two features that move together, so the points form a long, thin, tilted cloud. The usual x and y axes are a bad description of this cloud: to say where a point is, you need both numbers, and they mostly repeat each other.

Now rotate the axes so the first one runs along the length of the cloud and the second runs across it. Suddenly the first coordinate tells you almost everything, and the second is just a small wobble. Drop the second axis and you've halved the data while losing very little.

That's PCA. In more dimensions it finds a whole set of new axes, called principal components, with three properties:

1. The first component points in the direction where the data varies the most.
2. Each next component points in the direction of the most remaining variation, at right angles to all the earlier ones.
3. The components are sorted by how much of the data's variance they capture.

Keep the first k and you've gone from n features to k, losing as little variance as possible along the way.

Step 1: Center the data

self.mean = np.mean(X, axis=0)
X_centered = X - self.mean

PCA measures spread around the average, so we subtract each feature's mean first. Skip this and the first component tends to point from the origin toward wherever the data happens to sit, instead of along its spread.

The mean is saved on the object, because new data has to be centered with the training mean, not its own.

Step 2: SVD does the heavy lifting

U, S, Vt = np.linalg.svd(X_centered, full_matrices=False)

The Singular Value Decomposition splits the centered data matrix (m samples × n features) into three pieces:

full_matrices=False tells NumPy to skip the parts of U we'd never use. With 500 samples that's the difference between a 500 × 5 matrix and a 500 × 500 one.

Why SVD, and where the "eigen" comes from

The textbook recipe for PCA is different: build the covariance matrix and find its eigenvectors. Plug the SVD into the covariance formula and you can see both routes land in the same place:

That's exactly the form of an eigendecomposition. So the rows of Vt are the eigenvectors of the covariance matrix, and its eigenvalues are S² / (m − 1).

Step 3: How much does each component explain?

m = X.shape[0]
eigenvalues = (S ** 2) / (m - 1)
total_variance = np.sum(eigenvalues)
self.explained_variance_ratio_ = eigenvalues[:self.n_components] / total_variance

Each eigenvalue is the variance of the data along that component. Divide by the total and you get the share of all the variance each component captures. This number is how you decide how many components to keep.

Here's a neat fact that makes this work: the eigenvalues add up to exactly the sum of the original features' variances. Rotating the axes doesn't create or destroy any variance. It only moves it around, piling as much as possible onto the first few axes.

Step 4: Keep the top k, then project

self.components_ = Vt[:self.n_components, :]

SVD already sorts the components by importance, so the first n_components rows are the best directions. To describe a sample in the new coordinates, center it and take the dot product with each component:

def transform(self, X):
    X_centered = X - self.mean
    return X_centered @ self.components_.T

Each sample goes in as n numbers and comes out as k. fit_transform just calls fit and then transform on the same data.

That's the whole algorithm. Here's the class in one piece:

import numpy as np

class CustomPCA:
    def __init__(self, n_components):
        self.n_components = n_components

    def fit(self, X):
        self.mean = np.mean(X, axis=0)
        X_centered = X - self.mean
        U, S, Vt = np.linalg.svd(X_centered, full_matrices=False)

        eigenvalues = (S ** 2) / (X.shape[0] - 1)
        self.explained_variance_ratio_ = eigenvalues[:self.n_components] / eigenvalues.sum()
        self.components_ = Vt[:self.n_components, :]
        return self

    def transform(self, X):
        return (X - self.mean) @ self.components_.T

    def fit_transform(self, X):
        return self.fit(X).transform(X)

Does it find the hidden structure?

To test it I built a dataset where I know the answer: 500 samples and 5 features, made from only 3 independent sources.

np.random.seed(42)
feature1 = np.random.randn(500)
feature2 = feature1 * 2.5 + np.random.randn(500) * 0.2
feature3 = np.random.randn(500)
feature4 = feature3 * -1.2 + np.random.randn(500) * 0.5
feature5 = np.random.randn(500)
X = np.column_stack((feature1, feature2, feature3, feature4, feature5))

pca = CustomPCA(n_components=2)
X_reduced = pca.fit_transform(X)

Original shape: (500, 5)
Reduced shape: (500, 2)
Explained Variance Ratio by top 2 components: [0.65502027 0.24919894]
Total variance captured: 90.42%

Five columns down to two, and 90.4% of the variance survives. Fit all five components and the full picture looks like this:

The variances add up to 10.63, which is exactly the total of the five original feature variances. Nothing was lost by rotating; it was just reshuffled.

The drop after component 3 is the giveaway. Three components carry 99.1% of everything, and the last two carry almost nothing. PCA has worked out that this 5-column dataset really has 3 things going on, which is exactly how it was built. The last two components are just the small noise added to feature2 and feature4.

Looking inside the components

The numbers in components_ say how much each original feature contributes to each new axis:

PCA was never told about the groups, and it found all three:

• Component 1 is group A. It's almost entirely feature1 and feature2, with the same sign because they rise and fall together.
• Component 2 is group B. feature3 and feature4 get opposite signs, because feature4 was built with a negative multiplier.
• Component 3 is group C, which is just feature5 on its own.

One more property is worth seeing: the two new columns in X_reduced have a correlation of essentially zero (about 10⁻¹⁷). The original features were tangled together; the components never are. That's what "at right angles" buys you.

What you lose when you drop components

To see the cost of keeping only 2 components, map the reduced data back into the original 5 features. My class doesn't include an inverse_transform, but it's one line:

X_rebuilt = X_reduced @ pca.components_ + pca.mean

Compare X_rebuilt with X and the squared error is 9.58% of the data's total variance. That's exactly 100% − 90.42%. The "explained variance ratio" isn't an abstract score: it's literally how much of the data you can rebuild from the components you kept.

Checked against scikit-learn

On the demo data, CustomPCA(2) and sklearn.decomposition.PCA(2) agree:

Things that will bite you

Scale decides everything

PCA follows variance, so features with bigger numbers win. In the demo, feature2 has a variance of about 6 while the others sit between 0.9 and 1.7. That's why feature2 gets the biggest weight (−0.93) in component 1. It isn't more important; it's just bigger.

If your features use different units, say metres and dollars, the one with larger numbers will dominate for no good reason. Standardize each feature to mean 0 and standard deviation 1 before PCA, unless the raw scales genuinely mean something (pixel brightness, for example, already shares one scale).

Choosing how many components to keep

A common rule: keep the smallest k whose explained variance adds up to 90–95%. Another is to plot the shares and look for the elbow where they suddenly drop, like the fall from 8.7% to 0.9% above.

Since fit computes every eigenvalue anyway, fit once with n_components equal to the number of features to see the whole picture, then decide.

Big data needs a different SVD

np.linalg.svd decomposes the entire matrix, even though you usually only want the first few components. For large datasets, a randomized SVD that computes only the top k (such as sklearn.utils.extmath.randomized_svd) is far faster. I used exactly that in Eigenfaces From Scratch, where the data was 14,027 face images with 4,096 pixels each.

Takeaways

1. PCA rotates your axes so the first few carry as much of the data's variance as possible.
2. The algorithm is: center, SVD, square the singular values, keep the top rows of Vᵀ.
3. The rows of Vᵀ are the eigenvectors of the covariance matrix. SVD just gets there more accurately.
4. Explained variance tells you how much of the data you can rebuild from the components you keep.
5. Standardize first if your features use different units, and never read meaning into a component's sign.

Written by Abhi Sharma

Data and AI @ Google | I love to write, explain and build things.