What is PCA?

PCA lets you compress an image by keeping only its most important structure and discarding the rest. This post builds that up from eigenvectors and projections, then shows the reconstruction tradeoff in working NumPy code.

Why reduce dimensions at all?

You reduce dimensions because the math gets faster, the data becomes inspectable, and the low-variance directions you drop are usually noise anyway — the orthogonal new axes carry the structure.

The scatter below shows how a single component (the arrow) captures most of the variance in a two-dimensional dataset.

---

Try it on your own image—drag the k slider and watch reconstruction quality degrade as you drop components.

Mathematical Formulation

Let be a set of orthonormal vectors, and let be the matrix whose columns are these vectors, the Principal Component Matrix:

Now let be data points, each in , stacked as rows of the data matrix :

Because is orthonormal, any data point can be written as a linear combination of the principal components:

The coefficient is the projection of onto , just a dot product. This is the core operation PCA exploits: a change of basis followed by truncation. Projecting onto the standard basis just returns the original coordinates — PCA's power comes from rotating this basis to align with the directions of maximum variance.

---

PCA for Image Compression

Images map cleanly onto this framework. Given grayscale images of size , flatten each into a vector of length . Stack them as rows of .

To find the principal components, compute (the unnormalized covariance) and decompose it:

Note: The normalization factor does not affect eigenvectors, so drop it when you only care about basis directions.

Each eigenvector, reshaped to , is an eigenimage: a basis image from which any image in the dataset can be approximately reconstructed.

---

Worked Example: 2×2 Images

To make the algebra concrete, work through a minimal case: 4 images, each pixels.

Flatten each image row-wise and stack to form :

Compute :

The eigenvectors (columns) and eigenvalues of , sorted by descending eigenvalue:

The first column, corresponding to , captures the dominant structure. Reshaped, it becomes:

Reconstruction formula

Any image can be reconstructed from the eigenimages using the outer-product sum:

Using all eigenvectors gives exact reconstruction (up to floating-point error). Using fewer trades fidelity for compression.

---

Python Implementation

The one thing np.linalg.eig won't do for you is sort. Eigenvectors come back in arbitrary order, so without an explicit sort by descending eigenvalue, k=3 can outperform k=4.

import numpy as np

# Data matrix — each row is a flattened image
X = np.array([
    [1, 1, 0, 2],
    [0, 2, 1, 3],
    [1, 0, 0, 2],
    [2, 1, 0, 0],
], dtype=float)

# Un-normalised covariance
XTX = X.T @ X

# Eigendecomposition
eigenvalues, eigenvectors = np.linalg.eig(XTX)

# Sort by descending eigenvalue
order = np.argsort(eigenvalues)[::-1]
eigenvalues  = eigenvalues[order]
eigenvectors = eigenvectors[:, order]  # columns are eigenvectors

eigenvalues  → [23.208  4.916  1.798  0.078]

eigenvectors →
[[ 0.273  0.927  0.148  0.209]
 [ 0.453  0.077 -0.814 -0.356]
 [ 0.154 -0.197 -0.331  0.910]
 [ 0.834 -0.309  0.454 -0.043]]

def reconstruct(image_flat, eigenvectors, k):
    """Reconstruct an image using the top-k eigenvectors."""
    result = np.zeros_like(image_flat, dtype=float)
    for i in range(k):
        ev = eigenvectors[:, i]      # i-th eigenvector (column)
        alpha = ev @ image_flat      # projection coefficient
        result += alpha * ev
    return result

img1_flat = X[0]   # [1, 1, 0, 2]

for k in [4, 2]:
    rec = reconstruct(img1_flat, eigenvectors, k)
    print(f"k={k}:\n{rec.reshape(2, 2)}\n")

k=4 →  [[1.000e+00  1.000e+00]
        [8.882e-16  2.000e+00]]   ✓  exact (floating-point noise only)

k=2 →  [[0.901  1.304]
        [0.208  1.610]]

---

Takeaways

The reconstruction experiment confirms the theoretical picture: PCA finds a basis ordered by variance explained, and truncating that basis is a principled way to trade fidelity for storage. A few things to keep in mind when applying this in practice:

1. Mean-center your data. PCA on uncentred data conflates the mean direction with the first principal component. Subtract the pixel-wise mean across all images before computing .

2. Use SVD for numerical stability. np.linalg.svd on avoids squaring the condition number that np.linalg.eig introduces via .

3. Compression pays off when — storing eigenimages ( values) plus coefficients per image ( values) is cheaper than the raw pixels only in that regime.