PCA — Linear Algebra Perspective

PCA (Principal Component Analysis) is one of the most important algorithms in data science, and it is pure linear algebra.

Goal: find the orthogonal directions of maximum variance in high-dimensional data, and project data onto the top-k of these directions.

Algorithm:
1. Center data: X_c = X - mean(X)
2. Compute covariance matrix: Σ = X_cᵀ X_c / (n-1)
3. Eigendecompose: Σ = Q Λ Qᵀ
4. Sort eigenvectors by eigenvalue (descending)
5. Project: X_pca = X_c @ Q_k  (keep top k eigenvectors)

Alternatively via SVD (more stable):
1. Center data
2. SVD: X_c = U Σ Vᵀ
3. Principal components = columns of V
4. Projected data = U Σ (or equivalently X_c @ V)

Why it works:
The eigenvectors of the covariance matrix are exactly the directions that capture the most variance. This is the "directions where the data spreads out most."

In ML use cases:
• Dimensionality reduction before training
• Noise reduction
• Data visualization (2D/3D projection)
• Removing multicollinearity before linear models
• Feature compression for memory-constrained settings

1import numpy as np
2import matplotlib.pyplot as plt
3from sklearn.decomposition import PCA
4from sklearn.datasets import load_iris
5
6# ── PCA from scratch ──────────────────────────────────────────────
7np.random.seed(42)
8
9# Generate correlated 3D data
10n = 200
11t  = np.linspace(0, 2*np.pi, n)
12X  = np.column_stack([
13    np.cos(t) + 0.1*np.random.randn(n),
14    np.sin(t) + 0.1*np.random.randn(n),
15    0.1*np.random.randn(n)           # 3rd dim is pure noise
16])
17
18# Step 1: Center
19X_c = X - X.mean(axis=0)
20
21# Step 2: Covariance matrix
22Sigma = X_c.T @ X_c / (n - 1)
23
24# Step 3: Eigendecomposition
25vals, vecs = np.linalg.eigh(Sigma)
26
27# Step 4: Sort descending
28idx  = np.argsort(vals)[::-1]
29vals = vals[idx]
30vecs = vecs[:, idx]
31
32print("PCA from scratch:")
33print(f"  Eigenvalues:  {vals.round(4)}")
34for i, v in enumerate(vals):
35    print(f"  PC{i+1}: {v/vals.sum()*100:.1f}% variance")
36
37# Step 5: Project onto top-2 PCs
38X_pca = X_c @ vecs[:, :2]
39print(f"
40  Original shape: {X.shape}")
41print(f"  PCA shape:      {X_pca.shape}")
42
43# ── Compare with sklearn ──────────────────────────────────────────
44pca = PCA(n_components=2)
45X_sklearn = pca.fit_transform(X)
46print(f"
47  sklearn variance ratio: {pca.explained_variance_ratio_.round(4)}")
48print(f"  matches manual:         {vals[:2]/vals.sum()}")
49
50# ── PCA on Iris dataset ───────────────────────────────────────────
51iris   = load_iris()
52X_iris = iris.data          # (150, 4) — 4 features
53
54pca4 = PCA()
55pca4.fit(X_iris)
56cumvar = np.cumsum(pca4.explained_variance_ratio_)
57
58print("
59Iris PCA (4 features → 2D):")
60for i, (v, cv) in enumerate(zip(pca4.explained_variance_ratio_, cumvar)):
61    print(f"  PC{i+1}: {v*100:.1f}% (cumulative: {cv*100:.1f}%)")
62
63X_2d = pca4.transform(X_iris)[:, :2]
64print(f"
65  Projected to 2D: {X_2d.shape}")
66print(f"  2 PCs capture {cumvar[1]*100:.1f}% of variance")
67
68# ── Connection to SVD ─────────────────────────────────────────────
69X_c_iris = X_iris - X_iris.mean(axis=0)
70U, S, Vt = np.linalg.svd(X_c_iris, full_matrices=False)
71
72print(f"
73SVD principal components match PCA:")
74# Singular values² / (n-1) = eigenvalues of covariance
75svd_var = S**2 / (len(X_iris)-1)
76pca_var = pca4.explained_variance_
77print(f"  SVD eigenvalues: {svd_var.round(4)}")
78print(f"  PCA eigenvalues: {pca_var.round(4)}")
79print(f"  Match: {np.allclose(svd_var, pca_var)}")