🎯 What You'll Learn
- Understand the curse of dimensionality and why high-dimensional data is problematic for ML
- Grasp what principal components are geometrically: orthogonal directions of maximum variance in feature space
- Use
pca.explained_variance_ratio_to choose how many components to retain - Implement PCA with scikit-learn for visualization, preprocessing, and noise reduction
- Interpret loading plots and biplots to understand which original features drive each principal component
1 The Curse of Dimensionality
As the number of features (dimensions) in a dataset grows, several fundamental problems emerge that make machine learning harder:
- Data becomes exponentially sparse: In 1D, 10 points cover a line reasonably. In 10D, you'd need 10¹⁰ points for the same coverage. In practice, your 10,000 training samples become microscopically sparse in 100-dimensional space.
- Distances become meaningless: In high dimensions, the ratio of maximum to minimum pairwise distance approaches 1. All points are approximately the same distance from each other — making distance-based algorithms (KNN, K-Means, SVM with RBF) ineffective.
- Many features are correlated: A dataset of 100 economic indicators doesn't have 100 independent dimensions of variation. Perhaps 5–10 underlying factors (GDP growth, inflation, employment, etc.) drive most of the variation. The rest is redundant.
- Overfitting risk increases: Models have more parameters to overfit. A linear model on 100 features has 101 coefficients — each one can be tuned to noise.
import numpy as np
def distance_concentration_demo(n_dims_list, n_points=1000, n_trials=5):
"""Demonstrate that distances concentrate as dimensionality grows."""
np.random.seed(42)
print(f"{'Dimensions':>12} {'Mean Dist':>10} {'Std Dist':>10} {'Std/Mean (rel.)':>16}")
print("-" * 54)
for d in n_dims_list:
ratios = []
for _ in range(n_trials):
X = np.random.randn(n_points, d)
# Pairwise distances between first 100 points
dists = np.sqrt(((X[:100, None] - X[None, :100]) ** 2).sum(axis=2))
upper = dists[np.triu_indices(100, k=1)]
ratios.append(upper.std() / upper.mean())
mean_ratio = np.mean(ratios)
print(f"{d:>12} {'n/a':>10} {'n/a':>10} {mean_ratio:>16.4f}")
# std/mean decreases toward 0 as d increases → all distances similar
distance_concentration_demo([2, 5, 10, 50, 100, 500])
# As dimensions grow, std/mean → 0 (all pairwise distances are nearly equal)For tabular data with hundreds of correlated features, for image data (a 64×64 image has 4096 dimensions), or for text (vocabulary sizes of 50,000+), you cannot train good models directly on the raw features without first reducing dimensionality. PCA is the standard tool for linear dimensionality reduction.
2 What PCA Does
PCA (Principal Component Analysis) finds a new coordinate system for your data. The new axes — called principal components (PCs) — are chosen so that:
- PC1 is the direction along which the data has the maximum variance
- PC2 is the direction of the second-most variance, and is perpendicular (orthogonal) to PC1
- PC3 is orthogonal to both PC1 and PC2, and has the third-most variance
- ...and so on for all n_features PCs
By projecting data onto the top K PCs (where K << n_features), you retain the dimensions that contain the most information (variance) while discarding dimensions that are mostly noise or redundancy.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
np.random.seed(42)
# Create 2D data with high correlation (captures the concept simply)
n = 200
x1 = np.random.randn(n)
x2 = 0.8 * x1 + 0.2 * np.random.randn(n) # x2 is mostly x1 + noise
X_2d = np.column_stack([x1, x2])
# Scale first
scaler = StandardScaler()
X_2d_sc = scaler.fit_transform(X_2d)
# Apply PCA
pca_2d = PCA(n_components=2)
pca_2d.fit(X_2d_sc)
print("=== PCA on 2D Correlated Data ===")
print(f"PC1 explained variance ratio: {pca_2d.explained_variance_ratio_[0]:.4f} "
f"({pca_2d.explained_variance_ratio_[0]*100:.1f}%)")
print(f"PC2 explained variance ratio: {pca_2d.explained_variance_ratio_[1]:.4f} "
f"({pca_2d.explained_variance_ratio_[1]*100:.1f}%)")
print(f"\nPC1 direction: {pca_2d.components_[0].round(4)}")
print(f"PC2 direction: {pca_2d.components_[1].round(4)}")
print(f"PC1 · PC2 = {np.dot(pca_2d.components_[0], pca_2d.components_[1]):.10f}") # ≈ 0
# Visualize original data + PC directions
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].scatter(X_2d_sc[:, 0], X_2d_sc[:, 1], alpha=0.5, s=15)
origin = np.zeros(2)
scale = 2.0
for i, (color, label) in enumerate(zip(['red', 'blue'], ['PC1', 'PC2'])):
comp = pca_2d.components_[i] * np.sqrt(pca_2d.explained_variance_[i]) * scale
axes[0].annotate('', xy=origin + comp, xytext=origin,
arrowprops=dict(arrowstyle='->', color=color, lw=2))
axes[0].text(*(origin + comp * 1.15), label, color=color, fontsize=12, fontweight='bold')
axes[0].set_title('Original Feature Space\nwith PC Directions')
axes[0].set_aspect('equal')
axes[0].grid(alpha=0.3)
# Projected data
X_proj = pca_2d.transform(X_2d_sc)
axes[1].scatter(X_proj[:, 0], X_proj[:, 1], alpha=0.5, s=15, c='steelblue')
axes[1].set_xlabel('PC1 (max variance)')
axes[1].set_ylabel('PC2 (second variance)')
axes[1].set_title('Projected onto PC Space\nData is now uncorrelated')
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('pca_concept.png', dpi=150)3 The Math Intuition
You don't need to implement PCA from scratch, but understanding the math helps you use it correctly and debug problems. The key idea:
PCA finds the eigenvectors of the covariance matrix of the (scaled) data. The covariance matrix C = (1/n) XᵀX captures how each pair of features varies together. Its eigenvectors are the directions of maximum variance (the principal components), and the corresponding eigenvalues tell you how much variance each direction captures.
Let's make this concrete with the same kind of 2D correlated dataset used above. The scatter below shows 200 points with strong positive correlation, its centroid (the data's mean), and the two principal component directions — PC1 (longer arrow, direction of maximum variance) and PC2 (shorter arrow, orthogonal to PC1) — computed via a closed-form eigendecomposition of the 2×2 covariance matrix. No library needed: for a symmetric 2×2 matrix [[a, b], [b, d]], the eigenvalues have an exact quadratic-formula solution.
But don't just take PCA's word for it. Drag the angle slider to rotate a candidate axis (the dashed gray line) around the centroid, and watch the "variance captured" readout update live — it's simply the variance of the data after projecting every point onto that axis. Try to manually find the angle that maximizes captured variance before clicking "Snap to PC1". You'll find they're exactly the same direction — this is the entire idea behind PCA: PC1 isn't an arbitrary choice, it is by definition the variance-maximizing direction, and rotating to any other angle can only capture less.
Rotate the dashed axis and watch the variance readout — PC1 is whichever direction makes that number largest.
import numpy as np
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
np.random.seed(42)
n, d = 100, 4
X_raw = np.random.randn(n, d)
# Make some features correlated
X_raw[:, 1] = 0.8 * X_raw[:, 0] + 0.2 * np.random.randn(n)
X_raw[:, 3] = -0.7 * X_raw[:, 2] + 0.3 * np.random.randn(n)
# Scale
X = StandardScaler().fit_transform(X_raw)
# ── Manual PCA via eigendecomposition ──
cov_matrix = (X.T @ X) / n # covariance matrix: shape (d, d)
print("Covariance matrix shape:", cov_matrix.shape)
print("Covariance matrix (rounded):")
print(cov_matrix.round(3))
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
# eigh returns in ascending order — sort descending
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
explained_ratio = eigenvalues / eigenvalues.sum()
print(f"\nEigenvalues: {eigenvalues.round(4)}")
print(f"Explained variance ratio: {explained_ratio.round(4)}")
# ── Compare to sklearn PCA ──
pca_sk = PCA(n_components=4)
pca_sk.fit(X)
print(f"\nsklearn EVR: {pca_sk.explained_variance_ratio_.round(4)}")
# Should match (up to sign conventions in eigenvectors)In practice, sklearn's PCA uses Singular Value Decomposition (SVD) of the data matrix, which is numerically more stable than directly computing eigenvectors of the covariance matrix (especially for high-dimensional data). The results are mathematically equivalent. You never need to implement this yourself — just know that eigenvalues ↔ variance and eigenvectors ↔ PC directions.
4 Explained Variance Ratio: How Many Components?
The most important post-fit attribute is pca.explained_variance_ratio_: a 1D array where each entry is the proportion of total variance captured by that component. The cumulative sum tells you how much total variance you retain with the first K components.
The standard tool for visualizing this trade-off is the scree plot: a bar (or line) chart of explained variance per component, with a cumulative line overlaid. Below is a scree plot for the 4-feature correlated dataset from Section 3 (where feature 2 ≈ 0.8×feature 1 + noise, and feature 4 ≈ −0.7×feature 3 + noise). Because of those two correlated pairs, you'd expect roughly two "real" dimensions of variation hiding inside four features — watch where the bars drop off and where the cumulative line bends:
Individual explained variance (bars) and cumulative variance (line) for each principal component, computed from the Section 3 dataset's covariance matrix.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import fetch_california_housing
# A real high-dimensional dataset
housing = fetch_california_housing()
X_housing = housing.data
print(f"Housing dataset: {X_housing.shape}") # (20640, 8)
scaler = StandardScaler()
X_sc = scaler.fit_transform(X_housing)
pca = PCA() # no n_components → compute all
pca.fit(X_sc)
evr = pca.explained_variance_ratio_
cumulative_evr = np.cumsum(evr)
print("\nComponent Explained Cumulative")
print("-" * 40)
for i, (ev, cev) in enumerate(zip(evr, cumulative_evr)):
n_comp_str = f"PC{i+1:2d}"
bar = "█" * int(ev * 50)
print(f"{n_comp_str}: {ev:.4f} ({cev:.4f}) {bar}")
# How many components for 90%, 95%, 99% variance?
for threshold in [0.80, 0.90, 0.95, 0.99]:
n_comp = np.argmax(cumulative_evr >= threshold) + 1
print(f"\n{threshold*100:.0f}% variance → {n_comp} components "
f"(reduced from {X_sc.shape[1]})")
# ── Cumulative explained variance plot ──
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].bar(range(1, len(evr)+1), evr, color='steelblue', alpha=0.7, label='Individual')
axes[0].set_xlabel('Principal Component')
axes[0].set_ylabel('Explained Variance Ratio')
axes[0].set_title('Individual Explained Variance')
axes[0].legend()
axes[1].plot(range(1, len(cumulative_evr)+1), cumulative_evr,
'bo-', markersize=8, linewidth=2)
axes[1].axhline(y=0.95, color='red', linestyle='--', label='95% threshold')
axes[1].axhline(y=0.90, color='orange', linestyle='--', label='90% threshold')
axes[1].set_xlabel('Number of Components')
axes[1].set_ylabel('Cumulative Explained Variance')
axes[1].set_title('Cumulative Explained Variance')
axes[1].set_ylim([0, 1.05])
axes[1].legend()
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('pca_variance.png', dpi=150)The 95% threshold is a common starting point but should be adjusted based on your task. For visualization (2D/3D), you might keep only 2–3 components regardless. For a preprocessing pipeline, cross-validate n_components as a hyperparameter. For noise reduction, even 80–85% retention can work well because the discarded components often contain noise.
5 Full PCA Implementation with Scikit-learn
import numpy as np
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.datasets import load_digits
# High-dimensional dataset: 64 features (8x8 images of digits)
digits = load_digits()
X_digits, y_digits = digits.data, digits.target
print(f"Digits dataset: {X_digits.shape}") # (1797, 64)
# ── Method 1: Specify exact n_components ──
scaler = StandardScaler()
X_sc = scaler.fit_transform(X_digits)
pca_20 = PCA(n_components=20, random_state=42)
X_20 = pca_20.fit_transform(X_sc)
print(f"\nWith 20 components: {pca_20.explained_variance_ratio_.sum()*100:.1f}% variance retained")
print(f"Shape: {X_sc.shape} → {X_20.shape}")
# ── Method 2: Auto-select components for 95% variance ──
pca_95 = PCA(n_components=0.95, random_state=42)
X_95 = pca_95.fit_transform(X_sc)
print(f"\nFor 95% variance: {pca_95.n_components_} components selected")
print(f"Shape: {X_sc.shape} → {X_95.shape}")
# ── Method 3: PCA in a Pipeline (correct way for cross-validation) ──
pipe_raw = Pipeline([
('scaler', StandardScaler()),
('clf', LogisticRegression(max_iter=1000, random_state=42))
])
pipe_pca = Pipeline([
('scaler', StandardScaler()),
('pca', PCA(n_components=0.95, random_state=42)),
('clf', LogisticRegression(max_iter=1000, random_state=42))
])
cv_raw = cross_val_score(pipe_raw, X_digits, y_digits, cv=5, scoring='accuracy')
cv_pca = cross_val_score(pipe_pca, X_digits, y_digits, cv=5, scoring='accuracy')
print(f"\nWithout PCA: {cv_raw.mean():.4f} ± {cv_raw.std():.4f}")
print(f"With PCA: {cv_pca.mean():.4f} ± {cv_pca.std():.4f}")
# Key attributes after fitting
print(f"\nExplained variance ratio (first 10): {pca_95.explained_variance_ratio_[:10].round(4)}")
print(f"Components shape (loadings): {pca_95.components_.shape}")
print(f"Singular values (first 5): {pca_95.singular_values_[:5].round(2)}")
# ── Reconstruct data from PCA (inverse transform) ──
X_reconstructed = pca_95.inverse_transform(X_95)
reconstruction_error = np.mean((X_sc - X_reconstructed) ** 2)
print(f"\nMean reconstruction error (MSE): {reconstruction_error:.6f}")
print(f"This represents {(1 - pca_95.explained_variance_ratio_.sum())*100:.1f}% discarded variance")6 PCA for Visualization
Reducing to 2 or 3 dimensions allows direct visual inspection of high-dimensional data. PCA-based visualization can reveal cluster structure, outliers, and class separation that would be completely invisible in the original feature space.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_digits
digits = load_digits()
X, y = digits.data, digits.target
scaler = StandardScaler()
X_sc = scaler.fit_transform(X)
# ── 2D PCA Visualization ──
pca_2d = PCA(n_components=2, random_state=42)
X_2d = pca_2d.fit_transform(X_sc)
fig, axes = plt.subplots(1, 2, figsize=(15, 6))
colors = plt.cm.tab10(np.linspace(0, 1, 10))
for digit in range(10):
mask = y == digit
axes[0].scatter(X_2d[mask, 0], X_2d[mask, 1],
c=[colors[digit]], alpha=0.6, s=15, label=str(digit))
axes[0].set_xlabel(f'PC1 ({pca_2d.explained_variance_ratio_[0]*100:.1f}% variance)')
axes[0].set_ylabel(f'PC2 ({pca_2d.explained_variance_ratio_[1]*100:.1f}% variance)')
axes[0].set_title('MNIST Digits: 64D → 2D via PCA')
axes[0].legend(title='Digit', bbox_to_anchor=(1.02, 1), loc='upper left', fontsize=8)
# Annotate a few representative points
for digit in range(10):
mask = y == digit
centroid = X_2d[mask].mean(axis=0)
axes[0].annotate(str(digit), centroid, fontsize=10, fontweight='bold',
ha='center', va='center',
bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.7))
# ── 3D PCA ──
pca_3d = PCA(n_components=3, random_state=42)
X_3d = pca_3d.fit_transform(X_sc)
ax3d = fig.add_subplot(122, projection='3d')
for digit in range(10):
mask = y == digit
ax3d.scatter(X_3d[mask, 0], X_3d[mask, 1], X_3d[mask, 2],
c=[colors[digit]], alpha=0.4, s=10, label=str(digit))
ax3d.set_xlabel('PC1')
ax3d.set_ylabel('PC2')
ax3d.set_zlabel('PC3')
ax3d.set_title(f'3D PCA ({sum(pca_3d.explained_variance_ratio_)*100:.1f}% variance)')
plt.tight_layout()
plt.savefig('pca_digits_vis.png', dpi=150)
total_2d = pca_2d.explained_variance_ratio_.sum()
total_3d = pca_3d.explained_variance_ratio_.sum()
print(f"2D PCA: {total_2d*100:.1f}% variance captured")
print(f"3D PCA: {total_3d*100:.1f}% variance captured")7 PCA as Preprocessing: Noise Reduction and Speedup
Beyond visualization, PCA is a powerful preprocessing step. Low-variance components often capture noise rather than signal — discarding them can actually improve model performance, not just speed it up.
import numpy as np
import time
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC
from sklearn.pipeline import Pipeline
from sklearn.model_selection import GridSearchCV, cross_val_score
from sklearn.datasets import load_digits
X_digits, y_digits = load_digits(return_X_y=True)
# ── Time comparison: with vs without PCA ──
pipelines = {
'No PCA (64 dims)': Pipeline([('sc', StandardScaler()), ('svm', SVC(kernel='rbf', C=10))]),
'PCA 20 dims': Pipeline([('sc', StandardScaler()), ('pca', PCA(n_components=20)), ('svm', SVC(kernel='rbf', C=10))]),
'PCA 95% variance': Pipeline([('sc', StandardScaler()), ('pca', PCA(n_components=0.95)), ('svm', SVC(kernel='rbf', C=10))]),
}
for name, pipe in pipelines.items():
t0 = time.time()
scores = cross_val_score(pipe, X_digits, y_digits, cv=5, scoring='accuracy')
elapsed = time.time() - t0
print(f"{name}: accuracy={scores.mean():.4f}±{scores.std():.4f} time={elapsed:.2f}s")
# ── Cross-validating n_components as a hyperparameter ──
pipe_cv = Pipeline([
('sc', StandardScaler()),
('pca', PCA()),
('svm', SVC(kernel='rbf', C=10))
])
param_grid = {'pca__n_components': [5, 10, 20, 30, 40, 50]}
grid = GridSearchCV(pipe_cv, param_grid, cv=5, scoring='accuracy', n_jobs=-1)
grid.fit(X_digits, y_digits)
print(f"\nBest n_components: {grid.best_params_['pca__n_components']}")
print(f"Best CV accuracy: {grid.best_score_:.4f}")
# Results table
for params, mean_sc in zip(grid.cv_results_['params'], grid.cv_results_['mean_test_score']):
print(f" n_components={params['pca__n_components']:3d}: {mean_sc:.4f}")Avoid PCA when original feature interpretability is critical — in medical or legal settings where "feature X caused this prediction" matters. After PCA, you can only say "a combination of all features caused this prediction," which is often unacceptable. Also avoid PCA when features are already sparse (e.g., one-hot encoded categories) — PCA destroys the sparse structure and may increase memory usage. Use it primarily for numerical features with potential correlations.
8 Loading Plots and Biplots
The pca.components_ matrix (shape: n_components × n_features) contains the "loadings" — each row is a principal component, and each value tells you how much that original feature contributes to that component. Positive loading: feature increases → PC increases. Negative: inverse relationship.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris
iris = load_iris()
X_iris = iris.data
feature_names = iris.feature_names
scaler = StandardScaler()
X_sc = scaler.fit_transform(X_iris)
pca_iris = PCA(n_components=4)
pca_iris.fit(X_sc)
X_pca = pca_iris.transform(X_sc)
# ── Loading table ──
print("=== PCA Loadings (components_) ===")
print(f"{'Feature':<30} {'PC1':>8} {'PC2':>8} {'PC3':>8} {'PC4':>8}")
print("-" * 62)
for j, fname in enumerate(feature_names):
row = [pca_iris.components_[i, j] for i in range(4)]
print(f"{fname:<30} " + " ".join(f"{v:>8.4f}" for v in row))
evr = pca_iris.explained_variance_ratio_
print(f"\n{'Explained Variance':<30} " + " ".join(f"{v:>8.4f}" for v in evr))
print(f"{'Cumulative':<30} " + " ".join(f"{v:>8.4f}" for v in np.cumsum(evr)))
# ── Loading bar plots ──
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
for i in range(2):
loadings = pca_iris.components_[i]
colors = ['#e41a1c' if l < 0 else '#377eb8' for l in loadings]
axes[i].barh(feature_names, loadings, color=colors, alpha=0.8)
axes[i].axvline(0, color='black', linewidth=0.8)
axes[i].set_xlabel('Loading value')
axes[i].set_title(f'PC{i+1} Loadings ({evr[i]*100:.1f}% variance)\n'
f'Blue = positive, Red = negative')
axes[i].grid(axis='x', alpha=0.3)
plt.suptitle('Feature Contributions to Principal Components', fontsize=12)
plt.tight_layout()
plt.savefig('pca_loadings.png', dpi=150)
# ── Biplot: PC scatter + feature vectors ──
fig, ax = plt.subplots(figsize=(9, 8))
colors_species = ['#e41a1c', '#377eb8', '#4daf4a']
for sp in range(3):
mask = iris.target == sp
ax.scatter(X_pca[mask, 0], X_pca[mask, 1],
c=colors_species[sp], alpha=0.6, s=40, label=iris.target_names[sp])
# Draw feature loading arrows scaled for visibility
arrow_scale = 2.5
for j, fname in enumerate(feature_names):
ax.annotate('', xy=(pca_iris.components_[0, j] * arrow_scale,
pca_iris.components_[1, j] * arrow_scale),
xytext=(0, 0),
arrowprops=dict(arrowstyle='->', color='darkorange', lw=2))
offset_x = 0.1 * np.sign(pca_iris.components_[0, j])
offset_y = 0.1 * np.sign(pca_iris.components_[1, j])
ax.text(pca_iris.components_[0, j] * arrow_scale + offset_x,
pca_iris.components_[1, j] * arrow_scale + offset_y,
fname.replace(' (cm)', ''), fontsize=9, color='darkorange')
ax.set_xlabel(f'PC1 ({evr[0]*100:.1f}%)')
ax.set_ylabel(f'PC2 ({evr[1]*100:.1f}%)')
ax.set_title('Biplot: PCA Scatter + Feature Loading Vectors')
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('pca_biplot.png', dpi=150)
print("\nPC1 Interpretation: all features load positively and strongly.")
print("PC1 captures overall flower size — petal and sepal dimensions together.")
print("PC2 distinguishes sepal length from sepal width (opposite signs).")Real-World Spotlight: Eigenfaces — PCA on Facial Images
One of the most celebrated applications of PCA is Eigenfaces (Turk & Pentland, 1991): reducing high-dimensional face images (thousands of pixels) to a compact PCA representation that still enables accurate face recognition. Each principal component — an "eigenface" — captures a dimension of variation across the face dataset (lighting, pose, expression, identity features).
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report
from sklearn.datasets import fetch_lfw_people
# Load Labeled Faces in the Wild (LFW) dataset
# min_faces_per_person=70 keeps only subjects with enough samples
try:
lfw = fetch_lfw_people(min_faces_per_person=70, resize=0.4)
X_faces = lfw.data
y_faces = lfw.target
target_names = lfw.target_names
h, w = lfw.images.shape[1:]
n_samples, n_features = X_faces.shape
n_classes = len(target_names)
print(f"LFW Dataset: {n_samples} images, {n_features} features ({h}x{w} pixels)")
print(f"Number of people: {n_classes}")
print(f"Classes: {', '.join(target_names)}")
except Exception:
# Fallback: simulate the key metrics if LFW is unavailable
print("(Simulated LFW-like metrics — real LFW requires internet download)")
n_samples, n_features, h, w, n_classes = 1140, 1850, 50, 37, 7
X_faces = np.random.randn(n_samples, n_features)
y_faces = np.repeat(np.arange(n_classes), n_samples // n_classes)
target_names = [f'Person_{i}' for i in range(n_classes)]
print(f"Simulated: {n_samples} images, {n_features} features, {n_classes} classes")
# ── Split ──
X_tr, X_te, y_tr, y_te = train_test_split(
X_faces, y_faces, test_size=0.25, stratify=y_faces, random_state=42)
# ── Scale ──
scaler = StandardScaler()
X_tr_sc = scaler.fit_transform(X_tr)
X_te_sc = scaler.transform(X_te)
# ── PCA: try different n_components ──
print(f"\n{'n_comp':>8} {'Variance':>10} {'Train Acc':>11} {'Test Acc':>10}")
print("-" * 45)
for n_comp in [30, 50, 100, 150, 200]:
pca = PCA(n_components=n_comp, svd_solver='randomized', whiten=True, random_state=42)
X_tr_pca = pca.fit_transform(X_tr_sc)
X_te_pca = pca.transform(X_te_sc)
clf = SVC(kernel='rbf', C=1000, gamma=0.005)
clf.fit(X_tr_pca, y_tr)
tr_acc = clf.score(X_tr_pca, y_tr)
te_acc = clf.score(X_te_pca, y_te)
total_var = pca.explained_variance_ratio_.sum()
print(f"{n_comp:>8} {total_var:>10.4f} {tr_acc:>11.4f} {te_acc:>10.4f}")
# ── Visualize the first eigenfaces ──
pca_vis = PCA(n_components=12, whiten=True, random_state=42)
pca_vis.fit(X_tr_sc)
fig, axes = plt.subplots(3, 4, figsize=(12, 9))
for i, ax in enumerate(axes.flat):
if i < 12:
if n_features == h * w: # Only if shape is correct
eigenface = pca_vis.components_[i].reshape(h, w)
else:
eigenface = pca_vis.components_[i][:h*w].reshape(h, w) if h*w <= n_features else np.random.randn(h, w)
ax.imshow(eigenface, cmap='gray', interpolation='nearest')
ax.set_title(f'Eigenface {i+1}\n({pca_vis.explained_variance_ratio_[i]*100:.1f}%)', fontsize=9)
ax.axis('off')
plt.suptitle('Top 12 Eigenfaces (Principal Components of Face Images)', fontsize=12)
plt.tight_layout()
plt.savefig('eigenfaces.png', dpi=150)
# ── Financial portfolio example ──
print("\n--- Bonus: PCA on Financial Returns ---")
np.random.seed(42)
n_days, n_stocks = 500, 50
# Simulate 3 underlying market factors driving stock returns
factors = np.random.randn(n_days, 3)
loadings_sim = np.random.randn(n_stocks, 3)
noise_sim = 0.3 * np.random.randn(n_days, n_stocks)
returns = factors @ loadings_sim.T + noise_sim
scaler_fin = StandardScaler()
returns_sc = scaler_fin.fit_transform(returns)
pca_fin = PCA(n_components=0.90)
pca_fin.fit(returns_sc)
print(f"50 stocks → {pca_fin.n_components_} PCA factors for 90% variance")
print(f"Variance per factor: {pca_fin.explained_variance_ratio_[:5].round(4)} ...")The eigenfaces demonstration shows the fundamental compression power of PCA: 1850-dimensional face images can be represented by 150 eigenface coefficients while retaining enough information for highly accurate face recognition. Each eigenface is a "basis image" — a pattern of light and dark regions that captures a dimension of variation across faces. The first eigenfaces capture lighting variation (the most variable aspect of face images); later ones capture progressively finer identity features. Modern deep learning approaches have largely superseded eigenfaces for recognition accuracy, but PCA remains the standard for analysis and compression of image datasets.
✍️ Practice Exercises
- Load the Wine dataset (13 features, 3 classes). Apply PCA after scaling and plot the cumulative explained variance. How many components are needed for 95% variance? Visualize the data in 2D PCA space, color-coded by wine class. Do the three classes separate well?
- Train a
KNeighborsClassifieron the Digits dataset with and without PCA preprocessing (try n_components = 10, 20, 30, 40, 64). Plot test accuracy vs n_components. Where is the sweet spot? - Investigate reconstruction quality: apply PCA with 5, 10, 20, 40, 64 components to a digit image. Use
inverse_transformto reconstruct it and display all versions side-by-side. At what n_components is the reconstructed digit clearly recognisable? - Use
pca.components_to make a loading bar plot for the first 3 components of the Wine dataset. What features drive PC1? Do they make physical/chemical sense?
📚 Primary Source for This Lesson
scikit-learn: Decomposition — PCA — complete API reference including randomized PCA, incremental PCA for large datasets, and kernel PCA for non-linear dimensionality reduction.
For the mathematical foundations, see Jolliffe (2002) Principal Component Analysis (2nd ed.) — the definitive reference. For eigenfaces, see Turk & Pentland (1991) "Eigenfaces for Recognition" in the Journal of Cognitive Neuroscience.