Five hundred columns rarely mean five hundred ideas. Most wide datasets repeat themselves: correlated sensors, near-duplicate measurements, plain noise. Python PCA (principal component analysis) finds the few directions where the real variation lives and lets you drop the rest. This guide works the math with real numbers, builds PCA from scratch in NumPy, then does it properly in scikit-learn, including the two mistakes that silently break it.

“Code is like humor. When you have to explain it, it’s bad.”

Cory House

Last Updated: July 2026 | Tested on: Python 3.14.6, scikit-learn 1.9.0 | Difficulty: Advanced | Reading Time: 16 minutes

Start with the intuition, because the math makes a lot more sense once you can picture it. A dataset with 500 features is almost never 500 times richer than one with 10. Plenty of those columns are correlated, repeat each other, or are just noise. PCA finds the directions in your data that carry the most variation, then rewrites every row using only those directions. Keep the top handful and you have squeezed a wide table into a narrow one with barely any loss.

Here is the everyday picture. Think of a panoramic group photo at a wedding. The camera captures it in full 3D, but a flat print throws away depth and you still recognize every face. PCA does the same trick with numbers: it picks the viewing angle that keeps the most detail, then flattens the data onto it. Lose the depth you do not need, keep the faces you do.

One thing to be clear about up front. PCA is not a clustering algorithm and it is not a supervised model. It never looks at your labels. It is a pure transformation that rewrites your features in a new coordinate system where the first axis captures the most variance, the second captures the next most (at a right angle to the first), and so on down the line. The payoffs are faster training, less overfitting, some noise removal, and the ability to plot a 50-column dataset in 2D so you can actually see it.

Prerequisites

The PCA Pipeline

Read the pipeline top to bottom. You start with your raw table of N features. You standardize it so every column is on the same scale. You build the covariance matrix, which is just a grid of how strongly each feature moves with every other feature. Eigendecomposition of that grid hands you the principal axes: new directions, each ranked by how much variance it captures. You sort them, keep the top K, and project your data onto just those K.

The result is a narrower table that still carries most of the information. That variance ranking is the whole reason you can keep the top 2 or 3 components and still recognize the shape of a 50-column dataset, which is exactly why Python PCA is the go-to step before plotting or training on wide data.

The Math, With Real Numbers

Most Python PCA explanations throw eigenvectors at you in the abstract and hope it sticks. Let us not do that. We will take four students, two exam scores each, and turn the crank by hand. Math and physics scores tend to rise together (a strong student is usually strong in both), so the two columns are highly correlated. That redundancy is exactly what PCA loves to find.

Quick real-life picture before the numbers: think of rating a restaurant on “taste” and on “presentation.” In practice those two scores move together, so a single “overall quality” rating captures almost everything both numbers were telling you. PCA builds that combined rating for you automatically. Now let us watch it happen with real values.

📄 pca_handcalc.py: PCA on 4 students by hand, then checked in code

import numpy as np

# 4 students, 2 exam scores each (math, physics). Strongly correlated.
math_score    = np.array([60.0, 70.0, 80.0, 90.0])
physics_score = np.array([62.0, 68.0, 84.0, 86.0])
X = np.column_stack([math_score, physics_score])

# Step 1: standardize each column (subtract mean, divide by std)
mean = X.mean(axis=0)
std = X.std(axis=0)
X_std = (X - mean) / std
print(f"Mean per column: {np.round(mean, 3)}")
print(f"Std  per column: {np.round(std, 3)}")

# Step 2: covariance matrix of the standardized data
cov = np.cov(X_std, rowvar=False)
print("Covariance (= correlation) matrix:")
print(np.round(cov, 3))

# Step 3: eigenvalues and eigenvectors, largest first
vals, vecs = np.linalg.eigh(cov)
order = np.argsort(vals)[::-1]
vals, vecs = vals[order], vecs[:, order]
print(f"Eigenvalues: {np.round(vals, 3)}")

# Step 4: explained variance ratio
ratio = vals / vals.sum()
print(f"Explained variance ratio: {np.round(ratio, 3)}")
print(f"PC1 alone keeps {ratio[0]:.1%} of the variance")

▶ Output

Mean per column: [75. 75.]
Std  per column: [11.18  10.247]
Covariance (= correlation) matrix:
[[1.333 1.28 ]
 [1.28  1.333]]
Eigenvalues: [2.614 0.053]
Explained variance ratio: [0.98 0.02]
PC1 alone keeps 98.0% of the variance

Walking through the numbers: Both columns center on a mean of 75. After standardizing, the covariance matrix is almost all ones, because math and physics move together so tightly. Now look at the eigenvalues: 2.614 and 0.053. An eigenvalue is just the amount of variance captured along its direction, and they always sum to the total spread in the data (here 2.614 + 0.053 = 2.667). Divide each eigenvalue by that total and you get the explained variance ratio: 0.98 and 0.02. In plain words, one new axis (PC1) holds 98% of everything these two exams told you.

You could replace both scores with a single “overall academic strength” number and lose almost nothing. That single number is the first principal component.

The eigenvector that goes with that big eigenvalue points at roughly equal weight on both subjects, which makes sense: the direction of most variation is “both scores high together versus both low together.” The second eigenvector points the other way, “did better in one subject than the other,” and it barely matters here because almost nobody in this tiny group did. That is PCA in one sitting: find the directions, rank them by variance, keep the ones that count.

PCA From Scratch With NumPy

The four-student example fit in your head. Real data does not, so let us scale the same five steps up to 200 rows and 3 features, where one feature is deliberately built as a near-copy of another. Watch PCA catch the redundancy automatically.

📄 pca_from_scratch.py: the five PCA steps with NumPy

import numpy as np

rng = np.random.default_rng(42)
n = 200
f1 = rng.normal(0, 1, n)
f2 = rng.normal(0, 1, n)
f3 = 2 * f1 + rng.normal(0, 0.1, n)  # correlated with f1
X = np.column_stack([f1, f2, f3])

# Step 1: Standardize
X_std = (X - X.mean(axis=0)) / X.std(axis=0)

# Step 2: Covariance matrix
cov_matrix = np.cov(X_std, rowvar=False)
print("Covariance matrix:")
print(np.round(cov_matrix, 3))

# Step 3: Eigendecomposition
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

# Step 4: Explained variance
explained_var = eigenvalues / eigenvalues.sum()
cumulative = np.cumsum(explained_var)

print(f"\nExplained variance ratios:")
for i, (ev, cum) in enumerate(zip(explained_var, cumulative)):
    print(f"  PC{i+1}: {ev:.4f} ({ev:.1%}) cumulative: {cum:.1%}")

# Step 5: Project onto top 2 components
W = eigenvectors[:, :2]
X_pca = X_std @ W
print(f"\nOriginal shape: {X.shape}")
print(f"Reduced shape:  {X_pca.shape}")
print(f"Information retained: {cumulative[1]:.1%}")

▶ Output

Covariance matrix:
[[ 1.005 -0.071  1.003]
 [-0.071  1.005 -0.069]
 [ 1.003 -0.069  1.005]]

Explained variance ratios:
  PC1: 0.6693 (66.9%) cumulative: 66.9%
  PC2: 0.3302 (33.0%) cumulative: 99.9%
  PC3: 0.0006 (0.1%) cumulative: 100.0%

Original shape: (200, 3)
Reduced shape:  (200, 2)
Information retained: 99.9%

What happened here: Feature 3 was built as roughly two times feature 1 plus a sprinkle of noise, so the two columns are almost perfectly correlated (about 0.998, which is why the covariance matrix shows a 1.003 between them). PCA spotted that on its own. PC3 captures only 0.06% of the variance, which is its way of saying “feature 3 tells me almost nothing feature 1 did not already say.” The top two components together hold 99.9% of the variance, so dropping the third axis costs you essentially nothing.

That is the entire job of Python PCA at scale: find the redundancy you did not know was there and quietly remove it. Your own covariance values may differ by a hair depending on the random draw, but the story (one component near zero) stays the same.

Python PCA with scikit-learn

You will almost never hand-roll Python PCA in production. scikit-learn 1.9.0 (latest stable at the time of writing) gives you PCA in two lines, and it does the eigendecomposition the numerically stable way (via Singular Value Decomposition, or SVD) under the hood. Here is a real test you can feel: the built-in handwritten digits dataset has 1797 images, each an 8 by 8 grid, so 64 features per image. How many of those 64 do we actually need to classify the digit? Let us sweep the component count and watch both the variance kept and the cross-validated accuracy.

📄 sklearn_pca.py: how many components does the digits dataset really need

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_digits
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

digits = load_digits()
X, y = digits.data, digits.target       # 1797 images, 64 features each
X_scaled = StandardScaler().fit_transform(X)

print("Components | Variance | CV Accuracy")
print("-----------|----------|------------")
for nc in [10, 20, 30, 40, 64]:
    pca = PCA(n_components=nc, random_state=0)
    X_pca = pca.fit_transform(X_scaled)
    var = pca.explained_variance_ratio_.sum()
    acc = cross_val_score(
        LogisticRegression(max_iter=5000), X_pca, y, cv=5
    ).mean()
    note = "  <-- sweet spot" if nc == 30 else ""
    print(f"{nc:>10} | {var:>7.1%} | {acc:>10.3f}{note}")

▶ Output

Components | Variance | CV Accuracy
-----------|----------|------------
        10 |   58.9% |      0.850
        20 |   79.3% |      0.905
        30 |   89.3% |      0.911  <-- sweet spot
        40 |   95.1% |      0.914
        64 |  100.0% |      0.919

What happened here: Cutting from 64 features down to 30 keeps 89% of the variance and lands at 0.911 accuracy, versus 0.919 with the full 64. You threw away more than half the columns and lost less than one percentage point of accuracy. That is the trade PCA offers: a much smaller, faster, less noisy feature set for a tiny, often invisible, dip in performance. Notice the diminishing returns too.

Going from 30 to 40 to 64 components barely moves the accuracy needle, because those later components are mostly capturing pixel noise, not the shape of the digit. The popular "keep 95% of variance" rule would pick 40 components here, but 30 is already plenty for this classifier. Always sweep it on your own data rather than trusting a single magic threshold. (These numbers are reproducible: cross_val_score with an unshuffled 5-fold split is deterministic, so you should see exactly these values.)

When PCA Wins and When It Loses

Python PCA is a tool, not a reflex. Reach for it when these are true:

• You have many correlated features. Sensor readings, pixel values, gene expression columns, survey items that all measure the same trait. The more redundancy, the bigger the win.
• You want to plot high-dimensional data. Squash 50 columns to 2 components and you can finally see clusters and outliers on a scatter plot.
• Training is slow or overfitting. Fewer features means faster fits and less room for the model to memorize noise.

And it quietly loses in these cases:

• You need to explain the model. A principal component is a blend of all your original features, so "PC1 went up" does not map to "salary went up." If a regulator or a stakeholder needs a plain answer, keep your original columns and use feature selection instead.
• The structure is non-linear. PCA only finds straight-line directions. Data shaped like a spiral or a swiss roll defeats it. Reach for t-SNE (t-distributed Stochastic Neighbor Embedding) or UMAP (Uniform Manifold Approximation and Projection) there.
• You already have few features. Compressing 8 columns to 6 is rarely worth the loss of interpretability.

Common Mistakes

Mistake 1: Running PCA without standardizing first

This is the big one, and it is easy to prove. PCA hunts for the directions of maximum variance, and variance is measured in whatever units your columns happen to use. Picture a dataset with two columns: age in years (roughly 15 to 45) and salary in rupees (tens of thousands). Salary's raw numbers are a thousand times bigger, so its raw variance dwarfs age's. Skip scaling and PCA will conclude that salary is basically the only thing that varies.

Think of it like judging a singing contest where one mic is cranked to maximum and the other is barely on: the loud one wins every time, no matter who actually sang better. Standardizing turns every mic to the same volume.

❌ The mistake, and the fix, side by side

import numpy as np
from sklearn.decomposition import PCA

rng = np.random.default_rng(0)
n = 500
age = rng.normal(30, 5, n)               # ranges ~ 15 to 45
salary = rng.normal(60000, 15000, n)     # tens of thousands
X = np.column_stack([age, salary])

# WRONG: PCA on raw data. The huge-range salary column hijacks PC1.
pca_raw = PCA().fit(X)
print("No scaling, PC1 loadings [age, salary]:",
      np.round(pca_raw.components_[0], 4))
print("No scaling, explained variance ratio:",
      np.round(pca_raw.explained_variance_ratio_, 4))

# RIGHT: standardize first so both columns get a fair say.
X_std = (X - X.mean(axis=0)) / X.std(axis=0)
pca_std = PCA().fit(X_std)
print("Scaled,    PC1 loadings [age, salary]:",
      np.round(pca_std.components_[0], 4))
print("Scaled,    explained variance ratio:",
      np.round(pca_std.explained_variance_ratio_, 4))

▶ Output

No scaling, PC1 loadings [age, salary]: [0. 1.]
No scaling, explained variance ratio: [1. 0.]
Scaled,    PC1 loadings [age, salary]: [0.7071 0.7071]
Scaled,    explained variance ratio: [0.5026 0.4974]

Why this matters: Without scaling, PC1's loadings are [0, 1]: it points purely along salary and ignores age completely, and it claims 100% of the variance. Age has been erased. After standardizing, both columns weigh in equally ([0.7071, 0.7071]) and the variance splits roughly half and half. The fix is one line, or better, drop a StandardScaler into a Pipeline so you cannot forget it. The order is always StandardScaler, then PCA, then your model.

Mistake 2: Fitting PCA on the test set

Fit your scaler and your PCA on the training data only, then call .transform() on the test data. If you fit PCA on the whole dataset before splitting, information from the test set leaks into your components and your reported accuracy is a lie. The clean way is a Pipeline([("scaler", StandardScaler()), ("pca", PCA(n_components=30)), ("clf", LogisticRegression())]) so the scaler and PCA learn from training folds only, automatically, inside cross-validation.

Practice Exercises

1. Exercise 1: Load the Iris dataset (from sklearn.datasets import load_iris), standardize the 4 features, fit PCA, and print the explained variance ratio. How many components do you need to keep 95% of the variance? Then project to 2 components and confirm the shape is (150, 2).
2. Exercise 2: Write a function components_for_variance(X, threshold) that standardizes X, fits PCA, and returns the smallest number of components whose cumulative explained variance reaches threshold (for example 0.90). Test it on the digits dataset and check it returns roughly 30 for a 0.90 threshold.
3. Exercise 3: Build a full Pipeline of StandardScaler, PCA, and LogisticRegression on the digits dataset. Use GridSearchCV to search pca__n_components over [20, 30, 40] and report which value wins. This combines PCA with the model selection ideas from the train/test split post.

Conclusion

You now know what Python PCA actually does: standardize your features, build the covariance matrix, run eigendecomposition, rank the new axes by explained variance, and keep only the top few. You saw it three ways: by hand on four students, from scratch on redundant NumPy features, and in production on the digits dataset, where 30 components did nearly the same job as all 64. You also picked up the two habits that separate a working PCA from a broken one: always standardize first, and always fit on the training data only.

Next up is hyperparameter tuning, where you learn to search for the settings that make a model perform its best. For the full path from Python basics to deployed models, visit the Python + AI/ML tutorial series home.

Frequently Asked Questions

Interview Questions on PCA

If you can walk through these without peeking, you are ready for this topic in an interview.

Q: Why must you standardize features before running PCA?

PCA maximizes variance, and variance depends on the units each column is measured in. A salary column in rupees will have a raw variance thousands of times larger than an age column in years, so without scaling PCA concludes salary is the only thing that varies and ignores everything else. Standardizing (zero mean, unit variance) puts every feature on the same footing so each gets a fair say in the components. The one exception is when all your features already share the same units and scale, such as pixel intensities.

Q: What does an eigenvalue represent in PCA, and how does it relate to explained variance?

Each eigenvalue is the amount of variance the data has along its matching eigenvector (principal axis). The eigenvectors are the new directions and the eigenvalues rank them by importance. Divide an eigenvalue by the sum of all eigenvalues and you get the explained variance ratio for that component, which is why the ratios always add up to 1.

Q: How is PCA different from feature selection?

Feature selection keeps a subset of your original columns, so you can still point at "age" or "salary" and explain the model. PCA creates brand new features that are linear blends of all the originals, so after PCA no single number maps back to one real-world column. Choose feature selection when interpretability matters and PCA when you mainly want speed, denoising, or a 2D view of wide data.

Q: Scenario: you fit PCA on the entire dataset, then split into train and test. Cross-validated accuracy looks great, but production performance is worse. What went wrong?

Fitting PCA (and the scaler) on the full dataset before splitting leaks information from the test rows into the components, so your reported score is optimistic. The fix is to fit the scaler and PCA on the training folds only and merely .transform() the test data. In practice, wrap StandardScaler, PCA, and the model in a single Pipeline so cross-validation refits each step inside every fold and the leak cannot happen.

Q: Scenario: your PCA keeps 99% of the variance with 2 components, but your classifier's accuracy is poor. What do you check first?

High explained variance does not guarantee the retained directions are the ones that separate your classes, because PCA is unsupervised and never sees the labels. The class-discriminating signal may live in a low-variance direction that you just discarded. Check by keeping more components or by comparing model accuracy with and without PCA, and if the structure is non-linear consider t-SNE, UMAP, or a supervised method like LDA instead.

Q: How do you decide how many components to keep?

Look at the cumulative explained variance ratio and keep enough components to reach roughly 90 to 95 percent, or find the elbow where extra components stop adding much. For visualization, keep 2 or 3. The most reliable approach is to sweep the component count and watch your model's cross-validated score, since the right number depends on your data and task, not on a fixed threshold.

Series: Python + AI/ML Cookbook. Part 5: Machine Learning

Go deeper: scikit-learn documentation covers every edge case of this topic.

Related Posts

Previous: ML: Hierarchical Clustering & DBSCAN

Next: ML: Hyperparameter Tuning (Grid, Random, Bayesian)

Series Home: Python + AI/ML Tutorial Series