t-SNE, in depth: perplexity calibration, the Student-t map, the gradient, Barnes-Hut and what the plot does not show

Input similarities and perplexity

Start in the input space. For each point i, define how likely i is to pick j as a neighbour with a Gaussian centred on i:

p(j | i) = exp(-|x_i - x_j|^2 * beta_i) / sum_{k != i} exp(-|x_i - x_k|^2 * beta_i)
beta_i   = 1 / (2 * sigma_i^2)

The width sigma_i is different for every point, and that is the key design choice. It is chosen so that the distribution has a fixed perplexity, two raised to its entropy in bits, or equivalently entropy log(perplexity) in nats. Perplexity is a smooth count of effective neighbours: perplexity 30 means each point spreads its attention over about 30 others, whether it sits in a dense region or a sparse one. A binary search on beta per point hits the target entropy.

The conditional distributions are then symmetrised into one joint distribution, p_ij = (p(j|i) + p(i|j)) / (2N), which sums to 1 over all pairs and guarantees that outliers, whose rows would otherwise be ignored, still pull on someone.

The map: a heavy-tailed kernel

In the map, each point i has a 2-D position y_i. Similarity there uses a Student-t kernel with one degree of freedom, normalised over all pairs:

q_ij = (1 + |y_i - y_j|^2)^-1 / sum_{k != l} (1 + |y_k - y_l|^2)^-1

Why the heavy tail? A 2-D map has far less room than a 50-D space. Points at a moderate distance in the input cannot all fit at moderate distances in the plane, so a Gaussian map kernel squashes them together, the crowding problem that limited the earlier SNE method. The t kernel decays as a power law, so a moderately dissimilar pair can be placed much further apart in the map at little cost. That opens gaps between clusters, which is why t-SNE plots look so well separated.

The objective is the Kullback-Leibler divergence KL(P || Q) = sum p_ij log(p_ij / q_ij). KL is asymmetric in a way that matters: a large p_ij with small q_ij (neighbours torn apart) is expensive, while a small p_ij with large q_ij (strangers placed together) costs little. t-SNE protects neighbours first and treats long distances as a soft afterthought.

The gradient and a working implementation

The gradient has a clean form that reads like a spring system:

dKL/dy_i = 4 * sum_j (p_ij - q_ij) * (y_i - y_j) / (1 + |y_i - y_j|^2)

Each pair either attracts (p_ij greater than q_ij) or repels (p_ij smaller). The attraction only involves pairs with meaningful p_ij, which are sparse. The repulsion involves every pair through the normaliser of Q, and that all-pairs term is where the cost of t-SNE lives.

This is an exact implementation, small enough to read in one sitting, using the same schedule as scikit-learn 1.8.0: perplexity 30, early exaggeration 12 for the first 250 iterations with momentum 0.5, then momentum 0.8 to iteration 1,000, a learning rate of max(N / 12 / 4, 50), per-coordinate adaptive gains, and a PCA initialisation rescaled so its first coordinate has standard deviation 1e-4.

import numpy as np

def sq_dists(X):
    s = (X * X).sum(1)
    D = s[:, None] + s[None, :] - 2 * X @ X.T
    np.maximum(D, 0, out=D)
    return D

def calibrate(D, perplexity, tol=1e-5, iters=64):
    """Row-wise binary search on beta = 1/(2 sigma^2) so H(P_i) = log(perplexity)."""
    n = D.shape[0]
    P = np.zeros((n, n))
    target = np.log(perplexity)
    betas = np.ones(n)
    for i in range(n):
        d = np.delete(D[i], i)
        lo, hi, beta = 0.0, np.inf, 1.0
        for _ in range(iters):
            w = np.exp(-(d - d.min()) * beta)
            p = w / w.sum()
            H = -(p * np.log(np.maximum(p, 1e-300))).sum()
            if abs(H - target) < tol:
                break
            if H > target:            # too flat: sharpen
                lo = beta
                beta = beta * 2 if hi == np.inf else (beta + hi) / 2
            else:
                hi = beta
                beta = (beta + lo) / 2
        P[i, np.arange(n) != i] = p
        betas[i] = beta
    return P, betas

def tsne(X, perplexity=30.0, iters=1000, exag=12.0, exag_iters=250, seed=0):
    n = X.shape[0]
    P, betas = calibrate(sq_dists(X), perplexity)
    P = (P + P.T) / (2 * n)                       # symmetric joint p_ij, sums to 1
    P = np.maximum(P, 1e-12)
    lr = max(n / exag / 4, 50)                     # the "auto" rule
    Xc = X - X.mean(0)                             # PCA init, first coord std 1e-4
    U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
    Y = Xc @ Vt[:2].T
    Y = Y / Y[:, 0].std() * 1e-4
    vel = np.zeros_like(Y)
    gains = np.ones_like(Y)
    for t in range(iters):
        Pt = P * exag if t < exag_iters else P
        mom = 0.5 if t < exag_iters else 0.8
        W = 1.0 / (1.0 + sq_dists(Y))              # Student-t kernel, one degree of freedom
        np.fill_diagonal(W, 0.0)
        Q = np.maximum(W / W.sum(), 1e-12)
        G = 4.0 * ((Pt - Q) * W)                   # n x n coefficient matrix
        grad = G.sum(1)[:, None] * Y - G @ Y       # sum_j G_ij (y_i - y_j)
        same = np.sign(grad) == np.sign(vel)
        gains = np.where(same, gains * 0.8, gains + 0.2).clip(0.01)
        vel = mom * vel - lr * gains * grad
        Y = Y + vel
    return Y, betas

The optimisation schedule

Early exaggeration. Multiplying P by 12 for the first 250 iterations makes attraction dominate. Points that belong together collapse into tight groups while the map is still small and they can move freely; after exaggeration ends, the groups spread out and settle. In the run below the KL divergence was 1.91 at iteration 250 and 1.17 at iteration 1,000. KL during exaggeration is not comparable with KL after it, so judge convergence on the second phase only.

Learning rate. The old fixed default of 200 is too small for large N: the map is still expanding when the iteration budget runs out, and clusters look fragmented. scikit-learn's learning_rate="auto" scales with N, which is the practical fix. If a big map looks like scattered fragments, raise the learning rate or iterations before blaming the data.

Initialisation. Random initialisation lets the global arrangement of clusters land anywhere; different seeds give different layouts. PCA initialisation, the current scikit-learn default, starts the map in the data's main directions, so the coarse layout is more stable from run to run and carries a little more global structure.

Scaling past a few thousand points

The exact code builds N by N matrices, which is fine up to a few thousand points and hopeless beyond. Two approximations make t-SNE practical:

• Sparse P. Only about 3 times perplexity nearest neighbours carry non-negligible p(j|i). scikit-learn computes exactly min(N - 1, int(3 * perplexity + 1)) neighbours per point. Approximate nearest-neighbour indexes such as the structures in the live HNSW article make this step fast for large N.
• Barnes-Hut repulsion. Build a quadtree over the current map; a distant cell whose angular size is below a threshold is treated as one point at its centre of mass. This gives O(N log N) per iteration. scikit-learn's angle defaults to 0.5, and its documentation warns that error grows quickly above 0.8. Barnes-Hut in scikit-learn only supports up to 3 output dimensions. The live quadtree and octree article explains the tree.
• FFT interpolation. FIt-SNE and openTSNE interpolate the repulsive field on a grid and use FFTs, which scales roughly linearly in N and handles millions of points.

A sensible pipeline for wide data first reduces to around 50 dimensions with PCA, which removes noise and speeds up the neighbour search, then runs t-SNE on those scores.

Worked example: what the map distorts, measured

Here is the experiment that should change how you read t-SNE plots. Three Gaussian clusters in 50 dimensions: A has 200 points with standard deviation 1, B has 200 points with standard deviation 6 centred at +30 on every axis, and C has 100 points with standard deviation 1 at -30. Running the code above:

Cluster B is six times wider than A in the input and the same size in the map. That is not a bug: the per-point sigma adapted to B's density (10.70 against 1.82) precisely so every point sees about 30 neighbours, which equalises apparent density. Cluster C is twice as far from B as from A in the input, yet all three centroid distances in the map are within 8% of each other. The heavy tail only needs clusters to be far apart, not proportionally far apart.

Now the opposite test: 300 points of pure 10-D Gaussian noise, which has no clusters at all. Run k-means with k = 4 on the result. The silhouette score is 0.081 on the input, 0.386 on the perplexity-2 map and 0.363 on the perplexity-30 map. A plain 2-D Gaussian sample scores 0.313 under the same test, so part of that jump is just the drop to two dimensions, but the lesson stands either way: cluster-quality metrics computed on t-SNE coordinates measure the map, not the data.

Failure modes

• Reading sizes and gaps. As measured above, cluster area and between-cluster distance are not interpretable. Report them from the input space.
• Clustering on the map. Running k-means or DBSCAN on t-SNE output finds structure the optimiser helped create. Cluster in the input space (see the live k-means article) and use t-SNE only to colour and inspect the result.
• One perplexity. Low perplexity fragments real clusters into islands; high perplexity merges small groups. Run at least two values and trust only structure that survives both.
• Too few iterations. A map that is still expanding looks like disconnected fragments. Check that KL has flattened after exaggeration.
• Duplicates and unscaled features. Exact duplicate points make distances zero and break calibration for those rows; one feature in large units dominates every distance. Deduplicate and standardise first.
• Expecting a transform. Classic t-SNE has no function for new points; scikit-learn's TSNE offers only fit_transform. openTSNE can embed new points into an existing map. For a reusable projection, use PCA or a parametric model.
• Comparing runs pixel by pixel. Rotation, reflection and cluster placement change across seeds. Fix the seed and the initialisation for reproducible figures.

Trade-offs against other reducers

Choose t-SNE when the question is whether the representation separates things you care about, such as embeddings by class or cells by type. Do not choose it when you need to measure how far apart groups are, or to project new data in production.

What to do next

1. Standardise features, remove duplicates, and reduce to about 50 dimensions with PCA.
2. Run scikit-learn's TSNE with defaults (perplexity 30, PCA init, auto learning rate) and a fixed random_state; use openTSNE or FIt-SNE beyond about 100,000 points.
3. Repeat at a low and a high perplexity, for example 10 and 50, and keep only structure that appears in both.
4. Colour the map by known labels or by clusters computed in the input space, never by clusters computed on the map.
5. Report sizes and distances from the input space; caption the plot so readers know they are not meaningful in it.
6. Check that KL flattened after the exaggeration phase before publishing a figure.
7. Use PCA, UMAP's transform or a parametric model if new points must be placed later.