PCA: Geometric and Math Intuition

You've got a table of 50 customer-survey questions and you need a single satisfaction score per respondent. Averaging all 50 ignores that half of them are redundant — two questions asking the same thing inflate the influence of one idea. You need the one dimension that captures what all 50 questions are really measuring together.

PCA solves this by finding the direction in the data with the maximum spread — the axis along which your respondents differ the most. That becomes PC1: a single score that captures as much of the 50-question variance as possible. Everything else becomes PC2, PC3, and beyond, each capturing less.

This post builds the math from scratch on a tiny 4-sample dataset where every computation is traceable. 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))

X shape: (4, 2)
Correlation: 0.9939

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.

The Plan — Four Steps to PCA

We'll walk through PCA in four steps. At each step we compute something new, and by the end we'll have the full transform and know how much information we'd lose by compressing:

1. Center the data — subtract column means so variance is measured around zero
2. Compute the covariance matrix — summarize how features vary together
3. Maximize variance — compare candidate directions to see which captures the most spread
4. Project and reconstruct — compress to 1D and measure the loss

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}")

Column means: [2.5 4.5]

Centered data X_c:
  Sample 1: [1. 2.] → [-1.5 -2.5]
  Sample 2: [2. 4.] → [-0.5 -0.5]
  Sample 3: [3. 5.] → [ 0.5  0.5]
  Sample 4: [4. 7.] → [ 1.5  2.5]

✓ Step 1 complete. Columns now have zero mean. The centered values are symmetric around the origin — this is the reference frame for measuring variance.

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))

Covariance matrix C:
[[1.6667 2.6667]
 [2.6667 4.3333]]

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.

✓ Step 2 complete. Covariance matrix is [[1.667, 2.667], [2.667, 4.333]]. The large off-diagonal (2.667 relative to the diagonals) confirms hours and exam score move together — PCA can exploit this redundancy.

Step 3: Maximize Variance — Compare Candidate Directions

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 — which exceeds even this 5.667.

✓ Step 3 complete. Of the three directions tested, the diagonal [1/√2, 1/√2] wins with 5.667 variance. But it's not the maximum possible — that requires eigendecomposition (post 04).

Step 4: Project and Reconstruct

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))

Scores: [-2.8284 -0.7071  0.7071  2.8284]
Var(scores): 5.6667
u^T C u:     5.6667

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

Measure the Loss: Reconstruction

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))

Reconstruction errors: [0.7071 0.1768 0.1768 0.7071]
Mean squared reconstruction error: 0.25

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.

✓ Step 4 complete. Projected scores have variance 5.667 (matching uᵀCu), and reconstruction error is 0.25 squared per point. That loss is the price of compressing 2D to 1D.

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

Trace Table

Related Concepts

Backward: This post assumes you understand covariance (how two features vary together), the concept of variance (spread along a single axis), and the difference between correlation and causation. The 4-point 2D anchor is designed so every covariance value and eigenvector is hand-traceable.

Forward: The eigen-decomposition that computes PCA's principal components is covered in the next post (PCA via Eigen-Decomposition). From there, the PCA implementation post shows how to apply the transform, and the comparison with LDA shows the supervised counterpart that uses class labels instead of variance.

Honest Limitations

1. PCA maximizes variance, not separability. The direction with maximum variance is not necessarily the direction that separates classes. On a classification task, PCA can discard the very features that distinguish classes (e.g., two well-separated clusters with low within-cluster variance but high between-cluster variance along a low-variance direction).
2. Scale sensitivity. PCA on the raw 2D anchor (x₁ in [1,5], x₂ in [1,5]) treats both axes equally. If x₂ were in [0, 500] instead, it would dominate the first principal component regardless of correlation structure. Standardization is essential when features are not on the same scale.
3. Linear only. PCA finds linear directions of maximum variance. For data with non-linear structure (e.g., a circle or spiral), PCA's first principal component may point in a direction that is meaningless — no single linear direction captures the structure.

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?