Principal Component Analysis (PCA)

The Story That Explains PCA

Principal Component Analysis (PCA) is an unsupervised linear transformation technique. It re-expresses your data in a new coordinate system where the axes (principal components) are ordered by how much variance they explain. The first principal component (PC1) is the direction of maximum spread in your data; PC2 is the direction of maximum spread perpendicular to PC1, and so on.

Why Dimensionality Reduction Matters

Modern datasets routinely have hundreds or thousands of features. More features sounds better — but past a point it backfires in three serious ways:

The Four-Step PCA Algorithm

The Maths — Made Concrete with Numbers

Numerical Example: 2 Features → 1 Principal Component

Let's walk through PCA on a tiny dataset — 5 students with two exam scores. We want to find the single direction (PC1) that best summarises both scores.

Key Formulas

Visual Diagram — What PCA Actually Does

Scree Plot — Choosing How Many Components to Keep

Python Implementation — Step by Step

Step 1: Standardisation and Raw PCA from Scratch

import numpy as np

# ── Toy dataset: 6 students × 4 exam scores ──────────────────
X = np.array([
    [85, 88, 82, 90],
    [70, 72, 68, 74],
    [92, 95, 91, 96],
    [60, 58, 63, 55],
    [78, 80, 76, 82],
    [55, 52, 57, 50],
], dtype=float)

# Step 1 — Standardise (zero mean, unit variance)
mean = X.mean(axis=0)
std  = X.std(axis=0, ddof=1)
X_std = (X - mean) / std

# Step 2 — Covariance matrix
C = np.cov(X_std.T)
print("Covariance matrix shape:", C.shape)
print(np.round(C, 3))

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

# Step 4 — Explained variance ratio
evr = eigenvalues / eigenvalues.sum()
print("Eigenvalues:         ", np.round(eigenvalues, 4))
print("Explained var ratio: ", np.round(evr, 4))
print("Cumulative EVR:      ", np.round(np.cumsum(evr), 4))

# Step 5 — Project onto top-2 PCs
W  = eigenvectors[:, :2]
Z  = X_std @ W
print("Projected data (2 PCs):
", np.round(Z, 3))

Step 2: sklearn PCA — The Production Way

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

iris = load_iris()
X, y = iris.data, iris.target
feature_names = iris.feature_names

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

pca_full = PCA(n_components=None)
pca_full.fit(X_scaled)

print("Explained variance ratio per component:")
for i, ev in enumerate(pca_full.explained_variance_ratio_):
    print(f"  PC{i+1}: {ev:.4f}  ({ev*100:.2f}%)")
print(f"
95% variance threshold: {np.argmax(np.cumsum(pca_full.explained_variance_ratio_) >= 0.95) + 1} components")

pca_2d = PCA(n_components=2)
X_2d   = pca_2d.fit_transform(X_scaled)

print(f"
Original shape:  {X.shape}")
print(f"Reduced shape:   {X_2d.shape}")
print(f"Total variance retained: {pca_2d.explained_variance_ratio_.sum():.4f}")

print("
Component loadings (how much each feature contributes):")
for i, comp in enumerate(pca_2d.components_):
    print(f"  PC{i+1}:")
    for fname, loading in zip(feature_names, comp):
        print(f"    {fname:28s}: {loading:+.4f}")

Step 3: Choosing n_components with Variance Threshold

from sklearn.decomposition import PCA

pca_95 = PCA(n_components=0.95)
X_95   = pca_95.fit_transform(X_scaled)
print(f"Components needed for 95% variance: {pca_95.n_components_}")

pca_mle = PCA(n_components='mle')
X_mle   = pca_mle.fit_transform(X_scaled)
print(f"MLE auto-selected components:       {pca_mle.n_components_}")

X_reconstructed = pca_2d.inverse_transform(X_2d)
mse = ((X_scaled - X_reconstructed) ** 2).mean()
print(f"Reconstruction MSE (2 PCs):         {mse:.6f}")

Complete Real-World Example — Digits Dataset

The digits dataset contains 1,797 handwritten digit images (8×8 pixels = 64 features). We will use PCA to compress 64 dimensions to 2, visualise the digit clusters, and then compress to 30 dimensions and train a classifier.

import numpy as np
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

digits = load_digits()
X, y = digits.data, digits.target
print(f"Dataset: {X.shape[0]} images × {X.shape[1]} pixel features")

scaler   = StandardScaler()
X_scaled = scaler.fit_transform(X)

pca_30  = PCA(n_components=30, random_state=42)
X_30    = pca_30.fit_transform(X_scaled)
print(f"Variance retained (30 PCs): {pca_30.explained_variance_ratio_.sum():.4f}")
print(f"Features: 64 → 30 ({30/64*100:.0f}% of original)")

X_train, X_test, y_train, y_test = train_test_split(
    X_30, y, test_size=0.2, random_state=42)

for name, Xtr, Xte in [
    ("No PCA (64 features)",  *train_test_split(X_scaled, test_size=0.2, random_state=42)),
    ("PCA 30 components",      X_train, X_test)
]:
    lr = LogisticRegression(max_iter=1000, random_state=42)
    lr.fit(Xtr, y_train)
    acc = accuracy_score(y_test, lr.predict(Xte))
    print(f"  {name:25s}: {acc:.4f}")

PCA in a Full Machine-Learning Pipeline

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV, cross_val_score
from sklearn.datasets import load_wine

wine = load_wine()
X, y = wine.data, wine.target

pipe = Pipeline([
    ('scaler', StandardScaler()),
    ('pca',    PCA(n_components=0.95)),
    ('svm',    SVC(kernel='rbf', random_state=42))
])

cv_scores = cross_val_score(pipe, X, y, cv=5, scoring='accuracy')
print(f"CV Accuracy: {cv_scores.mean():.4f} ± {cv_scores.std():.4f}")

param_grid = {
    'pca__n_components': [0.90, 0.95, 0.99, 5, 8],
    'svm__C':            [0.1, 1.0, 10.0],
    'svm__gamma':        ['scale', 'auto'],
}

grid = GridSearchCV(pipe, param_grid, cv=5, scoring='accuracy', n_jobs=-1)
grid.fit(X, y)

print(f"Best CV accuracy: {grid.best_score_:.4f}")
print(f"Best params:      {grid.best_params_}")
print(f"PCs selected:     {grid.best_estimator_.named_steps['pca'].n_components_}")

Incremental PCA — When Data Doesn't Fit in RAM

from sklearn.decomposition import IncrementalPCA
import numpy as np

np.random.seed(42)
X_large = np.random.randn(100_000, 50)

ipca = IncrementalPCA(n_components=10, batch_size=500)
ipca.fit(X_large)

print("Explained variance ratio (top 10):")
print(np.round(ipca.explained_variance_ratio_, 4))
print(f"Total: {ipca.explained_variance_ratio_.sum():.4f}")

X_transformed = ipca.transform(X_large[:1000])
print(f"Transformed shape: {X_transformed.shape}")

PCA Diagram — The Full Data Flow

Key Hyperparameters

When to Use PCA — and When Not To

PCA vs Other Dimensionality Reduction Methods

Golden Rules