Why t-SNE Preserves Neighborhoods but Distorts Distances

t-SNE produces a two-dimensional embedding by defining a probability distribution over pairs of points in the original space, defining a second distribution over pairs in the embedding, and moving the embedded points until the two distributions are as close as one asymmetric divergence measure allows. Everything distinctive about the method, including the shapes of the plots it produces and the fact that distances between clusters carry no meaning, comes from the choice of those two distributions and of the divergence.

Two similarity distributions

For each point i, the high-dimensional distribution over its neighbours is a Gaussian centred on it,

$$p_{j|i} = \frac{\exp(-||x_i - x_j||^2 / (2 \sigma_i^2))}{\sum_{l \ne i} \exp(-||x_i - x_l||^2 / (2 \sigma_i^2))}$$

with the diagonal set to zero so that a point is not its own neighbour. Each point gets its own bandwidth sigma_i, fixed by a constraint described in the next section; the asymmetry between points is only removed at the end, by symmetrizing p_ij = (p_{j|i} + p_{i|j}) / (2n), at which point each pair has one probability used by both endpoints.

The low-dimensional distribution uses a different kernel, the Student t distribution with one degree of freedom,

$$q_{ij} = \frac{(1 + ||y_i - y_j||^2)^{-1}}{\sum_{l \ne m} (1 + ||y_l - y_m||^2)^{-1}}$$

and this is the choice that makes the picture readable. A Gaussian kernel in two dimensions would have to place a moderate number of points at moderate distances inside a fixed radius, but the area available at radius r grows only as r^2, whereas the volume of a ball in the original space grows as r^d. Matching the high-dimensional probabilities exactly would force the embedding to collapse toward the centre, the crowding problem the original presentation names. A kernel with a heavy tail assigns a more slowly decaying probability to distant pairs, which supplies enough repulsion to keep the embedding spread out. The single degree of freedom means there is no scale parameter to tune in the low-dimensional kernel.

Perplexity sets the neighbourhood size

The bandwidth sigma_i is not chosen by hand. It is found by binary search so that the conditional distribution has a prescribed perplexity, 2^H(P_i), where H is the Shannon entropy of that distribution in bits. Perplexity can be read as an approximate effective number of neighbours: a value of 30 means the Gaussian at that point is wide enough that roughly thirty neighbours carry appreciable mass. The search decouples the points, and it has no closed-form solution, so an embedding costs one small root find per point before any gradient step.

Values are usually chosen between 5 and 50, and implementations commonly default to 30. The choice changes the result substantially. A small perplexity produces many small, tight groups, because each point only wants its closest handful of neighbours nearby; a large perplexity produces a smoother map in which local structure is absorbed into larger masses. There is no value of perplexity that shows all structure at once, and inspecting several is part of using the method. Because each point’s bandwidth is fitted to its own neighbourhood distances, a dense region and a sparse region of the same dataset get kernels of different widths, which is a form of adaptivity that a single global radius cannot provide.

The objective is deliberately asymmetric

The embedding minimizes the Kullback-Leibler divergence from the high-dimensional to the low-dimensional distribution,

$$C = \sum_{i \ne j} p_{ij} \log \frac{p_{ij}}{q_{ij}}$$

and the direction of the divergence matters. The term for a pair contributes p log(p/q); when p is large and q is small the penalty is large and the points are pulled together, while when p is small and q is large the contribution is small because p multiplies it. The objective therefore concentrates on keeping true neighbours close and barely rewards pushing apart pairs that were never close to begin with. The gradient is a sum of attractive and repulsive terms weighted by p_ij - q_ij, with both scaled by the kernel’s derivative, so pairs that are currently too close relative to their target push harder than pairs that are too far.

That asymmetry is the reason the method is good at revealing local structure and the reason the resulting map cannot be read as a metric. It also means the global arrangement of the groups in the picture is weakly determined: which cluster ends up near which, and how large the gaps between them are, follow from the initialization and the optimization path rather than from the original distances.

Initialization, exaggeration, and why runs differ

The objective is non-convex, and the early iterations are dominated by repulsion, so a random starting configuration will spread points apart before the attractive terms take effect. Initializing from the leading linear components, as produced by Principal Component Analysis Through the Singular Value Decomposition, gives a deterministic start and preserves the coarse linear structure; initialization from another embedding, or from a random sample, is also used. Given the same initialization and the same parameters, the gradient descent is deterministic apart from numerical details, which is why implementation options for specifying a seed exist.

Early exaggeration multiplies all p_ij by a factor greater than one, commonly 4, for an initial phase of the optimization, which makes every point claim its neighbours much more strongly than it will later. The effect is that temporary clusters form quickly and nearly disjointly, and when the factor returns to one, the clusters relax into their final positions with wider gaps between them. The schedule is part of the algorithm’s output, not a reproduction of the data.

Optimization itself is momentum gradient descent, and the learning rate and momentum schedule both influence where the descent lands. The embedding is a local minimum of the divergence, and different local minima are different trades between how much of the neighbourhood structure is preserved in each region; two maps of the same data with the same perplexity and different seeds are both valid under the objective.

What the picture does and does not show

Distances between clusters are not data distances. They depend on the initialization and the exaggeration schedule, and the divergence barely penalizes placing two dissimilar groups near each other. The size of a cluster in the plot is not the number of points in it in any absolute sense, since the low-dimensional kernel has no scale parameter and cluster areas are not calibrated to anything: a group of fifty points and a group of five hundred can occupy similar areas. The density of points inside a cluster is not a density estimate either; the kernel flattens the differences that a density plot would show.

Structure that is absent from the plot is not necessarily absent from the data. Two groups can be merged by the optimization when the perplexity is too large, and a single group can be split into fragments when it is too small. Since a run yields one local minimum, comparing several runs is the cheapest way to separate stable structure from optimization artefact. The neighbourhood-graph machinery behind the method also connects it to Spectral Clustering and the Relaxation of Normalized Cut, where the same affinities are used to cut a graph rather than to position points.

The classical formulation has no out-of-sample mapping. The embedding is a set of coordinates fitted to one dataset, and there is no function to apply to a new point without re-running the optimization, because the objective couples every point to every other point through the normalizing constant of q. Variants that train a parametric map, and methods that approximate the position of new points by optimizing their coordinates against a frozen embedding, exist, but the standard use of t-SNE is a transformation of one fixed dataset into a picture.

Cost and the alternatives

Computing q naively requires all O(n^2) pair terms per iteration, both in time and in memory, which is why exact t-SNE does not scale past a few thousand points. The Barnes-Hut approximation computes the repulsive contribution with a spatial tree in the embedding space, evaluating distant groups of points as single summaries, which reduces the per-iteration cost to approximately O(n log n). Interpolation-based variants approximate the repulsion on a fixed grid with an FFT, removing the tree traversal from the inner loop. The attractive terms already involve only the neighbours with non-negligible p.

The alternative most often compared is UMAP, which builds a weighted neighbour graph, treats the edge weights as a fuzzy topological representation, and optimizes a cross-entropy between that graph and a low-dimensional analogue with a different family of low-dimensional weights. The two methods share the neighbourhood-graph foundation and the local-first bias, and they differ in the objective, the low-dimensional kernel family, and the amount of global structure the optimization tends to retain. Neither produces a metric embedding, and the honest use of both is as a visual index into structure that a separate, quantitative method then confirms. For group labels rather than a picture, the density-based definition in DBSCAN, Density Reachability, and the Parameter Cliff answers the question directly, and the purely probabilistic alternative is to model the density, as in Kernel Density Estimation and the Bandwidth Trade-off.

import numpy as np

def conditional_probs(D2, perplexity, tol=1e-5, max_iter=60):
    """D2: n by n squared distances. Rows of the result sum to 1."""
    n = len(D2)
    target = np.log(perplexity)          # target entropy of each row, in nats
    P = np.zeros((n, n))
    for i in range(n):
        d2 = D2[i].astype(float).copy()
        d2[i] = np.inf                   # a point is not its own neighbour
        lo, hi = 1e-20, 1e5              # bounds on beta = 1 / (2 sigma^2)
        beta = hi
        for _ in range(max_iter):
            beta = 0.5 * (lo + hi)
            w = np.exp(-beta * (d2 - d2.min()))
            w[i] = 0.0
            p = w / w.sum()
            H = -(p * np.log(np.maximum(p, 1e-300))).sum()
            if abs(H - target) <= tol:
                break
            if H > target:               # too diffuse: widen the kernel
                lo = beta
            else:
                hi = beta
        P[i] = p
    return P

The remaining honest statement is about purpose. A t-SNE plot is a hypothesis about which observations resemble one another, produced under one neighbourhood size and one local minimum. Confirming the hypothesis requires a method whose output is a number, such as the validity measures discussed in Cluster Validity, Stability, and Choosing the Number of Clusters, applied to the original representation rather than to the coordinates of the picture.