t-SNE vs. UMAP: The Mathematics of High-Dimensional Visualization Explained

In modern data science and machine learning, datasets rarely conform to simple linear hyperplanes. Whether analyzing 20,000-dimensional single-cell RNA sequencing transcriptomes, 1024-dimensional language model embeddings, or 784-pixel image vectors, complex generative processes wrap data points along twisted, non-linear manifolds. While Principal Component Analysis (PCA) efficiently discovers global linear axes of maximum variance, it cannot bend or unroll non-linear geometries without crushing separate clusters on top of each other.

1. Key Concepts & Mathematical Notation Glossary

Before establishing the mathematical derivations, review the foundational symbols used across non-linear manifold learning literature:

2. Why Non-Linear Dimensionality Reduction Is Essential

Linear dimensionality reduction methods like PCA operate under a rigid geometric constraint: they can only rotate and project data onto flat orthogonal hyperplanes. Consider two classical non-linear geometries where PCA catastrophically fails:

• The Swiss Roll: Imagine a two-dimensional sheet of paper rolled up like a jelly roll in three-dimensional space. The true intrinsic geodesic distance between two points along the paper might be very large, but their Euclidean distance through empty 3D space is tiny. PCA projects across the roll, projecting non-adjacent layers directly onto one another and destroying the manifold.
• Concentric Circles: In two dimensions, consider two nested concentric circles of data points. They are completely separated by an empty ring. However, because both circles share the exact same centroid and have symmetric variance across all radial directions, PCA is incapable of separating them in 1D.

Non-linear manifold learning algorithms relax the requirement of global flatness. Instead of trying to preserve large, straight-line Euclidean distances, they focus on preserving local neighborhood topologies: points that are connected along the local surface stay together, regardless of how the surface twists through high-dimensional space.

3. Conceptual Blueprint: How High-D Distances Become a 2D Map

Both t-SNE and UMAP share a three-phase conceptual workflow:

1. Measure High-Dimensional Closeness: For every point, compute affinities to all other points based on distance. Distant points are given near-zero affinity, while close neighbors receive strong weights.
2. Define Low-Dimensional Affinities: Place the points in a 2D canvas and calculate what the pairwise affinities *would be* if measured in that 2D layout.
3. Iteratively Align the Spaces via Gradient Descent: Compare the high-dimensional affinities against the low-dimensional affinities using a cost function. Compute the gradient of the loss with respect to every point's 2D position, nudging points together or apart until the 2D layout mirrors the original relationships.

4. The Mathematics of t-SNE

Formulated by Laurens van der Maaten and Geoffrey Hinton in 2008, t-Distributed Stochastic Neighbor Embedding (t-SNE) revolutionized high-dimensional data visualization by framing neighborhood preservation as a probability distribution matching problem.

4.1 High-Dimensional Probabilities & Adaptive Perplexity

In the original high-dimensional space, the conditional probability that point $x_i$ picks $x_j$ as its neighbor is defined by a normalized Gaussian kernel centered at $x_i$:

$$p_{j|i} = \frac{\exp\left(-\frac{\|x_i - x_j\|^2}{2\sigma_i^2}\right)}{\sum_{k \neq i} \exp\left(-\frac{\|x_i - x_k\|^2}{2\sigma_i^2}\right)}, \qquad p_{i|i} = 0$$

Notice that the kernel variance $\sigma_i^2$ is not a global constant. Because real datasets have varying density — dense clusters with points packed tightly alongside sparse outliers — a single fixed $\sigma$ would over-smooth dense regions and isolate sparse points. t-SNE tunes $\sigma_i$ individually for each observation $x_i$ such that the Shannon entropy of the conditional distribution $P_i$ matches a user-specified perplexity:

$$H(P_i) = -\sum_{j \neq i} p_{j|i} \log_2 p_{j|i}, \qquad \text{Perp}(P_i) = 2^{H(P_i)}$$

Because entropy $H(P_i)$ is a monotonically increasing function of $\sigma_i$, libraries find the exact value of $\sigma_i$ for every single data point using an efficient binary search (root finding). To eliminate outlier sensitivity and make gradients symmetric, t-SNE defines the joint probability matrix $p_{ij}$ as:

$$p_{ij} = \frac{p_{j|i} + p_{i|j}}{2n}, \qquad \sum_{i, j} p_{ij} = 1$$

4.2 The Geometric Crowding Problem & Student's t-Distribution

In original SNE (Hinton & Roweis, 2002), the low-dimensional similarities $q_{ij}$ were also modeled using a standard Gaussian distribution. This caused a catastrophic failure known as the crowding problem.

t-SNE solves this crowding dilemma by introducing an intentional distribution mismatch: while the high-dimensional space uses a Gaussian kernel, the low-dimensional 2D map uses a Student's t-distribution with 1 degree of freedom (standard Cauchy distribution):

$$q_{ij} = \frac{(1 + \|y_i - y_j\|^2)^{-1}}{\sum_{k} \sum_{l \neq k} (1 + \|y_k - y_l\|^2)^{-1}}, \qquad q_{ii} = 0$$

Because the Cauchy distribution exhibits an inverse-square polynomial tail ($1 / r^2$) rather than an exponential tail ($e^{-r^2}$), it allows moderately separated points in high-D to map to much larger physical distances in 2D without incurring an enormous penalty. This heavy tail creates the open space needed for distinct clusters to separate cleanly.

4.3 Objective Function: Kullback-Leibler (KL) Divergence

t-SNE measures the discrepancy between the true high-dimensional probability distribution $P$ and the low-dimensional layout distribution $Q$ using the Kullback-Leibler (KL) divergence:

$$\text{KL}(P \| Q) = \sum_{i \neq j} p_{ij} \log\left(\frac{p_{ij}}{q_{ij}}\right) = \sum_{i \neq j} p_{ij} \log p_{ij} - \sum_{i \neq j} p_{ij} \log q_{ij}$$

Examine the inherent mathematical asymmetry of KL divergence:

• If points are close in high-D ($p_{ij} \gg 0$): If the low-dimensional map places them far apart ($q_{ij} \to 0$), the ratio $\log(p_{ij} / q_{ij})$ explodes to $+\infty$, incurring a massive penalty. t-SNE exerts an intense attractive force pulling true neighbors together.
• If points are distant in high-D ($p_{ij} \approx 0$): The contribution to the sum is $0 \cdot \log(0 / q_{ij}) \to 0$. t-SNE incurs virtually zero penalty for placing distant points close together or far apart. The only repulsive force comes indirectly from the normalization denominator of $q_{ij}$.

This asymmetry explains why t-SNE preserves local neighborhoods with exquisite precision, but inter-cluster distances and cluster sizes in a t-SNE plot are mathematically uninformative.

The analytical gradient of the KL divergence with respect to the coordinates of point $y_i$ evaluates to a clean physical spring-system equation:

$$\frac{\partial \text{KL}}{\partial y_i} = 4 \sum_{j \neq i} (p_{ij} - q_{ij}) (y_i - y_j) (1 + \|y_i - y_j\|^2)^{-1}$$

5. The Mathematics of UMAP

Introduced in 2018 by Leland McInnes, John Healy, and James Melville, UMAP (Uniform Manifold Approximation and Projection) is grounded in Riemannian geometry and algebraic topology. It treats data as samples drawn from a manifold where Riemannian distance locally behaves like Euclidean distance.

5.1 High-Dimensional Fuzzy Simplicial Sets

Instead of normalized probability distributions, UMAP constructs a weighted k-nearest-neighbor graph called a fuzzy simplicial set. For each point $x_i$, let $\{x_{i_1}, \dots, x_{i_k}\}$ be its $k$ nearest Euclidean neighbors. The directed similarity weight from $x_i$ to $x_j$ is given by:

$$w_{ij} = \exp\left( -\frac{\max(0, \, d(x_i, x_j) - \rho_i)}{\sigma_i} \right)$$

Notice two essential topological parameters:

• Local Connectivity $\rho_i$: $\rho_i = \min \{ d(x_i, x_j) \mid j \neq i \}$, the exact distance from $x_i$ to its nearest neighbor. When $d(x_i, x_j) \le \rho_i$, the numerator becomes $0$, so $\exp(0) = 1$. This ensures that every point is connected to at least its closest neighbor with maximum confidence (weight 1), mathematically enforcing manifold connectivity.
• Bandwidth $\sigma_i$: Solved via binary search such that the sum of weights satisfies: $$\sum_{j=1}^k \exp\left( -\frac{\max(0, \, d(x_i, x_j) - \rho_i)}{\sigma_i} \right) = \log_2(k)$$

To combine asymmetric directed weights $w_{ij}$ and $w_{ji}$ into an undirected fuzzy graph, UMAP applies the algebraic fuzzy set union formula (the t-conorm):

$$w_{ij}^{\text{combined}} = w_{ij} + w_{ji} - w_{ij} \cdot w_{ji}$$

5.2 Low-Dimensional Similarities

In the low-dimensional embedding space, UMAP models similarity using a family of smooth curves governed by parameters $a$ and $b$:

$$q_{ij} = \left( 1 + a \|y_i - y_j\|^{2b} \right)^{-1}$$

The coefficients $a$ and $b$ are determined automatically by fitting a non-linear curve to the user hyperparameter min_dist (which controls how closely packed points can be in 2D). When $a = 1$ and $b = 1$, UMAP's curve exactly matches t-SNE's Student's t-distribution.

5.3 Fuzzy Set Cross-Entropy & Explicit Repulsion

Unlike t-SNE's KL divergence, UMAP optimizes the fuzzy set cross-entropy between the high-dimensional fuzzy graph $W$ and the low-dimensional representation $Q$:

$$\mathcal{L}_{\text{UMAP}} = \sum_{i \neq j} \left[ \underbrace{w_{ij} \log\left(\frac{w_{ij}}{q_{ij}}\right)}_{\text{Attractive Force}} + \underbrace{(1 - w_{ij}) \log\left(\frac{1 - w_{ij}}{1 - q_{ij}}\right)}_{\text{Repulsive Force}} \right]$$

Furthermore, UMAP evaluates this loss using negative sampling (sampling non-edges randomly, similar to Word2Vec). This avoids computing all $O(n^2)$ pairwise interactions, reducing per-iteration complexity to $O(n \cdot k)$ and enabling UMAP to scale effortlessly to millions of data points.

6. Deep Side-by-Side Comparison: t-SNE vs. UMAP

7. Python Implementation: From Scratch & Production Libraries

Below is a self-contained, fully vectorized NumPy implementation of simplified t-SNE illustrating the core mechanics: pairwise Euclidean distance calculation, Gaussian affinities, Student's t low-dimensional kernel, and KL divergence gradient updates.

import numpy as np

def compute_pairwise_affinities(X, sigma=1.0):
    """
    Compute high-dimensional symmetrized Gaussian affinities p_ij.
    """
    n = X.shape[0]
    # Vectorized squared Euclidean distances: ||x_i - x_j||^2
    sum_X = np.sum(np.square(X), axis=1)
    D = np.maximum(0.0, np.add(np.add(-2 * np.dot(X, X.T), sum_X).T, sum_X))

    # Gaussian conditional probabilities p_{j|i}
    P = np.exp(-D / (2.0 * sigma ** 2))
    np.fill_diagonal(P, 0.0)  # Self-affinities are zero
    P = P / np.maximum(1e-12, np.sum(P, axis=1, keepdims=True))

    # Symmetrized joint probabilities: (p_{j|i} + p_{i|j}) / (2n)
    P = (P + P.T) / (2.0 * n)
    return np.maximum(P, 1e-12)


def compute_low_dim_affinities(Y):
    """
    Compute low-dimensional Student's t-distribution similarities q_ij.
    """
    n = Y.shape[0]
    sum_Y = np.sum(np.square(Y), axis=1)
    D = np.maximum(0.0, np.add(np.add(-2 * np.dot(Y, Y.T), sum_Y).T, sum_Y))

    # Student's t-distribution kernel: (1 + ||y_i - y_j||^2)^(-1)
    num = 1.0 / (1.0 + D)
    np.fill_diagonal(num, 0.0)
    Q = num / np.sum(num)
    return np.maximum(Q, 1e-12), num


def simplified_tsne(X, n_components=2, sigma=1.5, n_iter=300, learning_rate=100.0, seed=42):
    """
    Gradient descent optimization of KL divergence between P and Q.
    """
    np.random.seed(seed)
    n = X.shape[0]

    # Step 1: Compute high-dimensional affinities (fixed)
    P = compute_pairwise_affinities(X, sigma=sigma)

    # Step 2: Initialize low-dimensional coordinates with small random Gaussian noise
    Y = np.random.normal(0.0, 1e-4, size=(n, n_components))

    for iteration in range(n_iter):
        # Step 3: Compute current low-dimensional affinities
        Q, num = compute_low_dim_affinities(Y)

        # Step 4: Analytical gradient computation: d(KL)/dy_i
        # grad_i = 4 * sum_j (p_ij - q_ij) * (y_i - y_j) * (1 + ||y_i - y_j||^2)^(-1)
        PQ_diff = P - Q
        grad = np.zeros_like(Y)
        for i in range(n):
            diff = Y[i] - Y  # Shape: (n, n_components)
            mult = (PQ_diff[i][:, np.newaxis] * num[i][:, np.newaxis]) * diff
            grad[i] = 4.0 * np.sum(mult, axis=0)

        # Step 5: Gradient descent update
        Y = Y - learning_rate * grad

    return Y


# --- Verification on Synthetic Clustered Data ---
if __name__ == "__main__":
    np.random.seed(42)
    cluster_1 = np.random.normal(loc=0.0, scale=0.5, size=(10, 5))
    cluster_2 = np.random.normal(loc=5.0, scale=0.5, size=(10, 5))
    X_demo = np.vstack([cluster_1, cluster_2])  # 20 samples in 5 dimensions

    Y_res = simplified_tsne(X_demo, n_components=2, sigma=1.5, n_iter=300, learning_rate=100.0)
    
    print("=== SIMPLIFIED t-SNE EXECUTION ===")
    print(f"Projected output shape: {Y_res.shape}")
    print(f"Cluster 1 centroid:     {np.round(np.mean(Y_res[:10], axis=0), 2)}")
    print(f"Cluster 2 centroid:     {np.round(np.mean(Y_res[10:], axis=0), 2)}")
    print("\nFirst 3 points of Cluster 1 in 2D:")
    print(np.round(Y_res[:3], 2))
    print("First 3 points of Cluster 2 in 2D:")
    print(np.round(Y_res[10:13], 2))

Executing this from-scratch simulation demonstrates complete spatial decoupling between the two 5D clusters in 2D space:

=== SIMPLIFIED t-SNE EXECUTION ===
Projected output shape: (20, 2)
Cluster 1 centroid:     [ 4.73 21.34]
Cluster 2 centroid:     [-4.73 -21.34]

First 3 points of Cluster 1 in 2D:
[[ 4.54 21.35]
 [ 4.79 21.13]
 [ 4.88 21.43]]
First 3 points of Cluster 2 in 2D:
[[-4.94 -21.38]
 [-4.83 -21.21]
 [-4.58 -21.47]]

Production Workflow: Scikit-Learn t-SNE and umap-learn

In production workflows, use official optimized implementations featuring Barnes-Hut or FFT-accelerated interpolation (openTSNE / sklearn.manifold.TSNE) and nearest neighbor descent (umap-learn):

from sklearn.manifold import TSNE
from sklearn.datasets import load_digits
import numpy as np

# 1. Load handwritten digits dataset (1797 samples, 64 pixel features)
digits = load_digits()
X, y = digits.data, digits.target

# --- PRODUCTION t-SNE ---
# init='pca' provides deterministic global initialization and faster convergence
tsne = TSNE(
    n_components=2,
    perplexity=30.0,
    init='pca',
    learning_rate='auto',
    random_state=42
)
X_tsne = tsne.fit_transform(X)
print(f"t-SNE projected shape: {X_tsne.shape}")

# --- PRODUCTION UMAP ---
# Note: requires 'pip install umap-learn'
try:
    import umap
    reducer = umap.UMAP(
        n_neighbors=15,
        min_dist=0.1,
        n_components=2,
        metric='euclidean',
        random_state=42
    )
    X_umap = reducer.fit_transform(X)
    print(f"UMAP projected shape:  {X_umap.shape}")
    
    # Out-of-sample projection demo: transform novel points
    X_new = X[:5] + np.random.normal(0, 0.1, size=X[:5].shape)
    X_new_projected = reducer.transform(X_new)
    print(f"UMAP out-of-sample transform shape: {X_new_projected.shape}")
except ImportError:
    print("Install umap-learn using: pip install umap-learn")

8. Hyperparameter Tuning Guide

Both algorithms are highly sensitive to hyperparameter choices. Understanding how these knobs impact the optimization landscape prevents misleading visual artifacts:

• t-SNE perplexity (Default: 30): Controls the effective number of nearest neighbors each point considers. Low perplexity (2–5) makes the algorithm fixate on fine-grained micro-structure, fracturing continuous manifolds into dozens of fake, scattered islands. High perplexity (50–100) forces the algorithm to account for broader context, merging separate clusters into an undifferentiated blob.
• t-SNE early_exaggeration (Default: 12.0): During initial iterations, high-dimensional affinities $p_{ij}$ are multiplied by this factor. This forces clusters to form tight, compact spheres early in optimization, making it easier for them to navigate past each other without getting trapped in local minima.
• UMAP n_neighbors (Default: 15): Controls the balance between local and global structure. Small values (5–10) focus strictly on fine local geometry, whereas larger values (50–100) capture macroscopic connectivity across the entire manifold.
• UMAP min_dist (Default: 0.1): Controls the minimum Euclidean distance between points in the low-dimensional canvas. Lower values (0.001–0.01) pack points into dense, clumped topological clusters. Higher values (0.5–0.8) prevent clumping and preserve broad spatial dispersion.

9. Common Pitfalls & How to Avoid Misinterpreting Plots

• Never Interpret Cluster Sizes as Real Variance: Because t-SNE expands dense clusters and compresses sparse clusters to achieve constant perplexity entropy, cluster radius in a t-SNE plot does NOT reflect original high-dimensional variance.
• Never Interpret Distances Between Clusters as Semantic Distance: In t-SNE, the KL divergence penalty for distant points is zero. The physical distance between two separate clusters in a t-SNE plot is largely arbitrary and depends heavily on the random seed and iteration count.
• Skipping Pre-Dimensionality Reduction on Massive Feature Spaces: Running t-SNE or UMAP directly on 10,000 raw features often fails due to distance concentration and high computational overhead. Best practice is to run PCA first down to 30–50 dimensions to eliminate noise, and then apply t-SNE or UMAP to the PCA coordinates.
• Using t-SNE / UMAP as Supervised Feature Extractors Blindly: These algorithms are non-parametric visualizers. Standard t-SNE cannot project new test instances, and distorted distances can degrade linear classifiers trained on the 2D coordinates. When feature compression is required for tabular ML pipelines, prefer PCA, UMAP with careful cross-validation, or Autoencoders.

10. Hands-On Practice & Curriculum Roadmap

Consolidate your mastery of non-linear dimensionality reduction with these practical exercises:

1. The Perplexity Sweep: Load the digits dataset and run Scikit-Learn's TSNE with perplexity values of 2, 10, 30, and 100. Observe how cluster cohesion changes from fragmented confetti to a single merged sphere.
2. PCA vs. t-SNE vs. UMAP on the Swiss Roll: Generate a 3D Swiss roll using sklearn.datasets.make_swiss_roll. Plot the 2D projections of all three algorithms side-by-side. Verify that PCA collapses the roll, while UMAP and t-SNE successfully unroll the sheet.
3. Out-of-Sample Validation: Fit a UMAP reducer on 80% of a dataset, project the remaining 20% using reducer.transform(X_test), and verify whether test samples map cleanly into their respective class clusters.