t-SNE

t-SNE (t-distributed stochastic neighbor embedding) places each item of a high-dimensional data set, such as an image, a patch or a network’s embedding of an image, at a point in a two- or three-dimensional map, so that items that are near neighbours in the original space are near neighbours in the map. Van der Maaten and Hinton (2008) introduced it as a variation of stochastic neighbor embedding (SNE; Hinton & Roweis, 2002) that is easier to optimize and crowds fewer points into the centre of the map. Several articles on this site report t-SNE plots from the literature: they compare layers of a network in Deep Features, projection-head outputs in SimCLR and simulated against real images in Embodied AI.

t-SNE is a visualization method, not a general dimensionality reduction. Unlike principal component analysis, it has no projection that can be applied to new points, its result depends on a non-convex optimization, and its map keeps neighbourhoods while discarding cluster sizes and most distances between clusters. This article covers the method as defined in the 2008 paper, the tree-based acceleration of van der Maaten (2014), scikit-learn 1.7.2’s implementation and how it departs from both, and an example on synthetic images that measures what a map keeps.

Problem

Given $n$ points $x_1, \dots, x_n$ in a high-dimensional space (pixel vectors, descriptors, network features), find points $y_1, \dots, y_n$ in $\mathbb{R}^2$ or $\mathbb{R}^3$ that can be drawn as a scatterplot and that show the data’s structure. Van der Maaten and Hinton argue that linear methods such as PCA and classical multidimensional scaling concentrate on keeping dissimilar points far apart, while for data near a low-dimensional nonlinear manifold it is usually more important to keep very similar points close together, which a linear map typically cannot do. Their aim was a single map that keeps the local structure and also reveals global structure such as clusters at several scales.

Inputs and Outputs

Input: $n$ vectors (or, in scikit-learn, a precomputed distance matrix), a perplexity $u$ and optimization settings.
Output: an $n \times 2$ (or $n \times 3$ ) array of map coordinates and the final value of the cost, the Kullback–Leibler divergence. The coordinates have no units, axes or orientation of their own: rotating or reflecting a map does not change its cost.

Intuition

Each point looks at the others through a Gaussian centred on it, and turns distances into the probabilities with which it would pick each other point as its neighbour. The width of each Gaussian is set separately so that every point has the same effective number of neighbours, the perplexity. In the map, the same is done with a much heavier-tailed kernel, a Student-t distribution with one degree of freedom. The map points are then moved by gradient descent until the two sets of probabilities agree as well as possible.

Two choices give the method its character. Because each Gaussian adapts to the local density, a tight cluster and a loose one receive the same neighbourhood structure, and the map has no way to show which was tighter. Because the map’s kernel has heavy tails, a pair that is moderately far apart in the original space can be shown much farther apart in the map at little cost, which leaves room to place each point’s close neighbours faithfully. This is the remedy for what van der Maaten and Hinton call the crowding problem: with ten intrinsic dimensions, eleven points can be mutually equidistant, and the volume available at moderate distances grows much faster than in a plane, so a two-dimensional map that keeps small distances must push moderately distant points much too far away. In SNE, which uses a Gaussian in the map as well, many small attractive forces from these too-distant points crush the map together in the centre and prevent gaps between natural clusters from forming.

Algorithm

The basic procedure (van der Maaten & Hinton, 2008, Algorithm 1, with the settings of their §3.4):

Compute squared distances $\lVert x_i - x_j \rVert^2$ between all pairs. The paper’s experiments first reduced the data to 30 dimensions with PCA (§4.2); the 2014 experiments used 50.
For each point $i$ , find by binary search the Gaussian width $\sigma_i$ whose conditional distribution $P_i$ has the chosen perplexity, and compute $p_{j|i}$ .
Symmetrize: $p_{ij} = (p_{j|i} + p_{i|j}) / 2n$ .
Draw the initial map from an isotropic Gaussian $\mathcal{N}(0, 10^{-4} I)$ .
For $T$ iterations: compute the Student-t affinities $q_{ij}$ of the current map, compute the gradient of $\mathrm{KL}(P \Vert Q)$ , and update the map with a momentum term.
During the first iterations, multiply all $p_{ij}$ by an early exaggeration factor.

The paper’s settings for all its maps were an exaggeration of 4 for the first 50 iterations, $T = 1000$ , momentum $\alpha(t) = 0.5$ for $t < 250$ and $0.8$ afterwards, and a learning rate starting at 100 and adapted after every iteration with the scheme of Jacobs (1988), which gradually increases the learning rate in directions in which the gradient is stable. Algorithm 1 prints the update as $Y^{(t)} = Y^{(t-1)} + \eta\, \partial C / \partial Y + \alpha(t)\,(Y^{(t-1)} - Y^{(t-2)})$ ; since $C$ is minimized, a descent step subtracts the gradient, and in our reading the plus sign is a slip.

Mathematical Formulation

Input affinities and perplexity

SNE converts distances into conditional probabilities (2008, Eq. 1):

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

The width $\sigma_i$ is chosen so that the perplexity of $P_i$ equals a value $u$ set by the user:

\mathrm{Perp}(P_i) = 2^{H(P_i)}, \qquad H(P_i) = -\sum_j p_{j|i} \log_2 p_{j|i}.

The paper describes the perplexity as a smooth measure of the effective number of neighbours, notes that it increases monotonically with $\sigma_i$ , which is why a binary search works, and says that SNE’s performance is fairly robust to its value, with typical values between 5 and 50. In dense regions the search finds a small $\sigma_i$ and in sparse regions a large one.

The joint probabilities are the symmetrized conditionals, $p_{ij} = (p_{j|i} + p_{i|j})/2n$ . The paper explains why it does not simply use one Gaussian over all pairs: for an outlier $x_i$ , every such $p_{ij}$ would be tiny, so its map position would hardly affect the cost and would be poorly determined. With the symmetrized definition $\sum_j p_{ij} > 1/2n$ for every point, so each point contributes significantly to the cost.

Output affinities, cost and gradient

In the map t-SNE uses a Student-t kernel with one degree of freedom, the Cauchy distribution (2008, Eq. 4):

q_{ij} = \frac{\left(1 + \lVert y_i - y_j \rVert^2\right)^{-1}}{\sum_{k \ne l} \left(1 + \lVert y_k - y_l \rVert^2\right)^{-1}}, \qquad q_{ii} = 0,

and minimizes the Kullback–Leibler divergence

C = \mathrm{KL}(P \,\Vert\, Q) = \sum_i \sum_{j} p_{ij} \log \frac{p_{ij}}{q_{ij}},

whose gradient, derived in the paper’s Appendix A, is (Eq. 5)

\frac{\partial C}{\partial y_i} = 4 \sum_j (p_{ij} - q_{ij})(y_i - y_j)\left(1 + \lVert y_i - y_j \rVert^2\right)^{-1}.

Each term is a spring along $y_i - y_j$ : it attracts when the map shows the pair farther apart than its $p_{ij}$ asks ( $q_{ij} < p_{ij}$ ) and repels when closer. Because the KL divergence is asymmetric, the cost of showing close neighbours far apart (small $q$ for a large $p$ ) is large, while showing distant points close together costs little; the cost therefore focuses on local structure.

The paper gives two reasons for the Cauchy kernel. For large distances $(1 + d^2)^{-1}$ approaches an inverse-square law, so the map’s joint probabilities are almost invariant to changes in the scale of the map for points far apart, and large clusters far apart interact like single points. And comparing the gradients of SNE, UNI-SNE (a variant with a uniform background term by Cook et al., 2007) and t-SNE as functions of the high- and low-dimensional distances (its Figure 1), it finds that only t-SNE strongly repels dissimilar points shown close together, and that this repulsion stays bounded rather than growing with distance as in UNI-SNE. It also notes that a Student-t distribution is an infinite mixture of Gaussians and that evaluating it needs no exponential.

SNE and symmetric SNE. SNE minimizes $\sum_i \mathrm{KL}(P_i \Vert Q_i)$ over the conditional distributions, with the same Gaussian kernel $\exp(-\lVert y_i - y_j \rVert^2)$ for every point in the map, and has the gradient $2\sum_j (p_{j|i} - q_{j|i} + p_{i|j} - q_{i|j})(y_i - y_j)$ . Symmetric SNE replaces the sum by one divergence between joint distributions, which gives the simpler gradient $4 \sum_j (p_{ij} - q_{ij})(y_i - y_j)$ ; the authors observed that it produced maps as good as asymmetric SNE and sometimes a little better. t-SNE is symmetric SNE with the Student-t kernel in the map. The paper describes SNE’s optimization as needing Gaussian noise added to the map in early iterations, a form of simulated annealing whose schedule had to be tuned, and argues that t-SNE’s long-range forces let it find good local optima without it.

Parameters

Perplexity $u$ . Sets how many neighbours each point’s distribution effectively covers. The paper’s experiments used 40 (2008, Table 1) and the 2014 experiments 50. Wattenberg, Viégas and Johnson (2016) show on small synthetic sets that it has a complex effect on the resulting pictures and that it should be smaller than the number of points, and suggest that getting the most from t-SNE may mean looking at plots at several perplexities. scikit-learn’s default is 30, and it raises an error if the perplexity is not less than the number of samples.
Early exaggeration. Multiplying all $p_{ij}$ by a factor makes almost all $q_{ij}$ , which still sum to 1, too small for their $p_{ij}$ , so the optimization models the large $p_{ij}$ by forming tight, widely separated clusters that can move past each other while the global arrangement is found (2008, §3.4). The 2008 paper used 4 for 50 iterations; van der Maaten (2014, §5.2) used 12 for the first 250 iterations, and reported that the trick becomes increasingly important as the data set grows.
Early compression, the paper’s other trick, adds an $\ell_2$ penalty on the map points’ distances from the origin early in the optimization; scikit-learn does not implement it.
Learning rate, momentum and iterations. See the paper’s settings above. Wattenberg and colleagues found that the most important thing is to iterate until the map is stable, and that maps stopped too early can show “pinched” shapes.
Output dimension and degrees of freedom. The paper fixes one degree of freedom for 2-D and 3-D maps and suggests more degrees of freedom when reducing to more than three dimensions, because the heavy tails would then hold a large part of the probability mass (§6.2).
The tree approximation’s $\theta$ , under Complexity.

Complexity

Computing $P$ and $Q$ and the exact gradient takes $O(n^2)$ time and memory per iteration, which the 2008 paper says makes the standard version infeasible for much more than about 10,000 points. Its §5 instead embeds a random subset of landmark points with affinities from random walks on a neighbourhood graph of all the points. On 6,000 MNIST landmarks with affinities from all 60,000 digits ( $k = 20$ neighbours), a 1-nearest-neighbour classifier had a 10-fold cross-validated error of 5.13% on the two-dimensional map, against 5.75% on the original 784-dimensional digits.

Van der Maaten (2014) approximates both halves of the cost:

Sparse input affinities. $p_{j|i}$ is computed only over the $\lfloor 3u \rfloor$ nearest neighbours of each point, found with a vantage-point tree in $O(uN \log N)$ time, and set to zero for all other points, so $P$ has $O(uN)$ nonzero entries.
Barnes-Hut forces. The gradient splits into attractive terms, $\sum_j p_{ij} q_{ij} Z (y_i - y_j)$ , cheap because $P$ is sparse, and repulsive terms, $-\sum_j q_{ij}^2 Z (y_i - y_j)$ with $Z = \sum_{k \ne l} (1 + \lVert y_k - y_l \rVert^2)^{-1}$ . A quadtree (an octree in 3-D) is built on the current map, and the points of a cell that is small and far from $y_i$ are replaced by their centre of mass, weighted by their number. The paper’s Eq. (4) prints the test as $r_{\text{cell}} / \lVert y_i - y_{\text{cell}} \rVert^2 < \theta$ , with $r_{\text{cell}}$ the cell’s diagonal; with $\theta = 0$ every interaction is computed exactly. The gradient then costs $O(N \log N)$ .

A dual-tree variant that summarizes cell–cell interactions had, “somewhat counterintuitively” in the author’s words, a slightly worse speed–accuracy trade-off. In the text’s summary of its MNIST experiment, Barnes-Hut t-SNE with $\theta = 0.5$ gave embeddings of the same quality as standard t-SNE in 1-nearest-neighbour error, and embedded all 70,000 digits in 751 seconds on a laptop (Intel Core i5 4258U), although the same paragraph says that standard t-SNE ( $\theta = 0$ ) was not run on the full set because it would take too long, so that comparison does not rest on a run of exact t-SNE on all 70,000 digits.

Implementation

scikit-learn 1.7.2’s sklearn.manifold.TSNE (read in its installed source, _t_sne.py, _utils.pyx, _barnes_hut_tsne.pyx and neighbors/_quad_tree.pyx) uses these defaults: n_components=2, perplexity=30.0, early_exaggeration=12.0, learning_rate="auto", max_iter=1000, init="pca", method="barnes_hut", angle=0.5. It differs from the 2008 paper in several ways:

Exaggeration. 12 for the first 250 iterations (_EXPLORATION_MAX_ITER = 250) at momentum 0.5, then momentum 0.8, the 2014 paper’s setting rather than the 2008 one. Gains are adapted per coordinate in Jacobs’s manner (+0.2 when the gradient keeps its sign, ×0.8 otherwise, floor 0.01).
Learning rate. "auto" sets it to $\max(N / E / 4, 50)$ , with $E$ the early exaggeration, which is 50 for the few hundred points of the example below. The user guide describes the heuristic of Belkina et al. (2019) as the sample size divided by the exaggeration; the docstring explains the extra factor of 4: scikit-learn’s learning rate is four times the one used by several other implementations, so its 200 corresponds to their 800. The default changed to "auto" in version 1.2, and n_iter was renamed max_iter in 1.5.
Initialization. "pca" (the default since 1.2) takes the first two principal components, computed with the randomized solver, and rescales them so that the first has standard deviation $10^{-4}$ , the standard deviation of scikit-learn’s own init="random"; the paper’s random start, $\mathcal{N}(0, 10^{-4} I)$ , has variance $10^{-4}$ , that is, standard deviation $10^{-2}$ . The user guide says PCA initialization mitigates the lack of explicitly preserved global structure.
Precision. Squared distances are cast to float32 for the binary search, which works in nats against $\ln u$ with tolerance $10^{-5}$ and at most 100 steps on the precision $\beta_i = 1/2\sigma_i^2$ ; the map is optimized in float32. With method="exact", joint probabilities are clipped from below at the float64 machine epsilon, $2.22 \times 10^{-16}$ (_joint_probabilities); the Barnes-Hut path’s sparse $P$ (_joint_probabilities_nn) has no such floor, and pairs outside the neighbour lists are zero.
Degrees of freedom. $\max(d - 1, 1)$ for a $d$ -dimensional map, citing van der Maaten’s parametric t-SNE paper of 2009: 1 for 2-D maps, as in the 2008 paper, but 2 for 3-D maps, where the 2008 paper uses 1.
Barnes-Hut. The neighbour count is $\min(n - 1, \lfloor 3u + 1 \rfloor)$ , one more than the paper’s $\lfloor 3u \rfloor$ . The summary test compares the squared cell width with the squared distance against angle squared, that is width/distance < angle, which the docstring calls the angular size of the cell; the paper’s printed Eq. (4) divides by the squared distance instead. Barnes-Hut is limited to at most three output dimensions; method="exact" computes the $O(N^2)$ gradient.
Stopping. Progress is checked every 50 iterations; the optimization stops early after n_iter_without_progress=300 iterations without a lower cost, or when the gradient norm falls below min_grad_norm=1e-7.

What a map keeps, measured on synthetic images

The script below has three blocks; all its data are synthetic images with known generating parameters.

A, a closed loop. 360 grayscale $21 \times 21$ images of a bright bar (Gaussian profile of standard deviation 1 px across, length 15 px with soft ends) through the image centre, at an angle $\theta$ drawn uniformly from $[0°, 180°)$ , with Gaussian pixel noise of standard deviation 0.05. A bar at $\theta$ and at $\theta + 180°$ is the same image, so the images lie on a closed curve whose true coordinates are $(\cos 2\theta, \sin 2\theta)$ ; the first printed line checks this. Each map is scored by its trustworthiness (below), by the winding number of the points, taken in order of $\theta$ , around the map’s centroid ( $\pm 1$ when the cyclic order survives, either sign since a map may be mirrored), and by the largest step between $\theta$ -neighbours in the map divided by the median step, which grows when the loop is broken. The row “true” scores the true coordinates themselves.
B, known sizes and distances. 240 images of $16 \times 48$ pixels in three clusters of 80: A shows a blob in the left third, B a ring in the middle, C a square on the right. The three templates occupy separate cells and are scaled to unit norm, so they are orthogonal; with contrasts 1, 1 and 7 the distances between cluster centres are $\sqrt{2}$ , $\sqrt{50}$ and $\sqrt{50}$ , a ratio $d_{AC}/d_{AB} = 5$ . Within each cluster the images vary only by a horizontal and a vertical illumination ramp with random amplitudes of standard deviation 0.01, 0.04 and 0.16, so the true cluster sizes (RMS radius about the centroid) are in the ratio 1 : 4 : 16. The row “true” measures the sample in pixel space. For the maps, the columns give the RMS radius ratios B/A and C/A, the centroid distance ratios, and the share of the joint probability $P$ that links cluster A to B and A to C.
C, the mechanics. A NumPy computation of $P$ by binary search on each $\beta_i$ is checked against scikit-learn’s _joint_probabilities, a private function that may change between versions; the gradient of Eq. 5 is checked against central finite differences of the cost; and the cost and gradient are checked against scikit-learn’s exact objective _kl_divergence.

Maps use TSNE with its defaults except the perplexity (5, 30 and 100) and random_state. The loop over seeds 0, 1 and 2 draws new data and sets the map’s seed. Trustworthiness (Venna & Kaski, 2001, as scikit-learn’s docstring defines it) penalizes each of a point’s $k$ nearest neighbours in the map that is not among its $k$ nearest in the original space, in proportion to its rank $r(i, j)$ there:

T(k) = 1 - \frac{2}{nk(2n - 3k - 1)} \sum_{i=1}^{n} \sum_{j \in \mathcal{N}_i^k} \max(0,\, r(i, j) - k),

so $T = 1$ when every map neighbour is a true neighbour. The script uses $k = 10$ . The NumPy code was written from the paper’s equations, not ported from scikit-learn. It ran with NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.7.2 on one CPU thread in about 9 seconds, and two runs gave identical output.

import os
os.environ["OMP_NUM_THREADS"] = os.environ["OPENBLAS_NUM_THREADS"] = os.environ["MKL_NUM_THREADS"] = "1"
import numpy as np
from scipy.spatial.distance import squareform
from sklearn.manifold import TSNE, trustworthiness
from sklearn.manifold._t_sne import _joint_probabilities, _kl_divergence  # private API

SEEDS, PERPS = [0, 1, 2], [5, 30, 100]

def bar_images(theta_deg, rng, noise=0.05, size=21, half_len=7.0, width=1.0):
    """21 x 21 images of a bar through the centre at angle theta (degrees), soft edges and ends."""
    c = (size - 1) / 2
    yy, xx = np.mgrid[0:size, 0:size] - c
    t = np.deg2rad(theta_deg)[:, None, None]
    v = xx * np.cos(t) + yy * np.sin(t)                  # along the bar
    u = -xx * np.sin(t) + yy * np.cos(t)                 # across the bar
    ends = np.exp(-np.maximum(np.abs(v) - half_len, 0) ** 2 / (2 * width ** 2))
    img = np.exp(-u ** 2 / (2 * width ** 2)) * ends
    if noise > 0:
        img = img + rng.normal(0, noise, img.shape)
    return img.reshape(len(theta_deg), -1)

def cluster_images(rng, n_per=80, kappa=(1, 1, 7), spread=(0.01, 0.04, 0.16)):
    """16 x 48 images: cluster k = kappa_k * template_k + a*G1 + b*G2, (a, b) ~ N(0, spread_k^2 I)."""
    H, W = 16, 48
    yy, xx = np.mgrid[0:16, 0:16] - 7.5
    r = np.hypot(yy, xx)
    cells = [np.exp(-r ** 2 / (2 * 3.0 ** 2)),                    # A: blob
             np.exp(-(r - 5.0) ** 2 / (2 * 1.0 ** 2)),            # B: ring
             ((np.abs(yy) < 4) & (np.abs(xx) < 4)).astype(float)] # C: 8 x 8 square
    T = np.zeros((3, H, W))
    for k in range(3):
        T[k, :, 16 * k:16 * (k + 1)] = cells[k]                  # each in its own cell
    T = T.reshape(3, -1)
    T /= np.linalg.norm(T, axis=1, keepdims=True)                # orthonormal templates
    gy, gx = np.mgrid[0:H, 0:W]
    G = np.stack([gx - gx.mean(), gy - gy.mean()]).reshape(2, -1).astype(float)
    G /= np.linalg.norm(G, axis=1, keepdims=True)                # two illumination ramps
    X, lab = [], []
    for k in range(3):
        ab = rng.normal(0, spread[k], (n_per, 2))
        X.append(kappa[k] * T[k] + ab @ G)
        lab += [k] * n_per
    return np.vstack(X), np.array(lab), T, G

def winding_and_gap(Y, theta):
    """Winding number of the theta-ordered closed polygon around the map centroid,
    and the largest step between theta-neighbours divided by the median step."""
    Z = Y[np.argsort(theta)] - Y.mean(0)
    ang = np.arctan2(Z[:, 1], Z[:, 0])
    d = np.diff(np.r_[ang, ang[0]])
    d = (d + np.pi) % (2 * np.pi) - np.pi
    step = np.linalg.norm(np.diff(np.r_[Z, Z[:1]], axis=0), axis=1)
    return d.sum() / (2 * np.pi), step.max() / np.median(step)

def p_numpy(D2, u, tol=1e-10, steps=200):
    """Joint P of t-SNE: binary search on beta_i = 1/(2 sigma_i^2) so that 2^H(P_i) = u."""
    n = len(D2)
    Pc, perp_err = np.zeros((n, n)), 0.0
    for i in range(n):
        d = np.delete(D2[i], i)
        d = d - d.min()                                  # cancels in the normalisation
        beta, lo, hi = 1.0, 0.0, np.inf
        for _ in range(steps):
            w = np.exp(-beta * d)
            p = w / w.sum()
            H = -np.sum(p[p > 0] * np.log(p[p > 0]))     # entropy in nats
            if abs(H - np.log(u)) < tol:
                break
            if H > np.log(u):                            # too many neighbours: narrow
                lo, beta = beta, (beta * 2 if hi == np.inf else (beta + hi) / 2)
            else:                                        # too few: widen
                hi, beta = beta, (beta + lo) / 2
        perp_err = max(perp_err, abs(2 ** (H / np.log(2)) - u) / u)   # 2^(bits) as in the paper
        Pc[i, np.arange(n) != i] = p
    return (Pc + Pc.T) / (2 * n), perp_err

def kl_grad(P, Y):
    """KL(P||Q) with the Student-t Q, and its gradient (paper Eq. 5)."""
    D2 = np.sum((Y[:, None] - Y[None]) ** 2, -1)
    W = 1.0 / (1.0 + D2)
    np.fill_diagonal(W, 0)
    Q = W / W.sum()
    m = P > 0
    C = np.sum(P[m] * np.log(P[m] / Q[m]))
    g = 4 * np.sum(((P - Q) * W)[:, :, None] * (Y[:, None] - Y[None]), axis=1)
    return C, g

# ---------------- A: a closed loop of oriented bars ----------------
z = bar_images(np.array([37.0, 217.0, 0.0, 179.0, 90.0]), None, noise=0)
print(f"A check: max|img(37)-img(217)| = {np.abs(z[0] - z[1]).max():.1e}; "
      f"distance 0-179 deg {np.linalg.norm(z[2] - z[3]):.3f}, 0-90 deg {np.linalg.norm(z[2] - z[4]):.3f}")
print("A: perp seed   T(k=10)  T_ref   winding  gap")
for seed in SEEDS:
    rng = np.random.default_rng(seed)
    th = rng.uniform(0, 180, 360)
    X = bar_images(th, rng)
    ref = np.c_[np.cos(np.deg2rad(2 * th)), np.sin(np.deg2rad(2 * th))]
    t_ref = trustworthiness(X, ref, n_neighbors=10)
    w, gap = winding_and_gap(ref, th)
    print(f"   true {seed:4d}   {t_ref:.4f}  {t_ref:.4f}  {w:+.3f}   {gap:5.2f}")
    for u in PERPS:
        Y = TSNE(perplexity=u, random_state=seed).fit_transform(X)
        w, gap = winding_and_gap(Y, th)
        print(f"   {u:4d} {seed:4d}   {trustworthiness(X, Y, n_neighbors=10):.4f}  {t_ref:.4f}  {w:+.3f}   {gap:5.2f}")

# ---------------- B: three clusters with known sizes and distances ----------------
print("B: perp seed  size B/A  C/A   dist AC/AB  BC/AB   P mass A-B  A-C")
for seed in SEEDS:
    rng = np.random.default_rng(seed)
    X, lab, T, G = cluster_images(rng)
    if seed == 0:
        print(f"B check: max|T T' - I| = {np.abs(T @ T.T - np.eye(3)).max():.1e}, "
              f"max|G G' - I| = {np.abs(G @ G.T - np.eye(2)).max():.1e}")
    D2 = np.sum((X[:, None] - X[None]) ** 2, -1)
    for u in [None] + PERPS:
        Y = X if u is None else TSNE(perplexity=u, random_state=seed).fit_transform(X)
        cen = np.array([Y[lab == k].mean(0) for k in range(3)])
        rms = [np.sqrt(np.mean(np.sum((Y[lab == k] - cen[k]) ** 2, 1))) for k in range(3)]
        dAB, dAC, dBC = (np.linalg.norm(cen[a] - cen[b]) for a, b in [(0, 1), (0, 2), (1, 2)])
        line = (f"   {'true' if u is None else u:>4} {seed:4d}  {rms[1] / rms[0]:7.3f} {rms[2] / rms[0]:7.3f}"
                f"  {dAC / dAB:8.3f} {dBC / dAB:7.3f}")
        if u is not None:
            P = squareform(_joint_probabilities(D2.astype(np.float32), u, 0))
            line += f"   {2 * P[np.ix_(lab == 0, lab == 1)].sum():.2e}  {2 * P[np.ix_(lab == 0, lab == 2)].sum():.2e}"
        print(line)

# ---------------- C: the mechanics against scikit-learn ----------------
rng = np.random.default_rng(0)
th = rng.uniform(0, 180, 360)
X = bar_images(th, rng)
D2 = np.sum((X[:, None] - X[None]) ** 2, -1)
for u in PERPS:
    Pn, perr = p_numpy(D2, u)
    Ps = squareform(_joint_probabilities(D2.astype(np.float32), u, 0))
    print(f"C1 perp {u:3d}: max|P_numpy - P_sklearn| / max P = {np.abs(Pn - Ps).max() / Ps.max():.1e}"
          f"  (tol 1e-04); row perplexity error {perr:.1e} (tol 1e-08)")
P, _ = p_numpy(D2, 30)
Y = np.random.default_rng(10).normal(size=(360, 2))
C, g = kl_grad(P, Y)
idx, h, fd = np.random.default_rng(11).choice(720, 20, replace=False), 1e-5, []
for k in idx:
    Yp, Ym = Y.ravel().copy(), Y.ravel().copy()
    Yp[k] += h
    Ym[k] -= h
    fd.append((kl_grad(P, Yp.reshape(-1, 2))[0] - kl_grad(P, Ym.reshape(-1, 2))[0]) / (2 * h))
print(f"C2: max|finite difference - Eq. 5| / max|grad| = {np.abs(np.array(fd) - g.ravel()[idx]).max() / np.abs(g).max():.1e} (tol 1e-06)")
C_sk, g_sk = _kl_divergence(Y.ravel().copy(), squareform(P, checks=False), 1, 360, 2)
print(f"C3: KL {C:.6f} vs sklearn {C_sk:.6f}, rel diff {abs(C - C_sk) / C:.1e}; "
      f"gradient rel diff {np.abs(g.ravel() - g_sk).max() / np.abs(g).max():.1e} (tol 1e-12)")

Output:

A check: max|img(37)-img(217)| = 6.1e-16; distance 0-179 deg 0.292, 0-90 deg 6.584
A: perp seed   T(k=10)  T_ref   winding  gap
   true    0   0.9995  0.9995  +1.000   11.04
      5    0   0.9986  0.9995  -1.000   56.05
     30    0   0.9996  0.9995  -1.000    4.05
    100    0   0.9994  0.9995  +1.000    4.20
   true    1   0.9996  0.9996  +1.000    9.47
      5    1   0.9985  0.9996  +0.000   71.49
     30    1   0.9997  0.9996  +1.000    4.69
    100    1   0.9994  0.9996  +1.000    3.56
   true    2   0.9996  0.9996  +1.000    8.91
      5    2   0.9985  0.9996  +0.000   64.09
     30    2   0.9997  0.9996  -1.000    3.44
    100    2   0.9994  0.9996  +1.000    3.20
B: perp seed  size B/A  C/A   dist AC/AB  BC/AB   P mass A-B  A-C
B check: max|T T' - I| = 4.4e-16, max|G G' - I| = 5.6e-15
   true    0    4.451  16.785     5.004   5.007
      5    0    1.352   1.058     1.004   0.962   2.84e-12  2.84e-12
     30    0    1.016   0.981     1.042   1.149   2.84e-12  2.84e-12
    100    0    1.353   0.454     2.291   2.298   4.20e-02  8.44e-03
   true    1    4.205  15.907     5.012   5.013
      5    1    0.923   1.033     1.000   0.947   2.84e-12  2.84e-12
     30    1    1.028   1.008     1.062   1.009   2.84e-12  2.84e-12
    100    1    1.303   0.186     2.291   2.297   4.20e-02  8.43e-03
   true    2    4.610  17.073     5.005   5.005
      5    2    0.809   1.076     0.936   1.882   2.84e-12  2.84e-12
     30    2    1.006   1.000     0.972   1.011   2.84e-12  2.84e-12
    100    2    1.370   0.463     2.293   2.294   4.20e-02  8.43e-03
C1 perp   5: max|P_numpy - P_sklearn| / max P = 6.2e-06  (tol 1e-04); row perplexity error 1.0e-10 (tol 1e-08)
C1 perp  30: max|P_numpy - P_sklearn| / max P = 9.6e-06  (tol 1e-04); row perplexity error 1.0e-10 (tol 1e-08)
C1 perp 100: max|P_numpy - P_sklearn| / max P = 5.6e-06  (tol 1e-04); row perplexity error 9.9e-11 (tol 1e-08)
C2: max|finite difference - Eq. 5| / max|grad| = 1.3e-08 (tol 1e-06)
C3: KL 2.727517 vs sklearn 2.727517, rel diff 3.3e-14; gradient rel diff 2.1e-16 (tol 1e-12)

The plan fixed four hypotheses and the checks’ tolerances before the script was first run.

The checks held. The data have the planned properties: a bar at 37° and at 217° gives the same image to rounding, the images at 0° and 179° are more than twenty times closer than those at 0° and 90°, the templates and ramps are orthonormal, and the sampled clusters have size ratios near 4 and 16 and a distance ratio near 5 in pixel space. Check C1 held: NumPy’s $P$ agrees with scikit-learn’s to within about $10^{-5}$ of the largest entry at all three perplexities, consistent, in our reading, with the library’s float32 search and its entropy tolerance of $10^{-5}$ nats, while the NumPy rows reach their perplexity to about $10^{-10}$ . Check C2 held: Eq. 5 matches the finite differences of the cost. Check C3 held: scikit-learn’s exact objective gives the same cost and gradient to rounding, although it stores only the $n(n-1)/2$ pairs $i < j$ , so its condensed $P$ sums to $1/2$ and its cost carries a factor of 2.

Neighbourhoods were kept (H1 held). The hypothesis was that at perplexities 30 and 100 the map’s trustworthiness is at least that of the true circle coordinates minus 0.02. Every such map scored between 0.9994 and 0.9997, within 0.0002 of the true coordinates’ 0.9995–0.9996, and at perplexity 30 above them on all three seeds. That a map can exceed the true coordinates follows from what is measured: trustworthiness compares with the neighbours of the noisy images, which t-SNE fits, not with the neighbours in $\theta$ . At perplexity 5, where nothing was predicted, the scores were 0.9985–0.9986.

The loop was kept at perplexities 30 and 100 (H2 held). The hypothesis was a winding number of $\pm 1$ for every seed at these perplexities, and all six maps gave it, with a largest step of 3.2 to 4.7 times the median, smaller than the true coordinates’ 8.9 to 11.0 because uniformly drawn angles leave uneven gaps on the true circle. At perplexity 5, for which nothing was predicted, the largest step was 56 to 71 times the median on all three seeds and the winding number was 0 for seeds 1 and 2: the loop was cut, and in two maps its pieces no longer went around a common centre. This is Wattenberg and colleagues’ point that topology may need more than one perplexity to read: here a map at one low perplexity would suggest an open curve or separate pieces where the data form a closed loop, although its trustworthiness stayed near 0.999, because each point’s few nearest neighbours were still right.

Cluster sizes were equalized at perplexities 5 and 30 but not at 100 (H3 held at 5 and 30, failed at 100). The hypothesis was that both map size ratios lie in $[0.5, 2]$ , against true ratios of about 4 and 16, at every perplexity and seed. At perplexities 5 and 30 all ratios lay between 0.809 and 1.352: a cluster sixteen times larger than another in pixel space appeared about the same size. This follows from the construction of $P$ : each $\sigma_i$ is searched until $P_i$ has $u$ effective neighbours, and with $u$ below the 80 images of a cluster that search scales each cluster’s neighbourhoods to its own spread, so $P$ for the three clusters differs little. At perplexity 100 H3 failed: C/A was 0.454, 0.186 and 0.463, below 0.5 on all three seeds, while B/A stayed inside the range (1.30–1.37). The largest cluster in pixel space came out the smallest in the map. Wattenberg and colleagues report the same appearance for one of three Gaussian clusters at perplexity 100.

Distances between clusters were not kept (H4 held). The hypothesis was that the map’s $d_{AC}/d_{AB}$ lies outside $[4, 6]$ around the true 5 at every perplexity and seed. At perplexities 5 and 30 it was 0.936–1.062: the three clusters were roughly equidistant. The printed $P$ mass, computed with the exact _joint_probabilities as in check C1, shows why. At these perplexities it is $2.84 \times 10^{-12}$ between A and B and between A and C, which is $2 \times 80 \times 80$ times $2.22 \times 10^{-16}$ : every cross-cluster $p_{ij}$ of the exact $P$ sits at the floor to which that function clips probabilities, so the exact cost would receive the same, negligible, information about the near pair A–B as about the far pair A–C. The maps themselves were fitted with the default Barnes-Hut method, whose sparse $P$ over each image’s $\lfloor 3u + 1 \rfloor$ nearest neighbours (16 and 91 at these perplexities) has no floor; in our reading it carries no more information about the cluster distances, since its cross-cluster entries are either absent or the same tiny Gaussian weights left unclipped. At perplexity 100 the neighbour lists hold all 239 other images, and the mass is 0.042 between A and B and 0.0084 between A and C, five times less, and the ratio rises to 2.29 on all three seeds, so it reflects the order of the distances but is still less than half the true value. Wattenberg and colleagues found that with 50 points per cluster a perplexity of 50 showed the global geometry and with 200 points per cluster none of their perplexities did; here no perplexity kept the distances, and the one that came closest, 100, lost the sizes.

Process. The four hypotheses and the 0.02 margin of the first were written before any script ran; a planning script that fitted no map then computed the reachable values, namely the true coordinates’ trustworthiness, which set the first hypothesis’s bars of 0.9775–0.9776, and the agreement of the checks, which set their tolerances at ten times the reached value rounded up to a power of ten. A pilot of the full script followed; the editor approved it unchanged, and the output shown is the first run, identical to the pilot. The plan departs from the batch’s description of the example in eight ways: perplexities 5, 30 and 100 were chosen; the winding number and largest step were added for the loop; the true coordinates’ trustworthiness was added as the reference; block B realizes sizes as illumination-ramp variation in the ratio 1 : 4 : 16 and distances as orthogonal templates with contrasts 1, 1 and 7, without pixel noise; the cross-cluster $P$ mass was added; check C3 against scikit-learn’s exact objective was added; the seed draws the data as well as the map’s seed; and “size ratios near 1” and “distance ratios not preserved” were made operational as the ranges $[0.5, 2]$ and outside $[4, 6]$ .

Properties and Behavior

Local structure first. The KL cost penalizes separating true neighbours heavily and joining non-neighbours lightly; with the heavy-tailed kernel, moderately distant points are pushed far out. Trustworthiness, which measures exactly what the cost favours, is high, as in the example; it cannot confirm that sizes, distances or global shape are faithful.

Density is equalized. Because every point gets the same perplexity, dense regions are expanded and sparse ones contracted. Wattenberg and colleagues show a cluster ten times as dispersed as another drawn the same size, and two parallel lines bending outwards because their middles are denser than their ends; they call the equalization a predictable feature, by design, rather than an ordinary distortion. The Embeddings article lists cluster sizes and gaps among the things a map does not show; its misconception that well-separated t-SNE clusters prove an embedding separates the classes is answered by measuring in the embedding space itself, for example with a nearest-neighbour classifier or a linear probe.

Distances between clusters depend on perplexity and cluster size. When the perplexity is small compared with the clusters, $P$ links them only by numbers at or below machine precision, as in the example, and their arrangement in the map is set by the optimization rather than by the data. Wattenberg and colleagues conclude that distances between well-separated clusters “may mean nothing”.

Noise can look like structure. For 500 points from a 100-dimensional unit Gaussian, Wattenberg and colleagues found apparent clusters at perplexity 2, and, at perplexity 30, an evenly spread cloud that they read as faithful to the near-spherical shell on which such points lie (the shell is discussed in Latent Space).

Non-convexity and runs. The cost is not convex, so maps can differ between runs and settings (2008, §6.2). The authors report that one choice of optimization parameters worked across their visualization tasks and that the quality of the optima did not vary much from run to run. scikit-learn’s user guide says it is legitimate to try several seeds and keep the map with the lowest KL divergence. Wattenberg and colleagues show five runs of a trefoil knot at perplexity 2 giving three different broken shapes as well as two circles, while five runs at perplexity 50 looked the same up to symmetry. The map’s KL divergence compares runs at one perplexity only, since $P$ changes with $u$ .

Scale of the map. In SNE, changing the variance of the map’s Gaussian only rescales the final map (2008, footnote 2); in our reading, t-SNE’s map coordinates likewise have no units, and the absolute spread of a map, and its axes, carry no meaning.

Limitations

The 2008 paper lists three weaknesses (§6.2):

General dimensionality reduction. It is unclear how t-SNE performs when reducing to more than three dimensions; the authors expect heavier tails to preserve local structure less well there, and suggest more degrees of freedom.
The curse of intrinsic dimensionality. t-SNE relies on local Euclidean distances, so on data with high intrinsic dimension and a highly varying manifold its implicit assumption of local linearity can fail. They suggest running t-SNE on a representation learned by a deep autoencoder, and note that intrinsically high-dimensional data cannot be fully represented in two or three dimensions at all.
Non-convex cost, with solutions that depend on optimization parameters and may differ between runs from random starts.

Further limitations follow from the method’s definition: t-SNE gives no mapping for points not in the fitted set (the paper names a parametric version as future work, and scikit-learn’s TSNE has fit_transform but no transform); the exact method is quadratic in $n$ ; and cluster sizes, inter-cluster distances, and the density of points are not to be read from a map, as the example and Wattenberg et al. show. Labelled colours on a map are a picture of neighbourhoods, not a test of a representation.

Variants

SNE (Hinton & Roweis, 2002), as van der Maaten and Hinton describe it: conditional distributions in both spaces, a Gaussian in the map, and optimization with momentum and annealing noise.
Symmetric SNE and UNI-SNE, described above through the 2008 paper.
Random-walk t-SNE (2008, §5): affinities between landmark points from random walks on a neighbourhood graph of the whole data set, so the undisplayed points shape the map.
Barnes-Hut and dual-tree t-SNE (van der Maaten, 2014): sparse nearest-neighbour affinities and tree approximations of the repulsive forces, $O(N \log N)$ per iteration; the Barnes-Hut version is scikit-learn’s default.
Heavier or lighter tails: scikit-learn’s choice of $\max(d - 1, 1)$ degrees of freedom for a $d$ -dimensional map follows van der Maaten’s parametric t-SNE paper (2009), as its source comment says.
UMAP, a neighbour-graph method with a different theory and optimization, is compared with t-SNE in Embeddings.

Use in Computer Vision

t-SNE maps are used to inspect learned image representations rather than to compute with them. The 2008 paper’s experiments already included images: on 6,000 MNIST digits t-SNE separated the classes almost perfectly, while Isomap and LLE showed large overlaps between classes and Sammon mapping separated only three of them somewhat, and on COIL-20, 72 views of each of 20 objects, it drew many objects’ viewpoint manifolds as closed loops, as the bars above form one. In later work, maps of network activations compare layers or training methods: Deep Features describes DeCAF’s plots, in which the first pooling layer showed no clear semantic clustering and a fully connected layer a clear one; SimCLR’s appendix compares the representation before and after the projection head; Embodied AI reports a map showing that simulated and real images of the same places fall into separate clusters; and the Nature deep-learning review shows a t-SNE projection of learned word vectors. Each of these is a qualitative picture; the quantities behind such claims are measured in the original space.