Dimensionality Reduction with PCA

Many datasets have dozens or hundreds of features, and many of them say almost the same thing. A student's height and arm span, a flat's area in square feet and its number of rooms, the 64 pixels of a handwritten digit: these features move together. Dimensionality reduction squeezes such data into fewer features while keeping most of the information.

The classic method is principal component analysis (PCA). In this lesson you will see what PCA does, decide how many components to keep, draw a 64-dimensional dataset in two dimensions, and use PCA inside a model pipeline.

The idea: directions of greatest spread

Imagine the heights and arm spans of 200 students plotted as points. They lie along a narrow diagonal band, because arm span is close to height. Describing each student by "position along the band" loses almost nothing; the small sideways spread is mostly noise.

PCA finds these directions automatically. The first principal component is the direction along which the data varies most. The second is the direction of greatest remaining variation at right angles to the first, and so on. Projecting the data onto the first few components gives new features that keep most of the variance (spread).

import numpy as np
from sklearn.decomposition import PCA

rng = np.random.default_rng(42)
height = rng.normal(160, 8, 200)              # cm
arm_span = height + rng.normal(0, 3, 200)     # nearly the same as height
X = np.column_stack([height, arm_span])

pca = PCA().fit(X)
print("Explained variance ratio:", pca.explained_variance_ratio_.round(3))
print("First component direction:", pca.components_[0].round(3))

X_1d = PCA(n_components=1).fit_transform(X)
print("Shape before:", X.shape, " after:", X_1d.shape)

Explained variance ratio: [0.957 0.043]
First component direction: [0.684 0.73 ]
Shape before: (200, 2)  after: (200, 1)

The first component keeps about 96% of the variance, so one number per student is almost as good as two. Its direction, roughly equal weights on height and arm span, means "overall body size". Each component is a weighted mix of the original features, and components_ holds those weights.

PCA looks for the largest variance, so a feature measured in big units would dominate. Unless all features share one unit (like pixels), scale them first with StandardScaler.

How many components to keep?

explained_variance_ratio_ tells you the share of the total variance each component keeps. Adding them up shows how much information a given number of components retains. You can also pass a fraction as n_components, and PCA keeps just enough components to reach it.

import numpy as np
from sklearn.datasets import load_wine
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

X, y = load_wine(return_X_y=True)
X_s = StandardScaler().fit_transform(X)

pca = PCA().fit(X_s)
cumulative = np.cumsum(pca.explained_variance_ratio_)
for k in [1, 2, 3, 5, 8, 13]:
    print(f"{k:>2} components keep {cumulative[k - 1]:.1%} of the variance")

pca95 = PCA(n_components=0.95).fit(X_s)
print("Components needed for 95%:", pca95.n_components_)

 1 components keep 36.2% of the variance
 2 components keep 55.4% of the variance
 3 components keep 66.5% of the variance
 5 components keep 80.2% of the variance
 8 components keep 92.0% of the variance
13 components keep 100.0% of the variance
Components needed for 95%: 10

Keeping all components keeps 100% (nothing is lost, the data is only rotated). For wine, 10 of 13 components are needed for 95%, so its features are not very redundant. Common choices are 90% to 99% for modelling, and 2 or 3 components for plotting.

Seeing high-dimensional data

The digits dataset has 64 pixel features per image. No one can picture 64 dimensions, but PCA can draw the data on a flat map with two components.

import matplotlib.pyplot as plt
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA

X, y = load_digits(return_X_y=True)
pca = PCA(n_components=2)
X_2d = pca.fit_transform(X)
print(f"2 of 64 dimensions keep {pca.explained_variance_ratio_.sum():.1%} of the variance")

plt.figure(figsize=(7, 6))
points = plt.scatter(X_2d[:, 0], X_2d[:, 1], c=y, cmap="tab10", s=8)
plt.colorbar(points, label="digit")
plt.xlabel("First principal component")
plt.ylabel("Second principal component")
plt.title("Handwritten digits after PCA")
plt.show()

2 of 64 dimensions keep 28.5% of the variance

Even with under a third of the variance, groups appear: zeros sit together, and so do fours and sixes, while some digits overlap. PCA draws straight-line projections; for prettier maps of clusters, people use t-SNE (sklearn.manifold.TSNE) or the separate UMAP library. Those methods are for pictures only: their output should not be fed to a model, and distances between far-apart groups on a t-SNE map mean little.

PCA inside a pipeline

PCA can make models faster and sometimes more accurate by removing noise. Like a scaler, it learns from data, so fit it on the training set only. The easiest way is to put it in a pipeline between the scaler and the model.

from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_digits(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.25, random_state=42, stratify=y)

for n in [5, 10, 20, 40]:
    model = make_pipeline(StandardScaler(), PCA(n_components=n, random_state=42),
                          LogisticRegression(max_iter=2000))
    model.fit(X_train, y_train)
    print(f"{n:>2} components -> test accuracy {model.score(X_test, y_test):.3f}")

full = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000)).fit(X_train, y_train)
print(f"all 64 pixels -> test accuracy {full.score(X_test, y_test):.3f}")

 5 components -> test accuracy 0.784
10 components -> test accuracy 0.858
20 components -> test accuracy 0.942
40 components -> test accuracy 0.960
all 64 pixels -> test accuracy 0.978

Twenty components (less than a third of the features) already get within a few points of the full model. On large datasets that saving in time and memory can matter a lot; on this small one, the full model is still slightly better. The number of components is a hyperparameter you can tune with cross-validation.

Strengths and limits

PCA is fast, has no labels to worry about, and is great for plotting, compressing and de-noising. But each component mixes all the original features, so the new features are hard to explain ("component 3 is 0.4 × alcohol minus 0.3 × colour..."). PCA also only finds straight-line structure, and it keeps the directions of largest variance, which are not always the ones that best separate the classes.

Common mistakes

• Not scaling first. A feature in large units takes over the first component. Scale unless all features share a unit.
• Fitting PCA on all the data before splitting. That leaks test information; fit it inside a pipeline on the training data.
• Reading components as original features. A component is a mix of features; check components_ to see the weights.
• Using t-SNE output as model input. t-SNE is for visualisation, and it cannot transform new samples.

Key points

• PCA finds new directions (principal components) ordered by how much variance they capture.
• explained_variance_ratio_ shows each component's share; n_components=0.95 keeps enough for 95%.
• Scale features before PCA unless they share one unit.
• Two or three components let you plot high-dimensional data; t-SNE and UMAP make nicer cluster maps.
• Put PCA in a pipeline so it is fitted on training data only, and tune the number of components.

from sklearn.datasets import load_breast_cancer
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_breast_cancer(return_X_y=True)
pca = make_pipeline(StandardScaler(), PCA(n_components=2)).fit(X)
print(f"Variance kept by 2 components: {pca[-1].explained_variance_ratio_.sum():.1%}")

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.25, random_state=42, stratify=y)
two = make_pipeline(StandardScaler(), PCA(n_components=2), LogisticRegression()).fit(X_train, y_train)
full = make_pipeline(StandardScaler(), LogisticRegression()).fit(X_train, y_train)
print(f"2 components: {two.score(X_test, y_test):.3f}")
print(f"30 features:  {full.score(X_test, y_test):.3f}")

Variance kept by 2 components: 63.2%
2 components: 0.944
30 features:  0.986