PCA: Geometric and Math Intuition

Principal Component Analysis finds a new coordinate system for your data — one aligned with the directions of maximum spread. The axes (principal components) are rotated to capture as much variance as possible, in decreasing order.

This post builds the full geometric and mathematical picture from scratch — centering, covariance, variance maximization, projection, and reconstruction — all by hand on a 4-sample dataset. Post 04 covers eigendecomposition; post 05 covers the sklearn implementation.

Anchor: 4-sample dataset — hours studied vs exam score. Strong correlation means the data almost lies on a line, so PCA should compress it dramatically.

import numpy as np

X = np.array([
    [1.0, 2.0],
    [2.0, 4.0],
    [3.0, 5.0],
    [4.0, 7.0],
])
feature_names = ['hours_studied', 'exam_score']
print("X shape:", X.shape)
print("Correlation:", np.corrcoef(X[:, 0], X[:, 1])[0, 1].round(4))

Nearly perfect correlation. In 2D, these 4 points almost lie on a straight line. PCA should find that line and express the data in 1 dimension.

What PCA Does (Geometric View)

PCA replaces the original axes (hours_studied, exam_score) with new axes (PC1, PC2) aligned to the data's natural orientation.

• PC1 = the direction with the most variance (longest spread of the projected points)
• PC2 = orthogonal to PC1, captures remaining variance

Step 1: Center the Data

PCA measures variance, which requires zero-mean data. Subtract the column means:

$\bar{X}=\left[\frac{1+2+3+4}{4},\frac{2+4+5+7}{4}\right]=[2.5,4.5]$

X_mean = X.mean(axis=0)
X_c = X - X_mean
print(f"Column means: {X_mean}")
print("\nCentered data X_c:")
for i, (orig, cent) in enumerate(zip(X, X_c)):
    print(f"  Sample {i+1}: {orig} → {cent}")

Step 2: Compute the Covariance Matrix

$C=\frac{1}{n-1}{X}_c^T{X}_c$

Computing each entry manually:

$\sum x_{1c}^2=(-1.5{)}^2+(-0.5{)}^2+(0.5{)}^2+(1.5{)}^2=2.25+0.25+0.25+2.25=5.0$

$\sum x_{2c}^2=(-2.5{)}^2+(-0.5{)}^2+(0.5{)}^2+(2.5{)}^2=6.25+0.25+0.25+6.25=13.0$

$\sum x_{1c}{x}_{2c}=(-1.5)(-2.5)+(-0.5)(-0.5)+(0.5)(0.5)+(1.5)(2.5)=3.75+0.25+0.25+3.75=8.0$

$C=\frac{1}{3}\left[\begin{array}{cc}5.0& 8.0\\ 8.0& 13.0\end{array}\right]=\left[\begin{array}{cc}1.667& 2.667\\ 2.667& 4.333\end{array}\right]$

C = np.cov(X_c, rowvar=False)   # rowvar=False: each column is a variable
print("Covariance matrix C:")
print(C.round(4))

Reading the matrix:

• C[0,0] = 1.667 = Var(x₁_c), std ≈ 1.29
• C[1,1] = 4.333 = Var(x₂_c), std ≈ 2.08
• C[0,1] = 2.667 = Cov(x₁_c, x₂_c) — positive and large → features move together

The correlation coefficient confirms the near-perfect linear relationship:

$r=\frac{C[0,1]}{\sqrt{C[0,0]}\cdot \sqrt{C[1,1]}}=\frac{2.667}{\sqrt{1.667}\cdot \sqrt{4.333}}=\frac{2.667}{1.291\times 2.082}=\frac{2.667}{2.688}=0.992$

Large off-diagonal entries mean features are redundant — they contain overlapping information. PCA will compress this redundancy into fewer dimensions.

The Variance Maximization Objective

PCA finds the unit vector u (the principal component direction) that maximizes the variance of projected scores:

$\text{Var}(\mathbf{u})={\mathbf{u}}^{T}C\mathbf{u}\text{subject to }\Vert \mathbf{u}\Vert =1$

For $\mathbf{u}=[a,b]$ with $a^2+b^2=1$:

${\mathbf{u}}^{T}C\mathbf{u}=a^2\cdot 1.667+2ab\cdot 2.667+b^2\cdot 4.333$

Comparing three candidate directions:

The diagonal direction captures more variance than either original axis. This is the key insight: rotating the coordinate system can reveal structure the original axes hide. Eigendecomposition (post 04) finds the optimal rotation.

Projection onto a Direction

Once we choose direction u, projecting data point ${\mathbf{x}}_c$ onto u gives a scalar score:

$\text{score}={\mathbf{x}}_c\cdot \mathbf{u}=x_{1c}\cdot a+x_{2c}\cdot b$

For u₃ = [1/√2, 1/√2]:

Verifying that the variance of scores equals ${\mathbf{u}}^{T}C\mathbf{u}$:

$\text{Var(scores)}=\frac{(-2.828{)}^2+(-0.707{)}^2+(0.707{)}^2+(2.828{)}^2}{3}=\frac{7.998+0.500+0.500+7.998}{3}=\frac{17.0}{3}=5.667✓$

u3 = np.array([1/np.sqrt(2), 1/np.sqrt(2)])
scores = X_c @ u3
print("Scores:", scores.round(4))
print("Var(scores):", scores.var(ddof=1).round(4))
print("u^T C u:    ", (u3 @ C @ u3).round(4))

The two quantities are identical — projecting onto u and computing score variance is exactly the same as computing ${\mathbf{u}}^{T}C\mathbf{u}$.

Reconstruction from Projection

Keeping only PC1 means approximating each original point from its score:

${\mathbf{x}}_{\text{reconstructed}}=\text{score}\times \mathbf{u}$

For Sample 1 (score = −2.828):

${\mathbf{x}}_{\text{recon}}=-2.828\times [0.707,0.707]=[-2.0,-2.0]$

True centered value: [−1.5, −2.5]. Reconstruction error:

$ϵ=\sqrt{(-1.5-(-2.0){)}^2+(-2.5-(-2.0){)}^2}=\sqrt{0.25+0.25}=0.707$

For Sample 4 (score = +2.828):

${\mathbf{x}}_{\text{recon}}=+2.828\times [0.707,0.707]=[+2.0,+2.0]$

True: [+1.5, +2.5]. Error = $\sqrt{0.25+0.25}=0.707$ (symmetric).

reconstructed = np.outer(scores, u3)
errors = np.sqrt(((X_c - reconstructed)**2).sum(axis=1))
print("Reconstruction errors:", errors.round(4))
print("Mean squared reconstruction error:", (errors**2).mean().round(4))

This reconstruction error is exactly the information lost by keeping only 1 PC. The optimal PC1 (from eigendecomposition) minimizes this reconstruction error — equivalently, maximizes captured variance.

PCA is NOT the Same as Linear Regression

Both PCA and linear regression fit a line through the data, but they minimize different things:

• PCA minimizes perpendicular (orthogonal) distance from each point to the line
• Linear Regression minimizes vertical distance (residuals in y direction)

For this dataset the lines look similar (near-perfect correlation makes them nearly identical). But in general — when the x-axis has much more spread than the y-axis, or when neither axis is clearly "input" vs "output" — the two lines diverge significantly.

The distinction matters for interpretation:

• Use regression when one variable is a prediction target
• Use PCA when all variables are inputs and you want to compress the space without designating one as special

Test Your Understanding

1. The covariance matrix for this dataset has off-diagonal entry C[0,1]=2.667. If you standardize each feature to have unit variance before computing C, what would the off-diagonal entries become? What matrix would you get, and what is it called?

2. We compared three directions and found uᵀCu = 1.667, 4.333, 5.667. The true PC1 captured by eigendecomposition has uᵀCu ≈ 5.94. Why can't uᵀCu exceed 5.94 for any unit vector — and what is the mathematical reason that 5.94 is the maximum?

3. Reconstruction error for sample 1 was 0.707 using u₃=[1/√2, 1/√2]. The true PC1 (eigenvector) would give a smaller reconstruction error. What would the reconstruction error be if you used the true PC1 for sample 1, given that Var(PC1)≈5.94 and Var(PC2)≈0.06? Show the calculation.

4. PCA minimizes perpendicular distance while linear regression minimizes vertical distance. If you rotate the axes 90 degrees (swap x₁ and x₂), the regression line changes slope but the PCA line stays the same. Why? What property of PCA makes it invariant to which variable you call "input" vs "output"?

5. The covariance matrix C is 2×2 for 2D data and 64×64 for the digits dataset. For a dataset with 10,000 genes and 200 patients (n=200, d=10,000), what is the shape of the covariance matrix? Can you invert it? What does this imply for PCA computation when n << d?