Principal component analysis

What it is and when you reach for it

Principal component analysis takes high-dimensional data and finds a small set of orthogonal directions, the principal components, that capture as much of its variance as possible, so that each point can be described by a handful of coordinates instead of hundreds. It is the workhorse of unsupervised linear dimensionality reduction: no labels, no nonlinearity, just a rotation of the axes onto the directions that matter followed by dropping the ones that do not. Among its neighbors, PCA is the linear special case of the autoencoder (a linear autoencoder with squared loss learns exactly the PCA subspace), the unsupervised cousin of linear discriminant analysis (which maximizes between-class rather than total variance), and the same mathematics as the truncated SVD that powers latent semantic analysis and classical recommender systems. You reach for it to visualize data in two or three dimensions, to compress features before a downstream model, to whiten inputs so every direction has unit variance, to denoise by discarding low-variance components, and as the first thing to try whenever "this has too many correlated columns" is the problem. When the structure you care about is linear, PCA is exact, fast, and has no hyperparameters to tune beyond how many components to keep.

The math

Setup and centering

Let X be the n × d data matrix, one sample per row. PCA always begins by centering: subtract the column mean μ so the data is mean-zero, X_c = X − 1μ^T. This is not optional bookkeeping. PCA is about variance, variance is measured around the mean, and skipping the centering makes the first "principal component" point roughly at the mean of the cloud rather than along its spread. With centered data the sample covariance is C = X_c^TX_c / (n − 1), a d × d symmetric positive-semidefinite matrix whose diagonal holds the per-feature variances and whose off-diagonals hold the covariances.

The variance-maximization view

Ask for the unit direction u along which the projected data has the largest variance. The projection of the centered data onto u is X_cu, and its variance is (X_cu)^T(X_cu) / (n − 1) = u^TCu. So the first component solves

maximize u^TCu subject to ‖u‖ = 1.

A Lagrange multiplier turns the constraint into u^TCu − λ(u^Tu − 1); differentiating and setting to zero gives Cu = λu. The maximizer is an eigenvector of the covariance, and since u^TCu = λ at a solution, the variance captured equals the eigenvalue. The largest eigenvalue's eigenvector is the first principal component, the next largest (orthogonal to it, which symmetric eigenvectors automatically are) is the second, and so on. PCA is exactly the eigendecomposition of the covariance matrix, read as "directions ranked by variance".

The reconstruction view, and why it is the same problem

Now ask a seemingly different question: among all k-dimensional subspaces, which one lets us reconstruct the data with the smallest squared error? Represent the subspace by an orthonormal basis W (d × k, W^TW = I). Projecting a centered point x onto the subspace and back gives the reconstruction WW^Tx, and we minimize

Σ_i ‖x_i − WW^Tx_i‖².

Expand one term using orthonormality of the columns of W: ‖x − WW^Tx‖² = ‖x‖² − ‖W^Tx‖². The first piece is fixed by the data, so minimizing reconstruction error is identical to maximizing Σ_i ‖W^Tx_i‖² = trace(W^TX_c^TX_cW), the projected variance. Minimizing what you throw away is the same as maximizing what you keep, because total variance is conserved. The two views are one problem, and their common solution, the Eckart-Young theorem, is that the optimal W is the top-k eigenvectors of the covariance, equivalently the top-k right singular vectors of X_c.

PCA via the SVD, and why it beats the covariance eigendecomposition

The singular value decomposition writes the centered data as X_c = UΣV^T, with U (n × r) and V (d × r) orthonormal and Σ diagonal with singular values σ₁ ≥ σ₂ ≥ … ≥ 0. Substitute into the covariance:

C = X_c^TX_c / (n − 1) = VΣ²V^T / (n − 1).

This is already the eigendecomposition of C: the right singular vectors V are the principal components, and the eigenvalues are λ_j = σ_j² / (n − 1). So you never need to form C at all, and you should not. The reason is the same conditioning argument that governs least squares: forming X_c^TX_c squares the condition number, κ(C) = κ(X_c)², so it squares the relative error and can drown small-but-real principal directions in rounding noise. Working directly on X_c with the SVD pays only κ(X_c). The size argument reinforces it: when d is large but data lives in far fewer dimensions (think 10,000-pixel images that really span a few hundred directions), C is a 10,000 × 10,000 matrix you never want to materialize, while a thin SVD of X_c works in the smaller of n and d. There is a classic small demonstration of the danger, the Lauchli matrix: for X with rows (1,1,1), (ε,0,0), (0,ε,0), (0,0,ε), the true singular values include σ ≈ ε, but if ε is near √(machine epsilon) then X^TX has off-diagonal-vs-diagonal structure of order ε², which rounds away entirely, and the covariance route reports that direction as having exactly zero variance while the SVD route recovers it. Same data, and only the algorithm that avoided squaring keeps the information.

Explained-variance ratios

Because total variance is the sum of the eigenvalues (the trace is conserved under the rotation), component j explains a fraction

r_j = λ_j / Σ_k λ_k = σ_j² / Σ_k σ_k²

of the variance, and the cumulative sum of these ratios tells you how much you keep by retaining the first k components. This is how you choose k in practice: keep enough components to reach, say, 95% cumulative explained variance, or look for the "elbow" in the scree plot of eigenvalues where they flatten into noise. The ratios are the honest report of the compression: 50 components at 99% cumulative variance means you can store 50 numbers per sample and reconstruct within 1% of the original variance.

A worked 2-D example

Take four points: (2, 0), (0, 2), (−2, 0), (0, −2). The mean is (0, 0), so the data is already centered. The covariance (dividing by n − 1 = 3) is

C = X^TX / 3, with X^TX = [[8, 0], [0, 8]], so C = [[8/3, 0], [0, 8/3]].

This is a scaled identity: every direction has equal variance 8/3, both eigenvalues are 8/3, and PCA correctly reports no preferred axis, a diffuse circular cloud. Now tilt the data into a correlated shape, (3, 1), (1, 3), (−3, −1), (−1, −3), still mean-zero. Then

X^TX = [[20, 12], [12, 20]],   C = [[20/3, 4], [4, 20/3]].

For a symmetric matrix [[a, b], [b, a]] the eigenvectors are always (1, 1)/√2 with eigenvalue a + b and (1, −1)/√2 with eigenvalue a − b. Here that is eigenvalue 20/3 + 4 = 32/3 ≈ 10.67 along the (1, 1)/√2 diagonal, and 20/3 − 4 = 8/3 ≈ 2.67 along the (1, −1)/√2 anti-diagonal. The first component points up the main diagonal exactly as the point cloud does, and the explained-variance ratios are (32/3)/(40/3) = 0.8 and (8/3)/(40/3) = 0.2. Keeping the first component alone retains 80% of the variance and collapses each point onto the diagonal line. You can read the whole result off the 2 × 2 covariance by eye, which is what makes this the example to keep in your head.

        x2
         |        * (1,3)
         |      /
         |    /  PC1 (1,1)/√2, variance 32/3  (80%)
         |  /
  -------+------- x1
       / |
     /   |        PC2 (1,-1)/√2, variance 8/3  (20%)
   * (-3,-1) ... points stretch along the diagonal
          

Implementation, twice

The class below mirrors the scikit-learn contract: fit centers the data and takes the SVD, transform projects new data onto the stored components, inverse_transform maps codes back to the original space, and explained_variance_ratio reports the spectrum. Two details make the output match sklearn exactly. The first is the sign convention: the SVD determines each component only up to sign, so a deterministic rule is needed to make runs reproducible; the standard fix (used by sklearn's svd_flip) forces the largest-magnitude entry of each left singular vector to be positive. The second is the divisor n − 1 for an unbiased variance estimate. Both frameworks use a thin, economy-size SVD so the work scales with min(n, d).

import torch

class PCA:
    """PCA via the SVD of the centered data.

    Never forms the covariance matrix: the right singular vectors
    of X_c are already its eigenvectors, and working on X_c pays
    cond(X_c) instead of cond(X_c)^2.
    """

    def __init__(self, n_components):
        self.k = n_components

    def fit(self, X):
        self.mean_ = X.mean(dim=0, keepdim=True)
        Xc = X - self.mean_
        # full_matrices=False -> thin SVD, work scales with min(n, d)
        U, S, Vh = torch.linalg.svd(Xc, full_matrices=False)

        # sign convention: make the largest-|.| entry of each left
        # singular vector positive, so results are deterministic
        idx = U.abs().argmax(dim=0)
        signs = torch.sign(U[idx, torch.arange(U.shape[1])])
        U, Vh = U * signs, Vh * signs.unsqueeze(1)

        n = X.shape[0]
        self.components_ = Vh[: self.k]              # (k, d) rows are PCs
        self.singular_values_ = S[: self.k]
        var = S ** 2 / (n - 1)                       # eigenvalues of cov
        self.explained_variance_ = var[: self.k]
        self.explained_variance_ratio_ = (var / var.sum())[: self.k]
        return self

    def transform(self, X):
        return (X - self.mean_) @ self.components_.T

    def inverse_transform(self, Z):
        return Z @ self.components_ + self.mean_

import jax.numpy as jnp
from functools import partial
import jax

class PCA:
    """PCA via the SVD of the centered data.

    Never forms the covariance matrix: the right singular vectors
    of X_c are already its eigenvectors, and working on X_c pays
    cond(X_c) instead of cond(X_c)^2.
    """

    def __init__(self, n_components):
        self.k = n_components

    def fit(self, X):
        self.mean_ = X.mean(axis=0, keepdims=True)
        Xc = X - self.mean_
        # full_matrices=False -> thin SVD, work scales with min(n, d)
        U, S, Vh = jnp.linalg.svd(Xc, full_matrices=False)

        # sign convention: make the largest-|.| entry of each left
        # singular vector positive, so results are deterministic
        idx = jnp.argmax(jnp.abs(U), axis=0)
        signs = jnp.sign(U[idx, jnp.arange(U.shape[1])])
        U, Vh = U * signs, Vh * signs[:, None]

        n = X.shape[0]
        self.components_ = Vh[: self.k]              # (k, d) rows are PCs
        self.singular_values_ = S[: self.k]
        var = S ** 2 / (n - 1)                       # eigenvalues of cov
        self.explained_variance_ = var[: self.k]
        self.explained_variance_ratio_ = (var / var.sum())[: self.k]
        return self

    def transform(self, X):
        return (X - self.mean_) @ self.components_.T

    def inverse_transform(self, Z):
        return Z @ self.components_ + self.mean_

When you need only the top few components of a huge matrix, a full SVD is overkill. Power iteration finds the leading direction by repeatedly multiplying a random vector by the covariance and renormalizing: each multiply amplifies the component along the dominant eigenvector by λ₁ and the next by λ₂, so the ratio (λ₂/λ₁)^t shrinks and the vector converges to the top principal component, quickly when there is a clear spectral gap and slowly when the top eigenvalues are close. The crucial trick is to never form C: compute X_c^T(X_cv) as two matrix-vector products, which keeps both the cost and the conditioning on X_c's terms. Randomized SVD, which is what the production libraries use for large matrices, is essentially this idea applied to a block of vectors at once.

import torch

def top_component(Xc, iters=100):
    """Leading principal component by power iteration.

    Applies the covariance implicitly as two matvecs, X_c^T (X_c v),
    so C is never formed. Returns (unit direction, its eigenvalue).
    Converges like (lambda_2 / lambda_1)^iters.
    """
    n, d = Xc.shape
    v = torch.randn(d, dtype=Xc.dtype)
    v /= v.norm()
    for _ in range(iters):
        w = Xc.T @ (Xc @ v)     # = (n - 1) * C v, no C materialized
        v = w / w.norm()
    eig = (v @ (Xc.T @ (Xc @ v))) / (n - 1)   # Rayleigh quotient
    return v, eig

import jax
import jax.numpy as jnp

def top_component(Xc, iters=100, key=None):
    """Leading principal component by power iteration.

    Applies the covariance implicitly as two matvecs, X_c^T (X_c v),
    so C is never formed. Returns (unit direction, its eigenvalue).
    Converges like (lambda_2 / lambda_1)^iters.
    """
    n, d = Xc.shape
    key = key if key is not None else jax.random.PRNGKey(0)
    v = jax.random.normal(key, (d,), Xc.dtype)
    v /= jnp.linalg.norm(v)

    def step(v, _):
        w = Xc.T @ (Xc @ v)     # = (n - 1) * C v, no C materialized
        return w / jnp.linalg.norm(w), None

    v, _ = jax.lax.scan(step, v, None, length=iters)
    eig = (v @ (Xc.T @ (Xc @ v))) / (n - 1)   # Rayleigh quotient
    return v, eig

Using it on a real shape of problem

A common shape: a thousand samples in fifty dimensions whose real structure is low-rank, say five latent factors plus noise. Fit PCA, read the explained-variance ratios, and confirm the first five components dominate while the rest are a noise floor. The reconstruction error using k components should equal the tail sum of the discarded eigenvalues, a fact worth asserting in a test.

import torch
torch.manual_seed(0)

n, d, rank = 1000, 50, 5
Z = torch.randn(n, rank)
loadings = torch.randn(rank, d)
X = Z @ loadings + 0.1 * torch.randn(n, d)   # 5 factors + small noise

pca = PCA(n_components=10).fit(X)
print(pca.explained_variance_ratio_)   # first ~5 large, then a floor
print(pca.explained_variance_ratio_[:5].sum())   # close to 1.0

Xr = pca.inverse_transform(pca.transform(X))
mse = ((X - Xr) ** 2).mean()
print(mse)   # small: kept 10 of 50 directions, discarded near-noise

v, eig = top_component(X - X.mean(0, keepdim=True))
# v aligns with the first PCA component up to sign:
print((v @ pca.components_[0]).abs())   # ~1.0

import jax
import jax.numpy as jnp

key = jax.random.PRNGKey(0)
kz, kl, ke = jax.random.split(key, 3)

n, d, rank = 1000, 50, 5
Z = jax.random.normal(kz, (n, rank))
loadings = jax.random.normal(kl, (rank, d))
X = Z @ loadings + 0.1 * jax.random.normal(ke, (n, d))  # 5 factors + noise

pca = PCA(n_components=10).fit(X)
print(pca.explained_variance_ratio_)   # first ~5 large, then a floor
print(pca.explained_variance_ratio_[:5].sum())   # close to 1.0

Xr = pca.inverse_transform(pca.transform(X))
mse = ((X - Xr) ** 2).mean()
print(mse)   # small: kept 10 of 50 directions, discarded near-noise

v, eig = top_component(X - X.mean(0, keepdims=True))
# v aligns with the first PCA component up to sign:
print(jnp.abs(v @ pca.components_[0]))   # ~1.0

What to expect: the first five ratios should sum to nearly all of the variance with the remaining components forming a low, flat noise floor (their exact values are seed- and machine-dependent), and the power-iteration direction should align with the first PCA component to within a sign flip, an inner product whose absolute value is essentially 1. If the ratios do not fall off sharply, the data is not as low-rank as assumed, which is itself the useful finding.

Applications

The most visible use of PCA is visualization: projecting genomes, embeddings, or sensor logs to two or three dimensions to see cluster structure at a glance. Population genetics famously found that the top two principal components of European genotype data reproduce the map of Europe, and PCA plots are a standard first look at single-cell RNA sequencing data. In classic computer vision, "eigenfaces" are the principal components of face images, an early face-recognition method. Whitening, transforming data so every principal direction has unit variance (divide each projected coordinate by its singular value), is a standard preprocessing step that decorrelates inputs and was part of the ZCA-whitening pipeline in image classification. PCA as a compression front-end speeds up and stabilizes downstream models: reducing to tens of dimensions before k-means or a Gaussian mixture removes the curse of dimensionality that wrecks distance-based methods in raw high-dimensional space, and it is a routine denoising step because low-variance components are often mostly noise.

The link to modern representation learning runs through the SVD. Truncated SVD of a term-document matrix is latent semantic analysis, the ancestor of dense word embeddings, and truncated SVD of a user-item rating matrix is the matrix-factorization recommender that won much of the Netflix Prize era. In deep learning, PCA is the diagnostic you run on learned embeddings: reducing a table of transformer activations or contrastive image features to a few components exposes what varies, and PCA of the weight or gradient space underlies optimization analyses and second-order methods. It also appears inside systems as a compression primitive, for example reducing the dimensionality of vectors before storing them in an approximate-nearest-neighbor index. Whenever a pipeline says "reduce dimensionality first", PCA is the default, and it is default precisely because it is exact, cheap, and hyperparameter-light.

Against the real libraries

scikit-learn's PCA is the same centered-SVD computation with the production niceties: it picks a solver by problem shape, using a full LAPACK SVD for small matrices and Halko's randomized SVD (svd_solver="randomized") for large ones, where computing only the top components you asked for is far cheaper than a full decomposition. It also accepts n_components as a variance fraction (pass 0.95 and it keeps enough components to reach 95%), offers whiten=True, and applies the same svd_flip sign convention the code above reproduces, so outputs match. For sparse or truly enormous data, scikit-learn's TruncatedSVD skips the centering step entirely and factorizes the raw matrix, which is what you want for term-document or user-item matrices where centering would destroy sparsity and change the meaning; that is the LSA and recommender setting. IncrementalPCA handles data that does not fit in memory by updating the decomposition in minibatches. The from-scratch version above is genuinely enough when the data fits in memory, you want the fit on-device and differentiable (the JAX version is), or you are folding PCA into a larger GPU pipeline where a round-trip to sklearn would dominate the cost.

The important library caveat is when not to use PCA at all. PCA is linear: it can only rotate and drop axes, so structure that lives on a curved manifold (a spiral, a swiss roll, well-separated clusters that a linear projection smears together) is invisible to it. For visualization where preserving local neighborhood structure matters more than global variance, the right tools are UMAP and t-SNE (in scikit-learn as TSNE), nonlinear methods that make cluster structure pop in 2-D far better than PCA. The standard, and recommended, practice is to combine them: run PCA first to reduce to, say, 50 dimensions and remove noise, then run UMAP or t-SNE on the result, which is both faster and cleaner than running the nonlinear method on raw high-dimensional data. Use PCA when you want an exact, invertible, variance-ranked linear reduction; reach for UMAP or t-SNE when the goal is a picture of nonlinear cluster structure and you accept that their axes have no quantitative meaning and their global distances are not to be trusted.

Verification is a one-line comparison: fit both on the same fixed seed and check that explained-variance ratios agree, np.allclose(pca.explained_variance_ratio_, sklearn_pca.explained_variance_ratio_), and that components match up to sign after applying svd_flip (which the code already does, so a plain allclose on the components should pass). If the components come out flipped, the sign convention is the culprit, not a bug in the math; if the ratios disagree, check the n vs n − 1 divisor and whether one side centered the data.

Traps and misconceptions

"Skip the centering, it barely matters." It matters completely. Uncentered PCA measures variance around the origin, not around the mean, so its first component points toward the data's center of mass rather than along its spread. If you deliberately want the uncentered version, that is truncated SVD and it is a different tool with different meaning; do not reach it by forgetting to subtract the mean.

"Standardize the features first, always." Whether to scale columns to unit variance before PCA is a real modeling choice, not a default. If features share units and their variances are meaningful (pixel intensities, say), scaling throws away signal. If they are in incomparable units (age in years, income in dollars), the unscaled covariance is dominated by whichever feature has the largest numbers, so you should standardize or the first component just finds the highest-variance unit. Decide deliberately.

"The top components are the useful ones." PCA ranks directions by variance, which is not the same as ranking them by usefulness for your task. A high-variance direction can be lighting or a nuisance covariate, and the signal you care about can hide in a low-variance component. PCA is unsupervised; it does not know your label, so do not assume the discarded components are worthless for prediction.

"Eigendecompose the covariance, that's the definition." It is a correct definition and a worse algorithm. Forming X^TX squares the condition number and can annihilate small principal directions in floating point, as the Lauchli example shows. Take the SVD of the centered data instead; it computes the same answer and keeps the precision you started with.

"PCA will reveal my clusters." Only if they are linearly separated along high-variance directions. PCA is a linear projection, so nonlinear cluster structure (concentric rings, manifolds, clusters distinguished by low-variance features) can be smeared into overlap. When the goal is to see clusters and the structure is nonlinear, that is exactly the case for UMAP or t-SNE, ideally after a PCA denoising pass.