Dimensionality Reduction

In one dimension, a "neighborhood" of size 0.1 around a point covers 10% of a $[0, 1]$ interval. In 100 dimensions, a hypercube of side 0.1 covers $10^{-100}$ of the unit hypercube. Everything is far from everything. This is the curse of dimensionality. Three consequences worth internalizing.

Most of the volume of a high-dim hypercube is in the corners. Pick a uniform random point in the $d$-dim unit hypercube. As $d$ grows, the point is almost certainly far from the center and close to a corner. Visualizing this on a 2D sphere does not generalize.

Distances concentrate. For points drawn from any reasonable distribution in high dimensions, the ratio of max-pairwise-distance to min-pairwise-distance approaches 1. Nearest neighbor and farthest neighbor become almost the same point.

Random projections preserve distances surprisingly well. Johnson-Lindenstrauss (1984) showed that you can project from $d$ to $k = O(\log n / \epsilon^2)$ dimensions with random Gaussian matrices and preserve all pairwise distances to relative error $\epsilon$. This is a remarkable result and the basis for random-projection-based dim reduction.

The pragmatic move: assume your data lives on a low-dimensional manifold embedded in the high-dimensional space, and find a representation that respects that manifold. PCA assumes the manifold is a linear subspace. t-SNE and UMAP assume it is a smooth nonlinear surface.

Principal Component Analysis finds the orthonormal directions in feature space along which the data has the most variance. Here "variance" means the spread of the point cloud along a direction (not the bias-variance sense): project every point onto a candidate axis, and the variance of those projections measures how stretched-out the data is along it.

Formally, given centered data $\mathbf{X} \in \mathbb{R}^{m \times d}$ (each column mean-zero), the first principal component is the unit vector $\mathbf{w}_1$ that maximizes:

subject to $\|\mathbf{w}_1\|_2 = 1$. By Lagrange multipliers (or by recognizing the Rayleigh quotient, a standard maximize-a-ratio trick whose name you can safely ignore here), the solution is the top eigenvector of $\mathbf{X}^T \mathbf{X}$. (For centered data, $\frac{1}{m}\mathbf{X}^T \mathbf{X}$ is the covariance matrix: its diagonal holds each feature's variance and its off-diagonal entries hold how features co-vary in pairs.) An eigenvector of a matrix is a direction the matrix only stretches, never rotates: $\mathbf{A}\mathbf{v} = \lambda\mathbf{v}$, and the stretch factor $\lambda$ is its eigenvalue.

The second principal component is the unit vector that maximizes the same quantity, subject to being orthogonal to the first. It is the second eigenvector. And so on. The $k$-th principal component is the $k$-th eigenvector of $\mathbf{X}^T \mathbf{X}$, ordered by decreasing eigenvalue.

Equivalently (and this is the production formulation): compute the SVD of $\mathbf{X}$:

The SVD factors any matrix into three pieces: $\mathbf{U}$ ($m \times m$) and $\mathbf{V}$ ($d \times d$) have orthonormal columns (unit-length and mutually at right angles), and $\boldsymbol{\Sigma}$ is diagonal, holding the non-negative singular values in decreasing order. The columns of $\mathbf{V}$ are the principal components. The singular values $\sigma_i = \Sigma_{ii}$ relate to the eigenvalues of $\mathbf{X}^T\mathbf{X}$ by $\lambda_i = \sigma_i^2$. If you instead eigendecompose the population covariance $\mathbf{X}^T\mathbf{X}/m$, its eigenvalues are $\sigma_i^2/m$; the sample-variance convention used by sklearn divides by $m-1$. Projecting onto the first $k$ components means taking each point's dot product with those $k$ orthonormal directions, which keeps only its coordinates along them:

pca = PCA(n_components=5)
Z = pca.fit_transform(X)              # mean-center + SVD + project
print(pca.explained_variance_ratio_)
X_recon = pca.inverse_transform(Z)   # lossy decompress

def fit(self, X):
    X = np.asarray(X, dtype=np.float64)
    self.mean = X.mean(axis=0)
    Xc = X - self.mean
    U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
    self.components = Vt[:self.n_components]
    N = X.shape[0]
    total_var = (s ** 2).sum() / max(1, N - 1)
    self.explained_variance = (s[:self.n_components] ** 2) / max(1, N - 1)
    self.explained_variance_ratio = self.explained_variance / total_var
    return self

def transform(self, X):
    return (X - self.mean) @ self.components.T

def inverse_transform(self, Z):
    return Z @ self.components + self.mean

Note (sign flips): the sign of each principal component is arbitrary; SVD does not pin it down. Different implementations may flip signs of components, which means X_2d from sklearn and from your scratch impl can differ in sign. The reconstructions and explained variances are identical; only the absolute coordinates differ.

Caveat (centering): PCA requires centering. Without it, you compute the principal axes of $\mathbf{X}^T \mathbf{X}$ rather than the centered covariance, which gives directions that pass through the origin rather than through the data centroid. Always subtract the mean.

Two heuristics.

Cumulative explained variance. Sort the eigenvalues decreasing, compute cumulative sums normalized to 1. The number of components you keep is the smallest $k$ such that the cumulative explained variance exceeds your threshold (commonly 0.95).

import numpy as np

pca = PCA().fit(X)
cumsum = np.cumsum(pca.explained_variance_ratio_)
d = np.argmax(cumsum >= 0.95) + 1
print(f"Need {d} components for 95% variance")

Or, in one line, sklearn does it for you:

pca = PCA(n_components=0.95)
X_reduced = pca.fit_transform(X)
print(X_reduced.shape)   # (m, d) where d is the auto-chosen count

The elbow. Plot explained variance per component. There is often a sharp "elbow" where additional components add little. Pick the elbow. This is more subjective than the variance-threshold method but is often the right call when downstream models are robust to the exact dimensionality.

Once you have the components, you can compress data by projecting and decompress by inverse-projecting. The decompression is lossy: you cannot recover the variance you threw away.

pca = PCA(n_components=154)   # 95% variance on MNIST
X_reduced = pca.fit_transform(X_train)
X_recovered = pca.inverse_transform(X_reduced)

For MNIST, $784 \to 154$ dimensions preserves 95% of variance and the recovered images are still recognizable. The compression ratio (5.1x) is modest by modern standards (autoencoders and tokenization both do better), but PCA's compression is linear, which means you can do it analytically in either direction and you have a closed form for the reconstruction error.

import numpy as np
recon_error = np.mean((X_train - X_recovered) ** 2)
# Equals sum of variances thrown away (the trailing eigenvalues).

PCA gives you the optimal linear reconstruction. No other linear projection of the same dimension achieves lower reconstruction error. This is the Eckart-Young theorem and it is one of the cleanest optimality results in ML.

PCA materializes the full $m \times d$ data matrix. For large $m$, this fails.

Incremental PCA processes the data in mini-batches, updating the components after each batch. Same end result, much less memory.

from sklearn.decomposition import IncrementalPCA
import numpy as np

n_batches = 100
inc_pca = IncrementalPCA(n_components=154)
for X_batch in np.array_split(X_train, n_batches):
    inc_pca.partial_fit(X_batch)
X_reduced = inc_pca.transform(X_train)

The math: each partial_fit call does a thin SVD of the new batch combined with the current component estimate. This converges to the full PCA as batches accumulate, with a small accuracy loss compared to running full PCA at once.

IncrementalPCA also supports numpy.memmap for datasets that exceed RAM:

filename = "mnist_train.mmap"
X_mm = np.memmap(filename, dtype="float32", mode="readonly", shape=(60000, 784))
inc_pca.fit(X_mm)

The Johnson-Lindenstrauss lemma says you can project from $d$ to $k = O(\log m / \epsilon^2)$ dimensions with a random matrix and preserve all pairwise distances within $(1 \pm \epsilon)$. The projection matrix has i.i.d. Gaussian entries with appropriate scaling.

This is a deterministic guarantee about a random projection. You do not train anything. You generate a random matrix and multiply.

from sklearn.random_projection import GaussianRandomProjection

# Project to a dimension that JL guarantees preserves distances within 10%
gauss_rp = GaussianRandomProjection(eps=0.1)
X_rp = gauss_rp.fit_transform(X)
print(X_rp.shape[1])   # auto-selected by JL bound

When does this beat PCA? Two cases.

Speed. Random projection is $O(m \cdot d \cdot k)$, much faster than PCA's $O(m \cdot d^2 + d^3)$ for large $d$. For NLP-style sparse features ($d$ in the millions), random projection trains in seconds where PCA takes hours.

Privacy. Random projection destroys the original feature semantics in a way that PCA does not (PCA's components are linear combinations that can sometimes be inverted with side information). Random projection is the basis for several differentially-private dim-reduction methods.

SparseRandomProjection is a variant where the projection matrix is sparse, which is even faster. Same JL guarantee.

PCA finds linear principal axes. If your data lives on a curved manifold, the linear axes do not capture it. Recall that the kernel trick lets you operate in implicit nonlinear feature spaces. Kernel PCA applies the trick to PCA.

Instead of computing eigenvectors of $\mathbf{X}^T \mathbf{X}$, compute eigenvectors of the kernel Gram matrix $\mathbf{K}_{ij} = K(\mathbf{x}_i, \mathbf{x}_j)$. The components are eigenvectors in the implicit feature space, and you can project new points via $K(\mathbf{x}, \mathbf{x}_i)$ for training points $\mathbf{x}_i$.

from sklearn.decomposition import KernelPCA

rbf_pca = KernelPCA(n_components=2, kernel="rbf", gamma=0.04)
X_reduced = rbf_pca.fit_transform(X)

For the Swiss Roll dataset (a 2D manifold curled up in 3D), Kernel PCA with the RBF kernel can "unroll" the manifold, while linear PCA only finds the principal axes of the bounding box.

Hyperparameters: kernel (linear / rbf / poly / sigmoid), gamma (for RBF), degree (for poly). Tune by cross-validation on the downstream task.

PCA and kernel PCA assume globally smooth structure. Locally Linear Embedding assumes only local smoothness: each point can be reconstructed as a linear combination of its $k$ nearest neighbors. The embedding preserves these local reconstruction weights.

The algorithm:

1. For each point $x_i$, find the $k$ nearest neighbors and solve for weights $W_{ij}$ such that $x_i \approx \sum_j W_{ij} x_j$ with $\sum_j W_{ij} = 1$.
2. Find a low-dim embedding $\mathbf{z}_i$ that preserves these weights: $\mathbf{z}_i \approx \sum_j W_{ij} \mathbf{z}_j$.

Step 2 reduces to an eigenvalue problem on the matrix $(\mathbf{I} - \mathbf{W})^T (\mathbf{I} - \mathbf{W})$, with the embedding being the eigenvectors corresponding to the smallest nonzero eigenvalues.

from sklearn.manifold import LocallyLinearEmbedding

lle = LocallyLinearEmbedding(n_components=2, n_neighbors=10, random_state=42)
X_unrolled = lle.fit_transform(X_swiss_roll)

LLE works well on smooth manifolds without holes or sharp folds. On real-world data with noise, it is sensitive to the choice of n_neighbors. Too few neighbors and the manifold falls apart; too many and you recover something close to PCA.

t-SNE (t-distributed Stochastic Neighbor Embedding, van der Maaten and Hinton 2008) is the algorithm responsible for almost every "look at how the clusters separate" plot you have seen in an ML paper.

The recipe at a high level:

1. In high-dim space, compute pairwise affinities $p_{ij}$ that represent the probability that point $j$ would be picked as $i$'s neighbor (Gaussian-weighted distance, normalized per point).
2. In low-dim space, define similar affinities $q_{ij}$ using a Student-t distribution (heavier tails than Gaussian).
3. Find the low-dim embedding that minimizes KL divergence $\sum_{ij} p_{ij} \log (p_{ij} / q_{ij})$ via gradient descent. KL divergence measures how far one probability distribution sits from another (here, the low-dim affinities $q$ from the high-dim affinities $p$); it is zero only when they match and grows as they diverge.

The Student-t tails (instead of Gaussian) are critical: they allow distant clusters to spread out without artificially clustering far points together. This is the "t" in t-SNE.

from sklearn.manifold import TSNE
from sklearn.datasets import fetch_openml

mnist = fetch_openml('mnist_784', as_frame=False)
X, y = mnist.data[:5000], mnist.target[:5000]

tsne = TSNE(n_components=2, perplexity=30, random_state=42, init="pca")
X_2d = tsne.fit_transform(X)

Hyperparameters that matter:

• perplexity (5 to 50). Effectively "how many nearest neighbors to consider". Smaller perplexity emphasizes local structure; larger emphasizes global. Try 5, 30, 100 and compare.
• learning_rate (200 to 1000). The Distill paper recommends 200 for most cases.
• init: "pca" is more stable than "random"; use it.

Caveat (the central one): t-SNE distances are not meaningful. The output preserves neighborhoods, not distances. Two clusters that look far apart in the t-SNE plot are not necessarily far apart in the original space. Cluster sizes in t-SNE plots are also meaningless. The Distill explorable "How to Use t-SNE Effectively" (01-explorables/distill-misread-tsne) is mandatory reading before you put a t-SNE plot in a paper.

The explorable shows: t-SNE can produce different cluster arrangements with different perplexities. t-SNE can show clusters that are not really there (especially at low perplexity on small datasets). t-SNE can fail to show clusters that are there (at high perplexity).

Treat t-SNE primarily as a visualization tool. Its stochastic, non-parametric embedding distorts global geometry, so feeding a fitted plot directly to a downstream classifier is usually a poor default and complicates out-of-sample inference. If you evaluate a parametric or otherwise reproducible embedding as a feature transform, fit it inside each training fold and compare it against simpler baselines without leakage.

UMAP (McInnes et al. 2018) is the modern replacement for t-SNE for most visualization tasks. It is faster, scales to millions of points, and (importantly) preserves more global structure.

The math is more involved than t-SNE: UMAP builds a topological representation of the data as a weighted graph, constructs a corresponding low-dim graph, and minimizes cross-entropy between the two. The output looks similar to t-SNE in spirit but with different distortion properties.

# pip install umap-learn
import umap
reducer = umap.UMAP(n_neighbors=15, min_dist=0.1, n_components=2, random_state=42)
X_umap = reducer.fit_transform(X)

Why UMAP often beats t-SNE:

• Speed. UMAP is roughly 10x faster than t-SNE for large datasets.
• Reproducibility. UMAP is deterministic given a seed; t-SNE is not (gradient descent on a non-convex objective).
• Global structure. UMAP preserves more of the macro-arrangement of clusters than t-SNE does. Two clusters that are "actually" far in the high-dim space tend to be far in the UMAP output.
• Embedding new points. UMAP supports umap.UMAP.transform(X_new) for embedding new data; t-SNE does not without re-fitting.

Why UMAP can still mislead: the same caveats apply in milder form. UMAP cluster sizes are not strictly meaningful. UMAP can still produce artifacts at extreme hyperparameter values. Treat both t-SNE and UMAP as visualization tools, not as feature extractors.

For exploratory analysis of a new high-dim dataset, the canonical pipeline in 2026:

1. Standardize the features (StandardScaler).
2. PCA to 50 components. This denoises and speeds up the next step.
3. UMAP to 2 components. For visualization.
4. Cluster the UMAP output (k-means or DBSCAN) to find groups.
5. Sanity check by running a downstream classifier on the PCA output (not the UMAP output, which throws away too much information).

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import umap

pipeline = Pipeline([
    ("scaler", StandardScaler()),
    ("pca", PCA(n_components=50)),
])
X_50 = pipeline.fit_transform(X)
X_2d = umap.UMAP(n_neighbors=15, random_state=42).fit_transform(X_50)

The PCA step before UMAP is the trick. UMAP on raw 4096-dim embeddings is slow and unstable. UMAP on 50-dim PCA output is fast and stable. The PCA preserves enough variance to keep the structure, while throwing away the noise that confuses UMAP.

Classical dim reduction (PCA, kernel PCA, LLE) is mostly obsolete as a feature extractor for modern ML. Neural networks learn their own representations end-to-end, and the learned representations dominate hand-engineered or PCA-based features for every perceptual task.

Where dim reduction still matters:

• Visualization (t-SNE, UMAP for embedding spaces). When you train a CLIP model and want to see what its image embeddings look like, you reach for UMAP.
• Tabular ML pre-processing (PCA). When you have 1000 highly-correlated features in a tabular dataset, PCA to 50 components plus XGBoost beats XGBoost on the raw 1000 features in many cases.
• Anomaly detection. PCA reconstruction error is a baseline anomaly score. Points that PCA cannot reconstruct well are outliers.
• Compression. PCA is the simplest form of dimensionality compression and is still useful when you need a fast, interpretable, linear compression.
• Interpretability. Mechanistic interpretability research uses PCA on activations to find interpretable directions in neural network feature spaces (22-anthropic-recent/2024-scaling-monosemanticity builds on this idea, generalizing to sparse autoencoders).

The single most important modern application is the last one: PCA on neural network internal states reveals interpretable structure, which is the foundation of mech-interp. The same algorithm Pearson invented in 1901 is still load-bearing.