10 PCA — First Principles¶

Goal. Implement PCA from scratch via SVD, verify against sklearn, and explore projection geometry, reconstruction, and failure cases.

Prerequisites. Linear algebra (eigenvalues, SVD), matrix calculus.

Theory. See theory.md for derivations of variance maximization, reconstruction optimality, and the SVD connection.

In [1]:

import random

import matplotlib.pyplot as plt
import numpy as np

%load_ext autoreload
%autoreload 2

SEED = 42
random.seed(SEED)
rng = np.random.default_rng(SEED)

import random import matplotlib.pyplot as plt import numpy as np %load_ext autoreload %autoreload 2 SEED = 42 random.seed(SEED) rng = np.random.default_rng(SEED)

1. Problem Setup — WHY¶

Generate a 2D correlated cloud. If we must describe it with one axis, which direction preserves the most spread?

In [2]:

# 2D correlated data
n_2d = 200
angle = np.pi / 5
R = np.array([[np.cos(angle), -np.sin(angle)],
              [np.sin(angle),  np.cos(angle)]])
raw = rng.normal(size=(n_2d, 2)) * [3.0, 0.5]  # stretched ellipse
X_2d = raw @ R.T

fig, ax = plt.subplots(figsize=(6, 5))
ax.scatter(X_2d[:, 0], X_2d[:, 1], alpha=0.4, s=20, edgecolor='k', linewidth=0.3)
ax.set_aspect('equal')
ax.set_xlabel('x₁')
ax.set_ylabel('x₂')
ax.set_title('2D correlated data — elongated cloud')
plt.show()

# 2D correlated data n_2d = 200 angle = np.pi / 5 R = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]]) raw = rng.normal(size=(n_2d, 2)) * [3.0, 0.5] # stretched ellipse X_2d = raw @ R.T fig, ax = plt.subplots(figsize=(6, 5)) ax.scatter(X_2d[:, 0], X_2d[:, 1], alpha=0.4, s=20, edgecolor='k', linewidth=0.3) ax.set_aspect('equal') ax.set_xlabel('x₁') ax.set_ylabel('x₂') ax.set_title('2D correlated data — elongated cloud') plt.show()

2. Mathematical Core — WHAT¶

Rotate a line through the cloud and measure "shadow spread" (variance of projection). The angle that maximizes shadow spread also minimizes reconstruction error.

In [3]:

# Shadow experiment: variance as a function of projection angle
Xc_2d = X_2d - X_2d.mean(axis=0)
angles = np.linspace(0, np.pi, 200)
variances = []
recon_errors = []
for a in angles:
    v = np.array([np.cos(a), np.sin(a)])
    proj_scores = Xc_2d @ v
    variances.append(np.var(proj_scores))
    recon = np.outer(proj_scores, v)
    recon_errors.append(np.mean(np.sum((Xc_2d - recon) ** 2, axis=1)))

best_angle = angles[np.argmax(variances)]
print(f'Max variance at angle = {np.degrees(best_angle):.1f}°')
print(f'Min recon error at    = {angles[np.argmin(recon_errors)] * 180 / np.pi:.1f}°')

# Shadow experiment: variance as a function of projection angle Xc_2d = X_2d - X_2d.mean(axis=0) angles = np.linspace(0, np.pi, 200) variances = [] recon_errors = [] for a in angles: v = np.array([np.cos(a), np.sin(a)]) proj_scores = Xc_2d @ v variances.append(np.var(proj_scores)) recon = np.outer(proj_scores, v) recon_errors.append(np.mean(np.sum((Xc_2d - recon) ** 2, axis=1))) best_angle = angles[np.argmax(variances)] print(f'Max variance at angle = {np.degrees(best_angle):.1f}°') print(f'Min recon error at = {angles[np.argmin(recon_errors)] * 180 / np.pi:.1f}°')

Max variance at angle = 36.2°
Min recon error at    = 36.2°

In [4]:

# 3-panel plot: projections at different angles
demo_angles = [0, best_angle, best_angle + np.pi / 2]
fig, axes = plt.subplots(1, 3, figsize=(14, 4.5))
for ax, a in zip(axes, demo_angles):
    v = np.array([np.cos(a), np.sin(a)])
    proj = Xc_2d @ v
    recon = np.outer(proj, v)
    ax.scatter(Xc_2d[:, 0], Xc_2d[:, 1], alpha=0.3, s=10)
    ax.quiver(0, 0, v[0] * 4, v[1] * 4, scale=1, scale_units='xy',
              angles='xy', color='crimson', width=0.008)
    for i in range(0, n_2d, 5):
        ax.plot([Xc_2d[i, 0], recon[i, 0]], [Xc_2d[i, 1], recon[i, 1]],
                color='gray', alpha=0.2, lw=0.5)
    ax.set_aspect('equal')
    ax.set_xlabel('x₁')
    ax.set_title(f'θ = {np.degrees(a):.0f}° | var = {np.var(proj):.2f}')
axes[0].set_ylabel('x₂')
plt.suptitle('Projection angle determines shadow spread and reconstruction error')
plt.tight_layout()
plt.show()

# 3-panel plot: projections at different angles demo_angles = [0, best_angle, best_angle + np.pi / 2] fig, axes = plt.subplots(1, 3, figsize=(14, 4.5)) for ax, a in zip(axes, demo_angles): v = np.array([np.cos(a), np.sin(a)]) proj = Xc_2d @ v recon = np.outer(proj, v) ax.scatter(Xc_2d[:, 0], Xc_2d[:, 1], alpha=0.3, s=10) ax.quiver(0, 0, v[0] * 4, v[1] * 4, scale=1, scale_units='xy', angles='xy', color='crimson', width=0.008) for i in range(0, n_2d, 5): ax.plot([Xc_2d[i, 0], recon[i, 0]], [Xc_2d[i, 1], recon[i, 1]], color='gray', alpha=0.2, lw=0.5) ax.set_aspect('equal') ax.set_xlabel('x₁') ax.set_title(f'θ = {np.degrees(a):.0f}° | var = {np.var(proj):.2f}') axes[0].set_ylabel('x₂') plt.suptitle('Projection angle determines shadow spread and reconstruction error') plt.tight_layout() plt.show()

In [5]:

fig, ax1 = plt.subplots(figsize=(7, 4))
ax1.plot(np.degrees(angles), variances, color='steelblue', label='projection variance')
ax1.set_xlabel('angle (degrees)')
ax1.set_ylabel('variance', color='steelblue')
ax2 = ax1.twinx()
ax2.plot(np.degrees(angles), recon_errors, color='crimson', label='recon error')
ax2.set_ylabel('mean recon error', color='crimson')
ax1.axvline(np.degrees(best_angle), color='black', ls=':', label=f'best = {np.degrees(best_angle):.0f}°')
ax1.legend(loc='upper left', fontsize=8)
ax2.legend(loc='upper right', fontsize=8)
ax1.set_title('Max variance ⟺ Min reconstruction error')
plt.tight_layout()
plt.show()

fig, ax1 = plt.subplots(figsize=(7, 4)) ax1.plot(np.degrees(angles), variances, color='steelblue', label='projection variance') ax1.set_xlabel('angle (degrees)') ax1.set_ylabel('variance', color='steelblue') ax2 = ax1.twinx() ax2.plot(np.degrees(angles), recon_errors, color='crimson', label='recon error') ax2.set_ylabel('mean recon error', color='crimson') ax1.axvline(np.degrees(best_angle), color='black', ls=':', label=f'best = {np.degrees(best_angle):.0f}°') ax1.legend(loc='upper left', fontsize=8) ax2.legend(loc='upper right', fontsize=8) ax1.set_title('Max variance ⟺ Min reconstruction error') plt.tight_layout() plt.show()

3. Solution Method — HOW¶

Three solver routes that should agree:

Route	What it computes	When to use
Covariance eigh	Eigendecomposition of $S = X_c^TX_c/n$	Small $d$
SVD	Thin SVD of $X_c$	General, numerically best
Power iteration	Dominant eigenvector only	Streaming, very large $d$

In [6]:

# 3D data for solver comparison
n_3d, d_3d = 300, 3
raw_3d = rng.normal(size=(n_3d, d_3d)) * [4.0, 1.5, 0.3]
R3 = np.linalg.qr(rng.normal(size=(d_3d, d_3d)))[0]
X_3d = raw_3d @ R3.T
Xc_3d = X_3d - X_3d.mean(axis=0)

# Route 1: covariance eigh
S = Xc_3d.T @ Xc_3d / n_3d
eigvals_cov, eigvecs_cov = np.linalg.eigh(S)
idx = eigvals_cov.argsort()[::-1]
eigvals_cov = eigvals_cov[idx]
V_cov = eigvecs_cov[:, idx]

# Route 2: SVD
U, sigma, Vt = np.linalg.svd(Xc_3d, full_matrices=False)
V_svd = Vt.T
eigvals_svd = sigma ** 2 / n_3d

# Route 3: power iteration (first component only)
def power_iteration(A, n_iter=100, seed=42):
    rng_pw = np.random.default_rng(seed)
    v = rng_pw.normal(size=A.shape[1])
    v /= np.linalg.norm(v)
    for _ in range(n_iter):
        v_new = A @ v
        v = v_new / np.linalg.norm(v_new)
    lam = float(v @ A @ v)
    return lam, v

lam_pw, v_pw = power_iteration(S)

# Check agreement (sign-invariant)
def subspace_alignment(u, v):
    """Absolute cosine between unit vectors."""
    return abs(float(u @ v))

print('Eigenvalue comparison:')
print(f'  eigh:  {eigvals_cov}')
print(f'  SVD:   {eigvals_svd}')
print(f'  power: {lam_pw:.6f} (1st component only)')
assert np.allclose(eigvals_cov, eigvals_svd, atol=1e-10)
assert abs(lam_pw - eigvals_cov[0]) < 1e-8

for j in range(d_3d):
    align = subspace_alignment(V_cov[:, j], V_svd[:, j])
    print(f'  PC{j+1} alignment (eigh vs SVD) = {align:.10f}')
    assert align > 1 - 1e-10
print(f'  PC1 alignment (eigh vs power) = {subspace_alignment(V_cov[:, 0], v_pw):.10f}')

# 3D data for solver comparison n_3d, d_3d = 300, 3 raw_3d = rng.normal(size=(n_3d, d_3d)) * [4.0, 1.5, 0.3] R3 = np.linalg.qr(rng.normal(size=(d_3d, d_3d)))[0] X_3d = raw_3d @ R3.T Xc_3d = X_3d - X_3d.mean(axis=0) # Route 1: covariance eigh S = Xc_3d.T @ Xc_3d / n_3d eigvals_cov, eigvecs_cov = np.linalg.eigh(S) idx = eigvals_cov.argsort()[::-1] eigvals_cov = eigvals_cov[idx] V_cov = eigvecs_cov[:, idx] # Route 2: SVD U, sigma, Vt = np.linalg.svd(Xc_3d, full_matrices=False) V_svd = Vt.T eigvals_svd = sigma ** 2 / n_3d # Route 3: power iteration (first component only) def power_iteration(A, n_iter=100, seed=42): rng_pw = np.random.default_rng(seed) v = rng_pw.normal(size=A.shape[1]) v /= np.linalg.norm(v) for _ in range(n_iter): v_new = A @ v v = v_new / np.linalg.norm(v_new) lam = float(v @ A @ v) return lam, v lam_pw, v_pw = power_iteration(S) # Check agreement (sign-invariant) def subspace_alignment(u, v): """Absolute cosine between unit vectors.""" return abs(float(u @ v)) print('Eigenvalue comparison:') print(f' eigh: {eigvals_cov}') print(f' SVD: {eigvals_svd}') print(f' power: {lam_pw:.6f} (1st component only)') assert np.allclose(eigvals_cov, eigvals_svd, atol=1e-10) assert abs(lam_pw - eigvals_cov[0]) < 1e-8 for j in range(d_3d): align = subspace_alignment(V_cov[:, j], V_svd[:, j]) print(f' PC{j+1} alignment (eigh vs SVD) = {align:.10f}') assert align > 1 - 1e-10 print(f' PC1 alignment (eigh vs power) = {subspace_alignment(V_cov[:, 0], v_pw):.10f}')

Eigenvalue comparison:
  eigh:  [15.24724587  2.22246661  0.09424772]
  SVD:   [15.24724587  2.22246661  0.09424772]
  power: 15.247246 (1st component only)
  PC1 alignment (eigh vs SVD) = 1.0000000000
  PC2 alignment (eigh vs SVD) = 1.0000000000
  PC3 alignment (eigh vs SVD) = 1.0000000000
  PC1 alignment (eigh vs power) = 1.0000000000

4. Implementation — BUILD¶

In [7]:

class PCAScratch:
    """Principal Component Analysis via SVD.

    Attributes after fit:
        components_:               (n_components, n_features) — principal directions
        explained_variance_:       (n_components,)
        explained_variance_ratio_: (n_components,)
        mean_:                     (n_features,)
    """

    def __init__(self, n_components, standardize=False):
        self.n_components = n_components
        self.standardize = standardize

    def fit(self, X):
        X = np.asarray(X, dtype=float)
        n, d = X.shape
        self.mean_ = X.mean(axis=0)
        Xc = X - self.mean_
        if self.standardize:
            self.scale_ = Xc.std(axis=0, ddof=0)
            self.scale_[self.scale_ == 0] = 1.0
            Xc = Xc / self.scale_
        else:
            self.scale_ = None
        U, sigma, Vt = np.linalg.svd(Xc, full_matrices=False)
        k = self.n_components
        self.components_ = Vt[:k]  # (k, d)
        self.explained_variance_ = (sigma[:k] ** 2) / n
        total_var = np.sum(sigma ** 2) / n
        self.explained_variance_ratio_ = self.explained_variance_ / total_var
        self._n_samples = n
        return self

    def transform(self, X):
        Xc = np.asarray(X, dtype=float) - self.mean_
        if self.scale_ is not None:
            Xc = Xc / self.scale_
        return Xc @ self.components_.T

    def inverse_transform(self, Z):
        Xc_hat = np.asarray(Z) @ self.components_
        if self.scale_ is not None:
            Xc_hat = Xc_hat * self.scale_
        return Xc_hat + self.mean_

    def fit_transform(self, X):
        return self.fit(X).transform(X)

class PCAScratch: """Principal Component Analysis via SVD. Attributes after fit: components_: (n_components, n_features) — principal directions explained_variance_: (n_components,) explained_variance_ratio_: (n_components,) mean_: (n_features,) """ def __init__(self, n_components, standardize=False): self.n_components = n_components self.standardize = standardize def fit(self, X): X = np.asarray(X, dtype=float) n, d = X.shape self.mean_ = X.mean(axis=0) Xc = X - self.mean_ if self.standardize: self.scale_ = Xc.std(axis=0, ddof=0) self.scale_[self.scale_ == 0] = 1.0 Xc = Xc / self.scale_ else: self.scale_ = None U, sigma, Vt = np.linalg.svd(Xc, full_matrices=False) k = self.n_components self.components_ = Vt[:k] # (k, d) self.explained_variance_ = (sigma[:k] ** 2) / n total_var = np.sum(sigma ** 2) / n self.explained_variance_ratio_ = self.explained_variance_ / total_var self._n_samples = n return self def transform(self, X): Xc = np.asarray(X, dtype=float) - self.mean_ if self.scale_ is not None: Xc = Xc / self.scale_ return Xc @ self.components_.T def inverse_transform(self, Z): Xc_hat = np.asarray(Z) @ self.components_ if self.scale_ is not None: Xc_hat = Xc_hat * self.scale_ return Xc_hat + self.mean_ def fit_transform(self, X): return self.fit(X).transform(X)

In [8]:

# Demo on 2D data
pca2 = PCAScratch(n_components=1).fit(X_2d)
Z = pca2.transform(X_2d)
X_rec = pca2.inverse_transform(Z)

fig, axes = plt.subplots(1, 2, figsize=(11, 5))

# Original with PC direction
axes[0].scatter(X_2d[:, 0], X_2d[:, 1], alpha=0.3, s=15)
pc = pca2.components_[0]
m = X_2d.mean(axis=0)
axes[0].annotate('', xy=m + pc * 4, xytext=m - pc * 4,
                 arrowprops=dict(arrowstyle='->', color='crimson', lw=2))
axes[0].set_title(f'PC1 direction | EVR = {pca2.explained_variance_ratio_[0]:.1%}')

# Reconstruction
axes[1].scatter(X_rec[:, 0], X_rec[:, 1], alpha=0.4, s=15, color='crimson', label='recon')
axes[1].scatter(X_2d[:, 0], X_2d[:, 1], alpha=0.15, s=10, color='steelblue', label='orig')
axes[1].set_title('1-component reconstruction')
axes[1].legend(fontsize=8)

for ax in axes:
    ax.set_aspect('equal')
    ax.set_xlabel('x₁')
    ax.set_ylabel('x₂')
plt.tight_layout()
plt.show()

# Demo on 2D data pca2 = PCAScratch(n_components=1).fit(X_2d) Z = pca2.transform(X_2d) X_rec = pca2.inverse_transform(Z) fig, axes = plt.subplots(1, 2, figsize=(11, 5)) # Original with PC direction axes[0].scatter(X_2d[:, 0], X_2d[:, 1], alpha=0.3, s=15) pc = pca2.components_[0] m = X_2d.mean(axis=0) axes[0].annotate('', xy=m + pc * 4, xytext=m - pc * 4, arrowprops=dict(arrowstyle='->', color='crimson', lw=2)) axes[0].set_title(f'PC1 direction | EVR = {pca2.explained_variance_ratio_[0]:.1%}') # Reconstruction axes[1].scatter(X_rec[:, 0], X_rec[:, 1], alpha=0.4, s=15, color='crimson', label='recon') axes[1].scatter(X_2d[:, 0], X_2d[:, 1], alpha=0.15, s=10, color='steelblue', label='orig') axes[1].set_title('1-component reconstruction') axes[1].legend(fontsize=8) for ax in axes: ax.set_aspect('equal') ax.set_xlabel('x₁') ax.set_ylabel('x₂') plt.tight_layout() plt.show()

5. Library Comparison¶

In [9]:

try:
    from sklearn.decomposition import PCA as SkPCA
except ModuleNotFoundError:
    print('scikit-learn not installed; skipping PCA comparison.')
else:
    for k in [1, 2]:
        ours = PCAScratch(n_components=k).fit(X_3d)
        skl = SkPCA(n_components=k).fit(X_3d)

        # Sign-invariant alignment per component
        for j in range(k):
            align = abs(float(ours.components_[j] @ skl.components_[j]))
            print(f'k={k}, PC{j+1}: alignment = {align:.10f}')
            assert align > 1 - 1e-8, f'Misaligned PC{j+1}'

        # Explained variance ratio
        evr_diff = np.max(np.abs(ours.explained_variance_ratio_ - skl.explained_variance_ratio_))
        print(f'k={k}, max |EVR diff| = {evr_diff:.2e}')
        assert evr_diff < 1e-10

    # Cross-check SVD route vs eigh route
    ours_svd = PCAScratch(n_components=3).fit(X_3d)
    assert np.allclose(ours_svd.explained_variance_, eigvals_cov, atol=1e-10)
    print('\nSVD route and eigh route give identical eigenvalues. ✓')

try: from sklearn.decomposition import PCA as SkPCA except ModuleNotFoundError: print('scikit-learn not installed; skipping PCA comparison.') else: for k in [1, 2]: ours = PCAScratch(n_components=k).fit(X_3d) skl = SkPCA(n_components=k).fit(X_3d) # Sign-invariant alignment per component for j in range(k): align = abs(float(ours.components_[j] @ skl.components_[j])) print(f'k={k}, PC{j+1}: alignment = {align:.10f}') assert align > 1 - 1e-8, f'Misaligned PC{j+1}' # Explained variance ratio evr_diff = np.max(np.abs(ours.explained_variance_ratio_ - skl.explained_variance_ratio_)) print(f'k={k}, max |EVR diff| = {evr_diff:.2e}') assert evr_diff < 1e-10 # Cross-check SVD route vs eigh route ours_svd = PCAScratch(n_components=3).fit(X_3d) assert np.allclose(ours_svd.explained_variance_, eigvals_cov, atol=1e-10) print('\nSVD route and eigh route give identical eigenvalues. ✓')

k=1, PC1: alignment = 1.0000000000
k=1, max |EVR diff| = 2.22e-16
k=2, PC1: alignment = 1.0000000000
k=2, PC2: alignment = 1.0000000000
k=2, max |EVR diff| = 2.22e-16

SVD route and eigh route give identical eigenvalues. ✓

6. Experiments and Failures — VERIFY¶

In [10]:

# Reconstruction error as k increases
d_big = 8
n_big = 500
X_big = rng.normal(size=(n_big, d_big)) * np.arange(1, d_big + 1)
R_big = np.linalg.qr(rng.normal(size=(d_big, d_big)))[0]
X_big = X_big @ R_big.T

errors = []
for k in range(1, d_big + 1):
    pca_k = PCAScratch(n_components=k).fit(X_big)
    X_rec = pca_k.inverse_transform(pca_k.transform(X_big))
    err = np.mean(np.sum((X_big - X_rec) ** 2, axis=1))
    errors.append(err)

fig, ax = plt.subplots(figsize=(6, 4))
ax.bar(range(1, d_big + 1), errors, color='steelblue', alpha=0.7)
ax.set_xlabel('number of components k')
ax.set_ylabel('mean squared reconstruction error')
ax.set_title('Reconstruction error decreases monotonically with k')
plt.show()

# At full rank, reconstruction should be perfect
assert errors[-1] < 1e-20, 'Full reconstruction should be exact'

# Reconstruction error as k increases d_big = 8 n_big = 500 X_big = rng.normal(size=(n_big, d_big)) * np.arange(1, d_big + 1) R_big = np.linalg.qr(rng.normal(size=(d_big, d_big)))[0] X_big = X_big @ R_big.T errors = [] for k in range(1, d_big + 1): pca_k = PCAScratch(n_components=k).fit(X_big) X_rec = pca_k.inverse_transform(pca_k.transform(X_big)) err = np.mean(np.sum((X_big - X_rec) ** 2, axis=1)) errors.append(err) fig, ax = plt.subplots(figsize=(6, 4)) ax.bar(range(1, d_big + 1), errors, color='steelblue', alpha=0.7) ax.set_xlabel('number of components k') ax.set_ylabel('mean squared reconstruction error') ax.set_title('Reconstruction error decreases monotonically with k') plt.show() # At full rank, reconstruction should be perfect assert errors[-1] < 1e-20, 'Full reconstruction should be exact'

In [11]:

# Orthonormality and zero-mean checks
k_check = 3
pca_check = PCAScratch(n_components=k_check).fit(X_big)
W = pca_check.components_  # (k, d)
Z_check = pca_check.transform(X_big)

# W W^T should be I_k
gram = W @ W.T
assert np.allclose(gram, np.eye(k_check), atol=1e-12), 'Components not orthonormal'
print(f'W·Wᵀ ≈ I_{k_check}: ✓')

# Score means should be ≈ 0
score_means = Z_check.mean(axis=0)
assert np.allclose(score_means, 0, atol=1e-12), 'Score means not zero'
print(f'Score means ≈ 0: ✓ (max abs = {np.max(np.abs(score_means)):.2e})')

# Orthonormality and zero-mean checks k_check = 3 pca_check = PCAScratch(n_components=k_check).fit(X_big) W = pca_check.components_ # (k, d) Z_check = pca_check.transform(X_big) # W W^T should be I_k gram = W @ W.T assert np.allclose(gram, np.eye(k_check), atol=1e-12), 'Components not orthonormal' print(f'W·Wᵀ ≈ I_{k_check}: ✓') # Score means should be ≈ 0 score_means = Z_check.mean(axis=0) assert np.allclose(score_means, 0, atol=1e-12), 'Score means not zero' print(f'Score means ≈ 0: ✓ (max abs = {np.max(np.abs(score_means)):.2e})')

W·Wᵀ ≈ I_3: ✓
Score means ≈ 0: ✓ (max abs = 2.29e-16)

In [12]:

# Failure: scale sensitivity
X_scale_test = X_2d.copy()
X_scale_bad = X_2d.copy()
X_scale_bad[:, 1] *= 8.0

pca_orig = PCAScratch(n_components=1).fit(X_scale_test)
pca_bad = PCAScratch(n_components=1).fit(X_scale_bad)
pca_std = PCAScratch(n_components=1, standardize=True).fit(X_scale_bad)

pc_angle = lambda c: np.degrees(np.arctan2(c[1], c[0]))
print('PC1 angle (original scale):', f'{pc_angle(pca_orig.components_[0]):.1f}°')
print('PC1 angle (x₂ scaled ×8) :', f'{pc_angle(pca_bad.components_[0]):.1f}°')
print('PC1 angle (standardized)  :', f'{pc_angle(pca_std.components_[0]):.1f}°')
print('\nScaling x₂ rotates PC1 toward that axis. Standardization corrects this.')

# Failure: scale sensitivity X_scale_test = X_2d.copy() X_scale_bad = X_2d.copy() X_scale_bad[:, 1] *= 8.0 pca_orig = PCAScratch(n_components=1).fit(X_scale_test) pca_bad = PCAScratch(n_components=1).fit(X_scale_bad) pca_std = PCAScratch(n_components=1, standardize=True).fit(X_scale_bad) pc_angle = lambda c: np.degrees(np.arctan2(c[1], c[0])) print('PC1 angle (original scale):', f'{pc_angle(pca_orig.components_[0]):.1f}°') print('PC1 angle (x₂ scaled ×8) :', f'{pc_angle(pca_bad.components_[0]):.1f}°') print('PC1 angle (standardized) :', f'{pc_angle(pca_std.components_[0]):.1f}°') print('\nScaling x₂ rotates PC1 toward that axis. Standardization corrects this.')

PC1 angle (original scale): -143.5°
PC1 angle (x₂ scaled ×8) : -99.0°
PC1 angle (standardized)  : -135.0°

Scaling x₂ rotates PC1 toward that axis. Standardization corrects this.

In [13]:

# Failure: outlier sensitivity
X_clean = X_2d.copy()
X_outlier = np.vstack([X_2d, [[15.0, -10.0], [-12.0, 8.0]]])

pca_clean = PCAScratch(n_components=1).fit(X_clean)
pca_out = PCAScratch(n_components=1).fit(X_outlier)

align = abs(float(pca_clean.components_[0] @ pca_out.components_[0]))
print(f'PC1 alignment (clean vs with outliers): {align:.4f}')
print(f'PC1 angle clean:   {pc_angle(pca_clean.components_[0]):.1f}°')
print(f'PC1 angle outlier: {pc_angle(pca_out.components_[0]):.1f}°')
print('\nJust 2 extreme points can rotate PC1 significantly.')

# Failure: outlier sensitivity X_clean = X_2d.copy() X_outlier = np.vstack([X_2d, [[15.0, -10.0], [-12.0, 8.0]]]) pca_clean = PCAScratch(n_components=1).fit(X_clean) pca_out = PCAScratch(n_components=1).fit(X_outlier) align = abs(float(pca_clean.components_[0] @ pca_out.components_[0])) print(f'PC1 alignment (clean vs with outliers): {align:.4f}') print(f'PC1 angle clean: {pc_angle(pca_clean.components_[0]):.1f}°') print(f'PC1 angle outlier: {pc_angle(pca_out.components_[0]):.1f}°') print('\nJust 2 extreme points can rotate PC1 significantly.')

PC1 alignment (clean vs with outliers): 0.9919
PC1 angle clean:   -143.5°
PC1 angle outlier: -150.8°

Just 2 extreme points can rotate PC1 significantly.

In [14]:

# Sample size stability: angle to population PC vs n
d_pop = 5
true_cov = np.diag([10, 3, 1, 0.5, 0.1])
R_pop = np.linalg.qr(rng.normal(size=(d_pop, d_pop)))[0]
true_cov = R_pop @ true_cov @ R_pop.T

# Population PC1
evals_pop, evecs_pop = np.linalg.eigh(true_cov)
pc1_pop = evecs_pop[:, -1]

sample_sizes = [10, 20, 50, 100, 200, 500, 1000]
n_trials = 50
median_angles = []
for ns in sample_sizes:
    angles_trial = []
    for _ in range(n_trials):
        Xs = rng.multivariate_normal(np.zeros(d_pop), true_cov, size=ns)
        pca_s = PCAScratch(n_components=1).fit(Xs)
        angle_deg = np.degrees(np.arccos(np.clip(abs(float(pca_s.components_[0] @ pc1_pop)), 0, 1)))
        angles_trial.append(angle_deg)
    median_angles.append(np.median(angles_trial))

fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(sample_sizes, median_angles, 'o-', color='steelblue')
ax.set_xscale('log')
ax.set_xlabel('sample size n')
ax.set_ylabel('median angle to population PC1 (degrees)')
ax.set_title('PC1 converges to population direction as n grows')
plt.show()

# Sample size stability: angle to population PC vs n d_pop = 5 true_cov = np.diag([10, 3, 1, 0.5, 0.1]) R_pop = np.linalg.qr(rng.normal(size=(d_pop, d_pop)))[0] true_cov = R_pop @ true_cov @ R_pop.T # Population PC1 evals_pop, evecs_pop = np.linalg.eigh(true_cov) pc1_pop = evecs_pop[:, -1] sample_sizes = [10, 20, 50, 100, 200, 500, 1000] n_trials = 50 median_angles = [] for ns in sample_sizes: angles_trial = [] for _ in range(n_trials): Xs = rng.multivariate_normal(np.zeros(d_pop), true_cov, size=ns) pca_s = PCAScratch(n_components=1).fit(Xs) angle_deg = np.degrees(np.arccos(np.clip(abs(float(pca_s.components_[0] @ pc1_pop)), 0, 1))) angles_trial.append(angle_deg) median_angles.append(np.median(angles_trial)) fig, ax = plt.subplots(figsize=(6, 4)) ax.plot(sample_sizes, median_angles, 'o-', color='steelblue') ax.set_xscale('log') ax.set_xlabel('sample size n') ax.set_ylabel('median angle to population PC1 (degrees)') ax.set_title('PC1 converges to population direction as n grows') plt.show()

7. Connections¶

• Theory: theory.md — derivations, SVD connection, identifiability
• Eigenvalues & SVD: foundations
• Dimensionality Reduction: synthesis
• Autoencoder: topic 17 — nonlinear PCA analogy

Takeaway¶

• PCA = eigenvectors of covariance = right singular vectors of centered data.
• Max variance ⟺ min reconstruction error.
• SVD route is numerically best; eigh and power iteration agree.
• Always consider standardization; watch for outliers and nonlinear structure.