Principal component analysis

Principal component analysis (PCA) finds the directions in which data varies most and re-expresses each point using only the top few of those directions. It compresses many correlated features into a handful of new ones, which makes data easier to visualise, faster to model, and less noisy.

The idea

Imagine a cloud of points shaped like a tilted cigar in 2D. Almost all the spread is along the cigar's long axis; very little is across it. If you rotate your axes to line up with the cigar, the first new axis captures nearly everything and the second can be dropped with little loss. PCA does that rotation in any number of dimensions.

The new axes are the principal components. They are:

• Ordered by how much variance they capture (first is largest).
• Perpendicular to each other, so the new features are uncorrelated.
• Linear combinations of the original features.

The maths, briefly

1. Centre the data (subtract each feature's mean). Scale too, if features have different units.
2. Compute the covariance matrix $\Sigma = \frac{1}{m} X^T X$.
3. Find its eigenvectors and eigenvalues. Each eigenvector is a component; its eigenvalue is the variance along it.
4. Keep the top $d$ eigenvectors as columns of $W$ and project: $Z = X W$.

import numpy as np

rng = np.random.default_rng(0)
t = rng.normal(size=100)
X = np.column_stack([t, 0.5 * t + rng.normal(0, 0.1, 100)])   # two highly correlated features
X = X - X.mean(axis=0)

cov = X.T @ X / len(X)
eigenvalues, eigenvectors = np.linalg.eigh(cov)
order = eigenvalues.argsort()[::-1]
print("variance per component:", eigenvalues[order].round(4))
print("first component direction:", eigenvectors[:, order[0]].round(3))

The first component captures almost all the variance, because the second feature is nearly a copy of the first.

With scikit-learn

from sklearn.datasets import load_breast_cancer
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

X, y = load_breast_cancer(return_X_y=True)
X_scaled = StandardScaler().fit_transform(X)

pca = PCA().fit(X_scaled)
ratio = pca.explained_variance_ratio_
print("first 5 components explain:", ratio[:5].round(3))
print("cumulative:", ratio.cumsum()[:10].round(3))

Thirty features, and the first two components already carry over 60 percent of the variance; the first ten carry over 95 percent.

Choosing the number of components

• Pass a fraction and PCA keeps enough components to reach it: PCA(n_components=0.95).
• Or look at the cumulative explained variance and pick the point of diminishing returns.
• For plotting, use 2 or 3.

from sklearn.datasets import load_breast_cancer
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

X, y = load_breast_cancer(return_X_y=True)
X_scaled = StandardScaler().fit_transform(X)
pca = PCA(n_components=0.95).fit(X_scaled)
print("components needed for 95%:", pca.n_components_)
Z = pca.transform(X_scaled)
print("new shape:", Z.shape)

Visualising high dimensional data

Project to two components and plot, colouring by label. Even though PCA never saw the labels, the classes often separate.

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

X, y = load_breast_cancer(return_X_y=True)
Z = PCA(n_components=2).fit_transform(StandardScaler().fit_transform(X))
for label, name in ((0, "malignant"), (1, "benign")):
    centre = Z[y == label].mean(axis=0)
    print(f"{name:10s} centre in PC space: {centre.round(2)}")

The two class centres sit far apart along the first component.

PCA before a model

Fewer, uncorrelated features can make distance based and linear models faster and sometimes more accurate. Put it in a pipeline so it is fit on training data only.

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

X, y = load_breast_cancer(return_X_y=True)
for n in (2, 5, 10, 30):
    pipe = make_pipeline(StandardScaler(), PCA(n_components=n), LogisticRegression(max_iter=1000))
    print(f"{n:2d} components: cv accuracy={cross_val_score(pipe, X, y, cv=5).mean():.3f}")

Reading the components

pca.components_ holds one row per component with a weight per original feature. Large weights show which original features the component mixes.

Limits

• PCA is linear. Curved manifolds need t-SNE, UMAP or autoencoders (those are for visualisation; they do not give a reusable transform in the same way).
• Components are hard to explain to non technical stakeholders ("0.3 of radius minus 0.2 of texture").
• Variance is not always the same as usefulness. A low variance direction can still be the one that separates the classes.

Practice

1. Run PCA on the iris data (scaled) and report how much variance the first two components explain.
2. Print the weights of the first component on the breast cancer features and name the three original features it relies on most.
3. Compare a k-NN classifier on the raw scaled breast cancer features against one on the top 5 principal components.