Dimensionality Reduction (PCA, t-SNE, UMAP)

Background

Curse of dimensionality (Bellman, 1961) shows distances become meaningless in high dimensions
PCA dates to Karl Pearson (1901) and Harold Hotelling (1933); t-SNE to van der Maaten & Hinton (2008)
Critical preprocessing step before many classical ML algorithms

The Concept

Dimensionality reduction compresses high-dimensional data into fewer dimensions while preserving as much structure as possible. PCA finds directions of maximum variance (linear). t-SNE and UMAP find non-linear low-dimensional embeddings that preserve local neighborhood structure.

How it Links

PCA is used as a preprocessing step for almost every classical ML algorithm
Autoencoders and VAEs are the deep learning analogs of PCA
Connects to matrix factorization, spectral methods, and representation learning in transformers

Concrete Example

A dataset with 100 pixel-based features for 10,000 images. PCA reduces it to 20 dimensions that capture 95% of the variance - making kNN, SVM, and logistic regression both faster and more accurate. The 20 principal components are combinations of pixel patterns (edges, gradients, textures).

Prerequisites

Math you need: linear algebra, Clustering, Bias, variance, over/under-fitting

Learning outcomes

Derive PCA from both the variance-maximization formulation and the Eckart-Young best-rank-k Frobenius approximation, and prove they are equivalent via the SVD.
Explain whitening / sphering and probabilistic PCA, and recognize when kernel PCA is the right generalization.
State the t-SNE objective: KL divergence between a Gaussian neighborhood in input space and a heavy-tailed Student-t neighborhood in embedding space, with perplexity as the local-neighborhood knob.
Explain UMAP's fuzzy simplicial complex framing and why it preserves more global structure than t-SNE in practice.
Pick a method given (linear vs non-linear, scale, need for inverse transform, downstream model vs visualization).

TL;DR

PCA [1][2] is the linear projection onto the top- $k$ eigenvectors of the centered data covariance, equivalently the best rank- $k$ approximation of $X$ in Frobenius norm by the Eckart-Young theorem [3]. The connection is the SVD $X = U \Sigma V^\top$ : the right singular vectors are loadings, $X V$ are scores, and the explained variance ratios are $\sigma_i^2 / \sum_j \sigma_j^2$ . Probabilistic PCA [4] recasts PCA as a Gaussian latent-variable model, enabling missing-data handling and Bayesian extensions. Kernel PCA [5] applies PCA in an RKHS implicitly through a kernel. For non-linear visualization, t-SNE [6] minimizes the KL divergence between Gaussian neighborhoods $p_{j \mid i}$ in input space and heavy-tailed Student-t neighborhoods $q_{ij}$ in embedding space; perplexity controls the effective neighborhood size. UMAP [7] is built on a fuzzy simplicial complex view of the manifold and tends to preserve global structure better than t-SNE while running an order of magnitude faster. Cluster sizes and inter-cluster distances on a t-SNE plot are not interpretable [8].

Intuition

Three mechanisms reduce dimension:

Linear projection (PCA, FA, ICA). Find a small set of orthogonal directions that capture the most variance (or the most independent signal). Output: a global linear map you can apply to new points.
Manifold preservation (t-SNE, UMAP, LLE, Isomap). Build a graph that reflects local neighborhoods in the original space, then place points in a low-dimensional space that preserves that graph. Output: a 2-D or 3-D plot for the eye, not for a downstream model.
Probabilistic latent variables (PPCA, factor analysis, VAE). Posit that observations are generated from low-dimensional latents through a parametric model and learn the parameters by maximum likelihood.

PCA in one sentence

PCA is the rotation of the data such that the first axis has the most variance, the second axis has the most variance subject to being orthogonal to the first, and so on. The Eckart-Young theorem says this is also the unique linear projection (in $\ell_2$ ) that loses the least information when you keep only the first $k$ axes [3].

Math

PCA from variance maximization: top eigenvectors of the covariance matrix

Setup

Let $X \in \mathbb{R}^{n \times d}$ have rows $x_i^\top$ that have already been centered ( $\bar{x} = 0$ ). Define the sample covariance

S = \frac{1}{n - 1}\, X^\top X.

We want a unit vector $w \in \mathbb{R}^{d}$ that maximizes the variance of the projected data $X w$ .

Variance of a projection

Since $X$ is centered, $\mathbb{E}[X w] = 0$ , so

\operatorname{Var}(X w) = w^\top \!\left(\frac{1}{n-1} X^\top X\right)\! w = w^\top S w.

The optimization is

\max_{w}\; w^\top S w \quad \text{s.t.}\quad \lVert w \rVert_2 = 1.

Lagrangian and the eigenvalue equation

Form $\mathcal{L}(w, \lambda) = w^\top S w - \lambda(w^\top w - 1)$ . Setting $\nabla_w \mathcal{L} = 0$ gives

S w = \lambda w.

So critical points are eigenvectors of $S$ , with the constraint enforcing $\lVert w \rVert = 1$ . Multiplying on the left by $w^\top$ : $w^\top S w = \lambda$ , so the variance equals the eigenvalue. The maximum is the largest eigenvalue, attained at its eigenvector $w_1$ .

Iterating: the deflation argument

The second principal component $w_2$ maximizes variance subject to being orthogonal to $w_1$ . The same Lagrangian argument shows $w_2$ is the eigenvector of $S$ with the second-largest eigenvalue. Inductively, the top- $k$ principal components are the top- $k$ eigenvectors of $S$ , and the explained variance of component $i$ is $\lambda_i / \sum_j \lambda_j$ .

PCA from Eckart-Young: best rank-k Frobenius approximation

Statement

Among all matrices $A \in \mathbb{R}^{n \times d}$ of rank at most $k$ , the one that minimizes $\lVert X - A \rVert_F$ is given by the truncated SVD: if $X = U \Sigma V^\top$ with singular values $\sigma_1 \ge \sigma_2 \ge \cdots \ge 0$ , then

A^\star = U_{(k)}\, \Sigma_{(k)}\, V_{(k)}^\top, \qquad \lVert X - A^\star \rVert_F^{2} = \sum_{i = k+1}^{r} \sigma_i^{2}.

This is the Eckart-Young theorem (also Eckart-Young-Mirsky, since Mirsky 1960 extended it to all unitarily invariant norms) [3].

Proof sketch via SVD

Any rank- $k$ matrix $A$ has $\operatorname{rank}(A) \le k$ , so its null space has dimension at least $d - k$ . Choose any orthonormal basis $\{w_1, \dots, w_{d-k}\}$ of $\operatorname{null}(A)$ and observe

\lVert X - A \rVert_F^{2} \ge \sum_{j=1}^{d-k} \lVert (X - A) w_j \rVert_2^{2} = \sum_{j=1}^{d-k} \lVert X w_j \rVert_2^{2}.

By a min-max characterization of singular values (the Courant-Fischer-Weyl theorem applied to the symmetric $X^\top X$ ), the right-hand side is at least $\sum_{j = k+1}^{r} \sigma_j^{2}$ , with equality when the $w_j$ are right singular vectors $v_{k+1}, \dots, v_d$ . The truncated SVD attains this bound. ∎

Equivalence with variance maximization

The covariance is $S = \tfrac{1}{n-1} X^\top X$ , so the eigenvectors of $S$ are exactly the right singular vectors of $X$ , and the eigenvalues of $S$ are $\sigma_i^{2} / (n - 1)$ . The two formulations therefore yield the same projection. The scores are $X V$ , equivalently $U \Sigma$ , so they can be computed without forming $X^\top X$ - which matters numerically when $d$ is large.

Probabilistic PCA: PCA as a Gaussian latent-variable model

The model

Probabilistic PCA [4] posits a $k$ -dimensional latent $z \sim \mathcal{N}(0, I_k)$ and an observed $x \in \mathbb{R}^{d}$ generated by

x = W z + \mu + \epsilon, \qquad \epsilon \sim \mathcal{N}(0, \sigma^2 I_d).

Marginalizing $z$ gives $x \sim \mathcal{N}(\mu, W W^\top + \sigma^2 I)$ .

Maximum likelihood

Fitting by maximum likelihood (taking the gradient of the log-likelihood through the marginal Gaussian) yields the closed-form solution

W_{\text{MLE}} = U_k\, (\Lambda_k - \sigma^2 I)^{1/2}\, R, \qquad \sigma^2_{\text{MLE}} = \frac{1}{d - k}\sum_{j = k+1}^{d} \lambda_j,

where $U_k$ has the top- $k$ eigenvectors of the sample covariance $S$ , $\Lambda_k$ holds the corresponding eigenvalues $\lambda_1 \ge \cdots \ge \lambda_k$ , and $R$ is any orthogonal $k \times k$ rotation. As $\sigma \to 0$ , this recovers exact PCA.

What probabilistic PCA buys

Missing data. EM handles missing entries naturally.
Bayesian extensions. Bayesian PCA places priors on the columns of $W$ to determine $k$ automatically [9].
Mixtures of PCAs. A mixture model where each component is a low-dimensional PPCA captures non-linear manifolds.

The SVD bridge: scores, loadings, whitening

For centered $X = U \Sigma V^\top$ with $\Sigma = \operatorname{diag}(\sigma_1, \dots, \sigma_r)$ :

Loadings. The columns of $V$ are the principal directions in input space.
Scores. $X V = U \Sigma$ gives the projection of each row onto the principal axes.
Explained variance. Component $i$ explains fraction $\sigma_i^2 / \sum_j \sigma_j^2$ of the total variance.
Whitening. $X V \Sigma^{-1} = U$ produces uncorrelated components with unit variance (“sphering”). Useful as a preprocessing step before ICA or some downstream linear models.

In practice, never form $X^\top X$ when $d$ is large; compute the thin SVD of $X$ directly. numpy.linalg.svd(X, full_matrices=False) is the right call.

Kernel PCA

Kernel PCA [5] applies PCA in an RKHS reached by a kernel $K$ . The trick: in the dual formulation of PCA, the principal directions live in the span of the (centered) data, so we need only the inner products $\phi(x_i)^\top \phi(x_j) = K(x_i, x_j)$ .

The recipe:

Form the Gram matrix $K_{ij} = K(x_i, x_j)$ .
Center in feature space: $\tilde{K} = K - \mathbf{1}_n K - K \mathbf{1}_n + \mathbf{1}_n K \mathbf{1}_n$ , where $\mathbf{1}_n$ is the $n \times n$ matrix with all entries $1/n$ .
Eigendecompose $\tilde{K} = U \Lambda U^\top$ .
The projection of $\phi(x)$ onto the $k$ -th principal component is $\sum_i U_{ik} \tilde{K}(x, x_i) / \sqrt{\lambda_k}$ .

Kernel PCA can capture non-linear structure that linear PCA cannot, but it costs $O(n^2)$ memory for the Gram matrix and does not scale beyond a few tens of thousands of points without approximations.

t-SNE: probability-preserving non-linear embedding

t-SNE [6] is a visualization algorithm: it places high-dimensional points in 2-D or 3-D so that local neighborhoods are preserved.

High-dimensional similarities

For each point $i$ , define a Gaussian neighborhood with point-specific bandwidth $\sigma_i$ :

p_{j \mid i} = \frac{\exp(-\lVert x_i - x_j \rVert^{2} / 2 \sigma_i^{2})}{\sum_{k \ne i} \exp(-\lVert x_i - x_k \rVert^{2} / 2 \sigma_i^{2})}, \qquad p_{i \mid i} = 0.

Symmetrize: $p_{ij} = (p_{j \mid i} + p_{i \mid j}) / (2 n)$ .

The bandwidth $\sigma_i$ is chosen by binary search to make the perplexity $2^{H(P_i)}$ match a target value (typically 5 to 50), where $H(P_i)$ is the Shannon entropy of $\{p_{j \mid i}\}_j$ in nats. Perplexity is the effective number of neighbors and is the only really important hyperparameter.

Low-dimensional similarities (Student-t)

For points $y_i \in \mathbb{R}^{2}$ in the embedding, use a heavy-tailed Student-t with one degree of freedom:

q_{ij} = \frac{(1 + \lVert y_i - y_j \rVert^{2})^{-1}}{\sum_{k \ne \ell} (1 + \lVert y_k - y_\ell \rVert^{2})^{-1}}.

Loss

t-SNE minimizes the KL divergence

\mathcal{L}(Y) = \operatorname{KL}(P \,\|\, Q) = \sum_{i \ne j} p_{ij} \log \frac{p_{ij}}{q_{ij}},

via gradient descent. Two design choices give t-SNE its character:

KL is asymmetric. Putting $P$ first penalizes mapping nearby points (large $p_{ij}$ ) to far-apart positions (small $q_{ij}$ ) much more than the reverse, so t-SNE preserves local structure aggressively.
Heavy-tailed $Q$ . The Student-t in low dim has fatter tails than the Gaussian in high dim, which provides “repulsion” that counteracts the crowding problem: in high dim there is more volume far from any given point, and a Gaussian in low dim cannot capture that without piling distant points on top of each other.

Caveats

t-SNE plots are tricky to read

Three things every t-SNE consumer must internalize:

Cluster sizes are not interpretable. A small visual cluster does not mean a tight high-dimensional one. The embedding scales the variance per cluster to be roughly equal.
Distances between clusters are not interpretable. Two clusters that look far apart on the plot may not be far apart in the original space.
Perplexity matters. Run t-SNE at several perplexities (5, 30, 50, 100) and at least two random seeds, and only trust patterns that are stable across runs.

When in doubt, read Wattenberg, Viégas, and Johnson 2016 in Distill - the visual examples make every pitfall obvious [8].

UMAP: fuzzy simplicial complexes

UMAP [7] takes a different theoretical route. The core idea: model each point’s local neighborhood as a fuzzy simplicial set on a Riemannian manifold, take the union of those local sets, and then optimize a low-dimensional layout whose own fuzzy simplicial set has minimum cross-entropy with the high-dimensional one.

In practice the UMAP algorithm reduces to:

For each point $i$ , find its $n_{\text{neighbors}}$ nearest neighbors.
Compute a smooth weight $w_{ij}$ for each edge using a local distance scaling that calibrates the effective neighborhood size (analog of t-SNE’s perplexity).
Symmetrize using a probabilistic t-conorm: $w_{ij} = w_{ij} + w_{ji} - w_{ij} w_{ji}$ .
Initialize the low-dimensional embedding (often via Laplacian eigenmaps).
Optimize cross-entropy between the high-dim and low-dim fuzzy graphs using stochastic gradient descent with negative sampling.

UMAP empirically preserves more global structure than t-SNE because the loss includes a “repel” term over sampled non-edges (rather than t-SNE’s implicit repulsion via $Q$ ), and because the Laplacian-eigenmap initialization already encodes global geometry. It is also typically 5x to 10x faster than t-SNE and scales to a million points routinely.

Comparison

Method	Linear?	Preserves global structure?	Preserves local structure?	Deterministic?	Has inverse transform?	Scales to 1M points?
PCA	yes	yes (variance is global)	no (only via top-k summary)	yes	yes (transpose of W)	yes (randomized SVD)
Kernel PCA	no (linear in feature space)	yes (in RKHS)	no	yes	no (RKHS preimage is hard)	no (n^2 Gram matrix)
t-SNE	no	no (cluster distances are not meaningful)	yes (the design goal)	no (random init + non-convex)	no	no (typical limit ~50k)
UMAP	no	better than t-SNE in practice	yes	no (random init)	parametric variants exist	yes (linear in n with n_neighbors fixed)

Explained variance

A standard PCA diagnostic is the cumulative explained-variance curve. Pick the smallest $k$ such that $\sum_{i \le k} \sigma_i^2 / \sum_j \sigma_j^2$ exceeds your tolerance (typically 0.9 to 0.99). The Plot below shows this curve on a synthetic dataset with three informative directions and seven noise directions, where the elbow at $k = 3$ is sharp.

Implementations

The numpy version implements PCA via thin SVD in about ten lines. The sklearn version covers the production stack including kernel PCA and t-SNE. The umap-learn sketch shows the API for the most common alternative to t-SNE.

python

import numpy as np

def pca(X, k):
  """PCA via thin SVD. Returns (scores, loadings, explained_var_ratio)."""
  Xc = X - X.mean(axis=0)              # center
  U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
  loadings = Vt[:k].T                  # d x k
  scores = U[:, :k] * s[:k]            # n x k
  var = (s ** 2) / max(len(X) - 1, 1)
  explained = var[:k] / var.sum()
  return scores, loadings, explained

rng = np.random.default_rng(0)
# 3 informative directions in a 10-d space, plus noise.
Z = rng.normal(size=(200, 3))
W = rng.normal(size=(3, 10))
X = Z @ W + 0.05 * rng.normal(size=(200, 10))

scores, loadings, explained = pca(X, k=3)
print("explained variance ratio:", np.round(explained, 3))
print("scores[:3]:")
print(np.round(scores[:3], 3))

PCA in 10 lines via thin SVD. Always center first; never form X.T @ X for large d.

Run a from-scratch SVD-based PCA on a 3-D synthetic dataset live and inspect the principal axes plus explained variances:

Idle (packages: numpy)

import numpy as np

rng = np.random.default_rng(0)

# Generate data along three known orthogonal directions plus small noise.
A = np.array([[1.0, 2.0, 0.5],
            [2.0, 0.5, -1.0],
            [0.5, -1.0, 2.0]])
Q, _ = np.linalg.qr(A)                 # orthonormal basis (true loadings)
true_var = np.array([4.0, 2.0, 0.5])   # variance per principal direction
n = 300
Z = rng.normal(size=(n, 3)) * np.sqrt(true_var)
X = Z @ Q.T + 0.05 * rng.normal(size=(n, 3))

# PCA via thin SVD
Xc = X - X.mean(axis=0)
U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
loadings = Vt.T
scores = U * s

var = (s ** 2) / (n - 1)
explained = var / var.sum()

print("singular values    :", np.round(s, 3))
print("explained variance :", np.round(var, 3))
print("explained ratio    :", np.round(explained, 3))
print("loadings (columns are principal axes):")
print(np.round(loadings, 3))

# Sanity: |loadings| should match the true Q up to sign and column order
# (PCA is unique only up to sign and ordering of equal-variance directions).
print("|loadings . Q| (per principal axis):")
print(np.round(np.abs(loadings.T @ Q).max(axis=1), 3))

(no output yet)

Real-world uses

Use case	Method	Why
Visualizing 1k-50k high-d points (embeddings, gene expression)	t-SNE or UMAP	Both preserve local structure and produce a 2-D plot. UMAP is the modern default; t-SNE remains common for legacy comparison.
Decorrelating features before a downstream linear model	PCA + whitening	Removes collinearity that would otherwise destabilize the closed-form / IRLS solution.
Compressing image features	PCA (or autoencoders for non-linear cases)	PCA is the optimal linear compression in MSE; cheap to compute and to invert at decode time.
Data exploration on million-point datasets	UMAP (or PCA + UMAP cascade)	UMAP scales linearly in n. t-SNE typically caps near 50k. Run PCA -> 50 dims -> UMAP -> 2 dims for the cleanest results.
Anomaly detection on tabular data	PCA reconstruction error	Project + reconstruct; points with high reconstruction residual lie off the principal subspace. Cheap and surprisingly effective.

Pitfalls

PCA assumes meaningful linear directions in the input space

If your features are on wildly different scales (years, dollars, kilograms), the largest-variance direction is “the one with the biggest unit”, not “the one with the biggest signal”. Always standardize first. PCA on raw features is one of the most common mistakes in industrial ML.

Never feed t-SNE coordinates to a downstream model

t-SNE is a visualization algorithm. The optimization is non-convex, non-deterministic, and explicitly does not preserve global distance. Two points that land close together in the embedding may be arbitrarily far apart in the original space, and the embedding changes with each random seed. Use PCA, UMAP (with care), or a learned encoder for any downstream modeling task.

UMAP is faster but still has hyperparameters that bite

UMAP’s n_neighbors plays the role of t-SNE’s perplexity: smaller values emphasize local structure, larger values pull global structure into focus. min_dist controls how tightly clusters pack. Defaults are reasonable but not universal; sweep both knobs and at least two random seeds before trusting any visual pattern.

PCA is sign- and order-ambiguous for equal-variance directions

If two singular values are equal (or numerically very close), the corresponding principal axes are determined only up to a rotation in their shared eigenspace. Sign flips in any column of $V$ are also free. Reproducing exact loadings across runs requires fixing a convention, often “make the largest-magnitude entry of each column positive”.

Pop quiz

Check yourself

Score: 0 / 5(0 answered)

1. Which statement about PCA via SVD is correct?
2. The Eckart-Young theorem says...
3. Why does t-SNE use a Student-t with 1 df in the embedding space?
4. What does perplexity control in t-SNE?
5. You have a million 50-d embeddings and want a 2-D scatter plot. Which method is most practical?

Interview Q&A

Derive PCA from the SVD.onsite

Setup. Let $X \in \mathbb{R}^{n \times d}$ be centered ( $\bar{x} = 0$ ). Compute the thin SVD $X = U \Sigma V^\top$ with $\Sigma = \operatorname{diag}(\sigma_1, \dots, \sigma_r)$ and $\sigma_1 \ge \cdots \ge \sigma_r \ge 0$ .

Variance angle. The covariance is $S = X^\top X / (n - 1) = V (\Sigma^2 / (n - 1)) V^\top$ , so the eigenvectors of $S$ are the columns of $V$ and the eigenvalues are $\sigma_i^2 / (n - 1)$ . The variance-maximizing direction is $v_1$ , with variance $\sigma_1^2 / (n-1)$ .

Best-approximation angle. Eckart-Young says the best rank- $k$ approximation in Frobenius norm is $X_k = U_k \Sigma_k V_k^\top$ , and the residual is $\sum_{i > k} \sigma_i^2$ .

Equivalence. Both formulations select the same top- $k$ right singular vectors and produce identical scores $X V_k = U_k \Sigma_k$ . The explained variance ratio of component $i$ is $\sigma_i^2 / \sum_j \sigma_j^2$ .

Practical note. Compute the thin SVD of $X$ directly. Forming $X^\top X$ squares the condition number and is $O(n d^2)$ memory.

Sketch the proof of the Eckart-Young theorem.research

Statement. For $X \in \mathbb{R}^{n \times d}$ with SVD $X = U \Sigma V^\top$ , the matrix $X_k = U_k \Sigma_k V_k^\top$ minimizes $\lVert X - A \rVert_F$ over all $A$ with $\operatorname{rank}(A) \le k$ , with $\lVert X - X_k \rVert_F^{2} = \sum_{i > k} \sigma_i^{2}$ .

Lower bound. Let $A$ be any rank- $k$ matrix. Its null space has dimension at least $d - k$ . Pick orthonormal $\{w_1, \dots, w_{d-k}\} \subset \operatorname{null}(A)$ . Then

\lVert X - A \rVert_F^{2} \ge \sum_{j=1}^{d-k} \lVert (X - A) w_j \rVert_2^{2} = \sum_{j=1}^{d-k} \lVert X w_j \rVert_2^{2}.

By a min-max characterization of singular values (the Courant-Fischer-Weyl theorem applied to $X^\top X$ ), the minimum of $\sum_{j} \lVert X w_j \rVert_2^{2}$ over $(d-k)$ -dimensional subspaces of $\mathbb{R}^{d}$ is attained when $\{w_j\} = \{v_{k+1}, \dots, v_{d}\}$ , the bottom right singular vectors, and equals $\sum_{i > k} \sigma_i^{2}$ .

Upper bound. $X_k = U_k \Sigma_k V_k^\top$ has rank exactly $k$ and $\lVert X - X_k \rVert_F^{2} = \sum_{i > k} \sigma_i^{2}$ by direct computation.

So $X_k$ attains the lower bound and the truncated SVD is optimal. ∎ [3][10]

Why are distances in a t-SNE plot not interpretable globally?onsite

t-SNE minimizes $\operatorname{KL}(P \,\|\, Q)$ where:

$P$ is built from Gaussian neighborhoods in the original space, with a per-point bandwidth chosen so that each point’s effective neighborhood has the target perplexity. Different points get different bandwidths.
$Q$ is built from a heavy-tailed Student-t in the low-dimensional space.

Two consequences:

Per-point variance scaling. Because each point’s $\sigma_i$ is set to a fixed perplexity, dense and sparse regions of the input space are mapped to roughly equal-sized clusters in the output. Cluster sizes therefore lose meaning.
Asymmetric KL. Putting $P$ first heavily penalizes putting nearby high-d points far apart in low-d but barely penalizes putting far high-d points anywhere. So local structure is preserved aggressively at the cost of global structure.

The combined effect: clusters look reasonable, but distances between clusters - and even cluster sizes - are artifacts of perplexity, init, and seed. Wattenberg, Viégas, and Johnson 2016 have side-by-side examples that make this concrete [8].

Derive kernel PCA in the dual formulation.research

Setup. A feature map $\phi : \mathcal{X} \to \mathcal{H}$ to an inner-product space. We want to do PCA on $\{\phi(x_i)\}_{i=1}^{n}$ without computing $\phi(x_i)$ explicitly.

Centered features. Let $\bar{\phi} = \tfrac{1}{n}\sum_i \phi(x_i)$ and define $\tilde{\phi}(x_i) = \phi(x_i) - \bar{\phi}$ . The covariance operator is

C = \frac{1}{n}\sum_{i=1}^{n} \tilde{\phi}(x_i) \tilde{\phi}(x_i)^\top.

Eigenvectors live in the span of features. Any eigenvector $v$ of $C$ with $\lambda \ne 0$ satisfies $C v = \lambda v$ , which after expanding shows $v \in \operatorname{span}\{\tilde{\phi}(x_i)\}_{i=1}^{n}$ . Write $v = \sum_i \alpha_i \tilde{\phi}(x_i)$ .

Kernel formulation. Substituting into $C v = \lambda v$ and taking inner products with $\tilde{\phi}(x_j)$ gives

\frac{1}{n}\, \tilde{K}^{2}\, \alpha = \lambda\, \tilde{K}\, \alpha \;\;\Longleftrightarrow\;\; \tilde{K}\, \alpha = n \lambda\, \alpha,

where $\tilde{K}_{ij} = \langle \tilde{\phi}(x_i), \tilde{\phi}(x_j) \rangle$ is the centered Gram matrix, computable from the original $K_{ij} = K(x_i, x_j) = \phi(x_i)^\top \phi(x_j)$ as $\tilde{K} = K - \mathbf{1}_n K - K \mathbf{1}_n + \mathbf{1}_n K \mathbf{1}_n$ .

Projections. The projection of $\phi(x)$ onto principal axis $k$ is

\langle v_k, \tilde{\phi}(x) \rangle = \sum_{i=1}^{n} \alpha_i^{(k)} \tilde{K}(x, x_i),

normalized by $\sqrt{\lambda_k}$ . The whole computation goes through inner products only - we never form $\phi$ [5].

Build it yourself

A weekend exercise to internalize the family:

Implement PCA from scratch in 10 lines via thin SVD (see the numpy CodeCompare block) and reproduce sklearn’s loadings up to sign on make_blobs.
Verify Eckart-Young numerically: compute $\lVert X - A \rVert_F$ for $A$ = (truncated SVD) versus $A$ = (random rank- $k$ matrix matching $X$ on the trace) and confirm the SVD residual is the smaller of the two, equal to the analytical $\sum_{i > k} \sigma_i^{2}$ .
Run t-SNE on MNIST at perplexities $\{5, 30, 50, 100\}$ with two seeds each and observe how the cluster layout shifts; cluster sizes and inter-cluster distances are not stable.
Compare UMAP at n_neighbors $\in \{5, 15, 50, 200\}$ on the same data and observe how local-vs-global structure changes.

Interactive Demo

PCA Projection Demo

PC1 explained: 97.8%

PC2 explained: 2.2%

Rotation angle: 0°Explained variance: 65.4%

PC1 direction

PC2 direction

Rotation direction

What to try: Rotate the projection line and watch the explained variance change - the maximum is the first principal component.

References

Pearson, K. (1901). On lines and planes of closest fit to systems of points in space. Philosophical Magazine, Series 6, 2(11), 559-572. doi:10.1080/14786440109462720
Hotelling, H. (1933). Analysis of a complex of statistical variables into principal components. Journal of Educational Psychology 24(6), 417-441. doi:10.1037/h0071325
Eckart, C. and Young, G. (1936). The approximation of one matrix by another of lower rank. Psychometrika 1(3), 211-218. doi:10.1007/BF02288367
Tipping, M. E. and Bishop, C. M. (1999). Probabilistic principal component analysis. Journal of the Royal Statistical Society, Series B, 61(3), 611-622. doi:10.1111/1467-9868.00196
Schölkopf, B., Smola, A., and Müller, K.-R. (1998). Nonlinear component analysis as a kernel eigenvalue problem. Neural Computation 10(5), 1299-1319. doi:10.1162/089976698300017467
van der Maaten, L. and Hinton, G. (2008). Visualizing data using t-SNE. Journal of Machine Learning Research 9, 2579-2605. https://www.jmlr.org/papers/v9/vandermaaten08a.html
McInnes, L., Healy, J., and Melville, J. (2018). UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. arXiv preprint. arXiv:1802.03426
Wattenberg, M., Viégas, F., and Johnson, I. (2016). How to Use t-SNE Effectively. Distill. doi:10.23915/distill.00002
Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer. ISBN 978-0-387-31073-2
Hastie, T., Tibshirani, R., and Friedman, J. (2009). The Elements of Statistical Learning, 2nd ed.. Springer. https://hastie.su.domains/ElemStatLearn/

Clustering (K-Means, DBSCAN, GMM)Recommender Systems