Principal Component Analysis

Representation Learning taught us that better descriptions of data can make learning easier. Embeddings gave us dense vectors where similar objects sit near each other.

PCA asks a different representation question: can we compress data while preserving its most important structure?

Imagine tracking Jessie's restaurant visits with hours since last meal, burger price, friends ordering burgers, distance from home, rating, fries quality, noise level, and whether Jessie is trying to diet. Real datasets can have hundreds, thousands, or millions of features.

Many features may be redundant. Hours since last meal, stomach rumbling intensity, and number of times Jessie says she is starving are different measurements, but they all point toward one hidden factor: hunger.

PCA tries to find the hidden directions of variation that explain the data with fewer numbers.

Suppose we plot restaurant visits using two features: hours since last meal and stomach rumbling intensity. The points form a diagonal cloud. When one feature increases, the other usually increases too.

The original representation is two-dimensional, but the data mostly varies along one direction: overall hunger. PCA rotates the coordinate system to find that direction.

If the first component captures almost all of the variation, then one coordinate can replace two features with almost no information loss.

The Algorithm

1. Center The Data

PCA starts by subtracting the mean of each feature. For a data matrix $X$, the centered data is:

Centering matters because PCA cares about variation around the average, not the absolute location of the cloud.

2. Find Directions Of Maximum Variance

Place a line through the centered data and project every point onto that line. PCA asks which line makes the projected points spread out as much as possible. That line is the first principal component.

If projections are spread out, the line preserves information. If projections collapse together, the line loses information.

3. Project Onto Top Components

For a point $x_i\in\mathbb{R}^d$ and first component $v_1$, the projected coordinate is:

If we keep $k$ components, we form:

Here $Z$ is the lower-dimensional representation.

Worked Example

Suppose Jessie has three visits measured by hunger score and stomach rumbling:

The feature means are 4 and 3, so:

The covariance matrix is:

The first eigenvector points along the diagonal:

The second eigenvector points perpendicular to the diagonal:

All variance lies along $v_1$ and none lies along $v_2$. One dimension is enough.

Why Eigenvectors?

The covariance matrix describes how the data varies. An eigenvector of the covariance matrix is a direction that does not get rotated by that transformation. It only gets stretched.

The eigenvalue $\lambda$ tells us how much variance exists along direction $v$. PCA chooses eigenvectors with the largest eigenvalues.

Explained Variance

If PCA returns eigenvalues $80$, $15$, and $5$, the total variance is $100$. The first two components preserve:

This is how PCA helps decide how many dimensions to keep.

Reconstruction

Compression loses information if we drop components. To reconstruct approximately:

If $k=d$, reconstruction is perfect. If $k<d$, reconstruction is approximate. The lost information is the variance in the discarded components.

Scaling Matters

PCA is sensitive to feature scale. If one feature is distance in metres and another is price in dollars, distance may dominate simply because its numbers are larger. We often standardize first:

This implementation follows the exact algorithm: center, compute covariance, decompose, sort components, project, and optionally reconstruct.

import numpy as np

def pca(X, k):
    # 1. center data
    mu = X.mean(axis=0)
    X_centered = X - mu

    # 2. covariance matrix
    cov = (X_centered.T @ X_centered) / (X_centered.shape[0] - 1)

    # 3. eigen decomposition
    eigenvalues, eigenvectors = np.linalg.eigh(cov)

    # 4. sort eigenvalues descending
    idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    # 5. keep top k components
    W = eigenvectors[:, :k]

    # 6. project
    Z = X_centered @ W
    explained_variance_ratio = eigenvalues / eigenvalues.sum()
    return Z, W, mu, eigenvalues, explained_variance_ratio

def reconstruct(Z, W, mu):
    return Z @ W.T + mu