Book contents
Contents
raw Math
RAW Book Machine Learning Dimensionality Reduction

Introduction to t-Distributed Stochastic Neighbor Embedding

Robert Eisele

Imagine that every photograph, sound recording, or scientific measurement could be placed somewhere on a map. Similar objects would appear near one another; different ones would drift apart. The map would not reproduce the objects themselves. It would reveal the relationships between them.

The difficulty is that the objects we want to understand frequently have hundreds or thousands of numerical properties, while a map has only two or three coordinates. Something must be lost. The interesting question is therefore not how to preserve everything, but what we should choose to preserve.

This question leads, step by step, to t-distributed Stochastic Neighbor Embedding, usually shortened to t-SNE. Rather than introducing its equations as a collection of rules, we will construct them from the problems they are meant to solve.

An Object Is a Point Once We Choose Its Coordinates

Consider a grayscale image that contains 28 by 28 pixels. Each pixel can be assigned a brightness between zero and one. If we write those brightness values in a long row, the image becomes a vector of 784 numbers. A second image becomes another vector of 784 numbers. What we would ordinarily call a collection of images is, mathematically, a cloud of points in a space with 784 coordinate axes.

The same idea works beyond images. A sensor reading might contain acceleration, temperature, pressure, and hundreds of frequency measurements. A cell might be described by the expression levels of thousands of genes. A piece of text might be represented by a learned vector of numerical features. The meaning of a coordinate changes, but the geometric representation remains the same.

If there are \(n\) observations and each has \(d\) numerical features, we write the input data as

\[\mathcal X=\{\mathbf x_1,\ldots,\mathbf x_n\},\qquad \mathbf x_i\in\mathbb R^d.\]

Here \(\mathbf x_i\) denotes one observation, not one coordinate. It is a whole vector. The subscript \(i\) simply tells us which observation we mean.

We would like to replace each \(\mathbf x_i\) with a smaller vector \(\mathbf y_i\) that we can actually display:

\[\mathcal Y=\{\mathbf y_1,\ldots,\mathbf y_n\},\qquad \mathbf y_i\in\mathbb R^s,\qquad s\in\{2,3\}.\]

The most common choice is \(s=2\), because the result can be drawn directly as a scatter plot. \(s=3\) is also used: a three-dimensional plot remains visually accessible when we can rotate the view. It can represent an extra degree of geometric freedom, although a stationary screen still shows it through a two-dimensional projection. The original t-SNE paper explicitly describes both two- and three-dimensional maps. 2

784 coordinates 0.00.10.9...0.4 2D location
Each image is one high-dimensional point. t-SNE assigns it a location in a new two- or three-dimensional map; the drawing does not recover its 784 original coordinates.

Our first temptation might be to delete all but two or three coordinates of the original vectors. But what if the information distinguishing one handwritten digit from another is spread across hundreds of pixels? Removing most coordinates could destroy precisely the relationships we hoped to see.

A more sophisticated possibility is to create new axes as combinations of the old ones. Principal component analysis (PCA) does this by finding directions along which the data varies most. It provides an especially useful linear view, but a curved or folded arrangement cannot always be unfolded by choosing different straight axes. PCA and t-SNE therefore answer different questions: one prioritizes explained variance in a linear projection; the other prioritizes a probabilistic representation of neighborhoods. Neither is universally superior.

Why Many Coordinates Need Not Mean Many Degrees of Freedom

A thousand coordinates do not necessarily imply a thousand independent ways in which an object can change. Consider a camera photographing a toy robot mounted on a turntable. Every frame might contain hundreds of thousands of pixel values. Yet if the lighting, camera, and robot remain fixed, changing the angle of the turntable is enough to produce all the images. In that idealized experiment, a single number—the angle—controls a very high-dimensional observation.

We can express the idea as a function. Let \(\theta\) be the angle of the turntable and let \(\mathbf x(\theta)\in\mathbb R^d\) be the vector of pixel intensities produced at that angle. As \(\theta\) changes smoothly, the image vector traces a curved path through \(d\)-dimensional space. The path can have a single intrinsic degree of freedom even though it lives inside a space with \(d\) coordinate axes. With two independently moving joints, the family of images might instead trace a two-dimensional surface. This is the intuition behind the manifold hypothesis: observations often occupy a much smaller structured subset of all mathematically possible feature combinations. The hypothesis is a useful model, not a universal law. 11

Ambient space: three coordinatesIntrinsic coordinate: one angle curved path0θThis schematic is not a t-SNE result.
A curved family of observations may have one intrinsic parameter despite living in three coordinates.

The distinction between ambient and intrinsic dimension: an object can require many coordinates to describe, but few independent parameters to vary.

This suggests an attractive ambition. If the data really lies near a curved low-dimensional surface, perhaps we can unfold that surface and draw it on a sheet of paper. However, different methods mean different things by unfold. Isomap, for example, estimates long-range distances by following paths through a neighborhood graph in an attempt to recover geodesic structure. t-SNE does not do that. It does not attempt to preserve every distance along a manifold. Instead, it seeks a map in which locally similar observations still tend to be neighbors. 11, 2

There is also a distinction between the straight-line distance through the ambient space and the distance along a curved surface. Two observations on different turns of a spiral can be close in ordinary Euclidean distance even though traveling along the spiral from one to the other would take a much longer route. Which distance is meaningful depends on what we are studying. Local distances are not automatically trustworthy merely because they are local; noise, scaling, irrelevant coordinates, and the representation itself can all distort them.

This gives us a practical rule. Dimensionality reduction does not discover meaning from nothing. It begins with a definition of similarity and asks which part of the resulting geometry can be made visible. Before we search for a low-dimensional map, we must decide what the coordinates in the original space represent and what their distances mean.

The Geometry We Choose Before t-SNE

Imagine a little collection of doodle monsters: one has a single enormous eye and two horns, another has three eyes, and a third looks like an octopus that has forgotten how many tentacles it owns. Each drawing consists of a sequence of strokes: move the pen, draw a curve, lift the pen, and begin somewhere else. Two people could draw essentially the same three-eyed creature with different numbers of strokes and in different orders. Their raw descriptions need not even contain the same number of coordinates. What, then, should their distance be?

We could rasterize both drawings and compare their pixels. That is easy to define mathematically, but the resulting distance would often be sensitive to translation, stroke width, and other details that do not determine whether both drawings depict the same kind of ridiculous creature. We have a problem before dimensionality reduction begins: we have not yet decided what similarity should mean.

One alternative is to learn a common representation. In the original Sketch-RNN research, Ha and Eck designed a recurrent generative model for drawings represented as sequences of pen strokes. An encoder can turn such a variable-length sequence into a fixed-length latent vector. Each coordinate in this vector is learned from data rather than hand-selected as a pixel or a stroke position. [9]

Imagine an encoder that maps each sketch to a vector with 128 numbers. Now every sketch is a point in the same 128-dimensional space. We can compare these points, choose a distance measure, and finally apply PCA or t-SNE to draw a two- or three-dimensional map. But the map is a picture of the learned representation, not a direct picture of the original strokes. If the encoder learned to regard two different drawing styles as similar, those sketches may become neighbors. It might instead emphasize the number of eyes or the shape of the tentacles; only examining the actual learned representations could tell us.

Doodle monstersLearned featuresNeighborhood map z ∈ R¹²⁸ Conceptual drawing: no encoder was trained and no t-SNE output was measured for these monsters.
The crucial step precedes dimensionality reduction. These monsters and their displayed clusters are fictional, not observations from a trained model.

This distinction is especially important when interpreting embeddings of neural networks. A convincing t-SNE plot does not establish that a network has discovered the categories we hoped for. It shows relationships between its internal codes according to the metric we supplied. To understand those relationships, we must return to the original drawings and ask what the nearby examples actually have in common. The same principle applies to medical measurements, text embeddings, and sensor recordings: t-SNE organizes the representation it receives; it does not choose the meaning of its coordinates.

The Information We Cannot Afford to Lose

Let us start with something familiar: the distance between two points. In a plane we obtain it from the Pythagorean theorem. If two points differ by \(\Delta x\) horizontally and \(\Delta y\) vertically, then their squared distance is \((\Delta x)^2+(\Delta y)^2\). There is nothing specifically two-dimensional about this rule. Every additional perpendicular axis contributes one more squared difference.

For vectors with \(d\) coordinates, the resulting Euclidean distance is

\[d_{ij}=\|\mathbf x_i-\mathbf x_j\|=\sqrt{\sum_{r=1}^{d}(x_{ir}-x_{jr})^2}.\]

If we could find new points \(\mathbf y_i\) such that every distance between every pair stayed the same, the new map would be geometrically faithful. But that is generally impossible when we reduce the number of dimensions. A genuinely three-dimensional arrangement cannot be placed in a plane while preserving all its pairwise distances. Lowering the dimensionality removes geometric freedom.

So we must choose. We might try to maintain all distances approximately, or we might give some relationships priority over others. Suppose we are comparing handwritten digits. Knowing that two particular images of the digit 3 are close neighbors is useful. Knowing their exact distance from a rather different 8 is usually less important. We would therefore prefer a map that keeps close neighbors together, even if it distorts some long-range distances.

Original relationshipsMap preserves neighbors
Keeping nearby points together does not require every cluster to retain its shape, orientation, size, or exact distance to other clusters.

This turns the problem around. We no longer demand identical geometry. We ask whether the neighborhood of each observation can be recognized after the transformation. Once this becomes the objective, several decisions that initially look arbitrary begin to follow naturally.

One qualification matters from the beginning: what counts as a neighbor depends on the original features and the distance metric. If one coordinate is measured in thousands and another in fractions, Euclidean distance may mostly reflect the first. Feature scaling and the choice of metric are part of the scientific question, not mere technical housekeeping.

From a Hard Decision to a Soft Neighborhood

Imagine standing at one point \(\mathbf x_i\) and asking which observations are its neighbors. We could choose the five closest points and declare everyone else irrelevant. This gives a sharp boundary: the fifth point is a neighbor, while the sixth is not, even when their distances differ by almost nothing.

A softer description would be more informative. We could say that a point at distance 1 is very plausible as a neighbor, a point at distance 2 is less plausible, and a point at distance 10 is extremely unlikely. Crucially, these preferences should change smoothly as distances change.

What mathematical function could express that idea? We need a weight \(w(r)\) with three properties: it is positive; it is largest at \(r=0\); and it decreases smoothly as \(r\) increases. The Gaussian-shaped function

\[w(r)=\exp(-\beta r^2),\qquad \beta>0,\]

satisfies all three. Squaring the distance makes the weight depend only on separation, not direction, and the exponential changes a squared distance into a positive, rapidly decaying value. The Gaussian is a useful design choice, not a theorem that similarity must have this form.

The parameter \(\beta\) controls how quickly attention fades. It is conventional to write \(\beta=1/(2\sigma^2)\), giving the familiar Gaussian expression. For the selected point \(\mathbf x_i\), and allowing its width to be chosen separately, the weight of a candidate \(\mathbf x_j\) becomes

\[w_{j\mid i}=\exp\!\left(-\frac{d_{ij}^2}{2\sigma_i^2}\right).\]

The vertical bar reminds us that we are considering candidate \(j\) from the viewpoint of point \(i\). So far these are merely positive weights, not probabilities. They do not necessarily add up to one.

To turn them into probabilities, divide each weight by the total weight of all possible neighbors. This elementary normalization is the entire reason for the denominator in the next equation:

\[p_{j\mid i}=\frac{w_{j\mid i}}{\sum_{k\ne i}w_{k\mid i}}=\frac{\exp\!\left(-d_{ij}^2/(2\sigma_i^2)\right)}{\displaystyle\sum_{k\ne i}\exp\!\left(-d_{ik}^2/(2\sigma_i^2)\right)}.\]

We exclude the point itself, so \(p_{i\mid i}=0\). Everything else sums to one: \(\sum_{j\ne i}p_{j\mid i}=1\). We can now read \(p_{j\mid i}\) as the probability of selecting \(j\) from the other observations by a random rule that favors nearer points. This is the central construction of the earlier Stochastic Neighbor Embedding (SNE) method of Hinton and Roweis. 1

It is worth being precise about the word probability. Nothing here says that two observations belong to the same labeled class. The probabilities are assigned by our mathematical neighbor-selection rule, without knowing any class labels. That is why t-SNE can be used on unlabeled data.

Why exponential weights are a natural choice

At this point the exponential might seem pulled out of thin air. There is a deeper route to it. Suppose we know two things about our neighbor selection: the probabilities must add to one, and on average we prefer short squared distances. We do not want to impose any other preference that our assumptions have not justified. A standard way to express that restraint is to choose the probability distribution with the largest entropy subject to a fixed mean squared distance.

Write the unknown probabilities as \(r_j\) and the squared distances from our selected point as \(d_j^2\). We require \(\sum_jr_j=1\) and \(\sum_jr_jd_j^2=R\) for some prescribed average \(R\). We maximize the entropy \(-\sum_jr_j\ln r_j\) while enforcing these two conditions. Introducing multipliers \(a\) and \(b\), we differentiate

\[\mathcal L=-\sum_jr_j\ln r_j-a\Big(\sum_jr_j-1\Big)-b\Big(\sum_jr_jd_j^2-R\Big).\]

Setting the derivative with respect to each \(r_j\) to zero gives \(-\ln r_j-1-a-bd_j^2=0\). Rearranging and normalizing yields

\[r_j=\frac{e^{-bd_j^2}}{\sum_ke^{-bd_k^2}}.\]

So an exponential of the negative squared distance appears when we combine a preference for nearby points with a maximum-entropy assumption. Setting \(b=1/(2\sigma_i^2)\) gives precisely the Gaussian-shaped neighborhood probabilities used in SNE and t-SNE. This is a motivation, not a proof that a Gaussian is the only sensible choice. We deliberately chose squared Euclidean distance and a particular uncertainty principle; changing those assumptions would lead to another similarity kernel. The important point is that the formula now has an intelligible origin.

A numerical example

Suppose that, from one selected point, three other points lie at distances \(1\), \(2\), and \(4\). Set \(\sigma_i=1\). The unnormalized weights are

\[\begin{aligned}w_1&=e^{-1/2}\approx0.6065,\\w_2&=e^{-4/2}\approx0.1353,\\w_3&=e^{-16/2}\approx0.000335.\end{aligned}\]

Their sum is approximately \(0.7422\). Dividing each weight by this sum gives roughly \(81.7\%\), \(18.2\%\), and \(0.05\%\). One observation dominates the neighborhood; the distant candidate barely contributes. Nothing discontinuous happened: we obtained a graded judgment from ordinary distances.

Interactive experiment: A probabilistic neighborhood

Change the Gaussian bandwidth. Every bar shows the probability with which the selected blue point would pick another point as its neighbor. The effective number of neighbors is calculated from the entropy, not set directly.

Perplexity: 1.61
\u0078ᵢ

One Neighborhood Radius Cannot Fit Everyone

We have quietly introduced a new problem. What width \(\sigma_i\) should the Gaussian have? A small value makes our selected point notice only very close neighbors. A large value gives more distant points meaningful probability. But the same numerical width behaves differently in a dense group and in a sparse one.

Imagine two towns. In the center of the first, ten houses fit within a short walk. In the countryside surrounding the second, ten houses may be separated by several kilometers. If we say that every house should pay attention to its local community, a fixed radius makes little sense. We should choose a radius adapted to local density.

That suggests specifying not the radius itself, but roughly how many other points should matter. We need a smooth measure of the number of plausible neighbors—one that still works when the probabilities are unequal.

The starting point is uncertainty. Suppose a random choice selects one neighbor with certainty. There is no uncertainty about the result. If it chooses uniformly among four neighbors, each outcome is possible and the uncertainty is greater. Shannon entropy measures precisely this distinction. For the conditional distribution around \(\mathbf x_i\), we write

\[H(P_i)=-\sum_{j\ne i}p_{j\mid i}\log_2p_{j\mid i}.\]

This formula deserves a derivation of its meaning. A choice with probability \(p\) carries \(-\log_2 p\) bits of surprise. Rare choices are more surprising than common ones. Averaging that surprise over all possible choices gives \(-\sum p\log_2p\), which is Shannon entropy.

Now imagine that exactly \(k\) neighbors are equally likely. Every one has probability \(1/k\). Substituting this simple case gives

\[H=-\sum_{j=1}^{k}\frac1k\log_2\!\frac1k=\log_2k.\]

We want to recover \(k\), so we undo the base-two logarithm. This gives the quantity t-SNE calls perplexity:

\[\boxed{\operatorname{Perp}(P_i)=2^{H(P_i)}.}\]

We did not invent a new neighborhood count by decree: we asked what number of equally probable neighbors would produce the same uncertainty. If four neighbors are equally likely, perplexity is four. If one neighbor has probability \(0.8\) and two others have \(0.1\) each, perplexity is about \(1.89\), reflecting a neighborhood dominated by one candidate.

That number need not be an integer. It is an effective neighbor count, not a hard selection of exactly that many observations.

Dense region: small σSparse region: large σ
Equal perplexity does not imply equal physical radius or a fixed set of exactly \(k\) neighbors.

In practice, the user selects a target perplexity \(K\), and t-SNE varies each \(\sigma_i\) until \(\operatorname{Perp}(P_i)\approx K\). It commonly does this through binary search. If a Gaussian is too narrow, the probability concentrates on very few neighbors and perplexity is low. Broadening it spreads probability across more candidates and raises perplexity. The width is an intermediate quantity determined by the target, not a synonym for perplexity.

There are limits: for a finite dataset of \(n\) observations the effective neighbor count cannot exceed \(n-1\). Larger target perplexities make similarities more global; smaller ones emphasize finer local variations. Neither is automatically more truthful. 2

A Neighbor Relationship Has Two Directions

Here is a subtle consequence of our construction. Even though Euclidean distance is symmetric, the resulting conditional probability generally is not:

\[d_{ij}=d_{ji},\qquad p_{j\mid i}\ne p_{i\mid j}.\]

Why? Because every point has its own width \(\sigma_i\) and its own normalizing denominator. A point in a sparse region might regard a nearby member of a dense group as an unusually important neighbor. That member of the dense group may have many equally close alternatives and assign the first point relatively little probability.

But the map we want to draw has only one distance between any two points. It seems natural to give each pair a single, symmetric similarity. We could average the two conditional probabilities. Yet we also want all pair similarities together to form one probability distribution.

Notice the bookkeeping. For each of the \(n\) points, the conditional probabilities sum to one. Across all selected points they therefore sum to \(n\). If we include both directions, their combined sum is \(2n\). Consequently, averaging the directions and normalizing globally leads to

\[\boxed{p_{ij}=\frac{p_{j\mid i}+p_{i\mid j}}{2n},\qquad p_{ii}=0.}\]

Now \(p_{ij}=p_{ji}\) and \(\sum_{i\ne j}p_{ij}=1\). The factor \(2n\) is not an arbitrary correction. It is what makes the symmetric pair similarities add up to one. We have arrived at the joint probability distribution \(P\) used in symmetric t-SNE. 2

The original data will not move during optimization. Once these \(p_{ij}\) values have been computed, they remain our fixed description of the relationships we are trying to preserve.

Why Copying the Gaussian Does Not Work

We now need a second distribution, describing the relationships in the map. An obvious idea is to place another Gaussian kernel around each map point and choose the new positions so that its probabilities match the original ones. This is closely related to the earlier SNE approach. Why not stop there?

To see the difficulty, forget probabilities for a moment and think about space. A circle of radius \(r\) in two dimensions has area \(\pi r^2\). A sphere in three dimensions has volume proportional to \(r^3\). More generally, the volume of a ball in \(d\) dimensions is proportional to \(r^d\).

Suppose we double the radius of a ball. Its volume increases by a factor of \(2^d\). In two dimensions the larger ball has four times the area; in ten dimensions it has \(1024\) times the volume. The shell between radius \(r\) and \(2r\) therefore contains \(2^d-1\) times as much volume as the original ball: only three times as much in two dimensions, but \(1023\) times as much in ten.

When a ball's radius doubles 2 dimensions ×4 volume 10 dimensions ×1024 Bar length indicates the exponent d in 2^d (not a geometric drawing of a 10D ball).
With twice the radius, volume scales by 2 to the power of the dimension. High-dimensional neighborhoods have much more room around them than a plane can reproduce.

This is not a literal count of dataset points, nor a complete proof of what an embedding will do. It exposes the geometric mismatch: high-dimensional space has a very different balance of nearby and moderately distant locations than a two-dimensional map. When many observations have to occupy the relatively small amount of space surrounding a local neighborhood in the map, they tend to crowd one another.

A Gaussian kernel in the map makes the problem worse because its weights become tiny very quickly with distance. To make moderately dissimilar points have extremely small Gaussian affinities, we may need to push them far away. With many pairs competing for room, the resulting arrangement becomes difficult to optimize. Van der Maaten and Hinton identified this as the crowding problem. 2

We need a similarity rule for the map that still gives high weights to close neighbors, but does not make the weight of a moderately distant pair vanish almost immediately.

Why the Letter t Appears

Let us specify what would help. We want a positive, smoothly decreasing function of distance. It should be strong at short distances, so real neighbors can pull together. But it should have a heavier tail than a Gaussian: even fairly separated points should retain some non-negligible weight.

A particularly simple choice is

\[v(r)=\frac{1}{1+r^2}.\]

At zero distance this is one. At distance one it is one half. At distance three it is one tenth. For large \(r\), it falls approximately as \(1/r^2\). Compare that with a Gaussian of unit scale, which falls as \(e^{-r^2/2}\). At distance three the Gaussian weight is about \(0.011\), while the new function still gives \(0.1\). The two functions tell us substantially different stories about distant points.

0101234567 Gaussian exp(−r²/2)Student-t kernel 1/(1+r²)
Both curves are rescaled to 1 at zero distance, so the comparison is between affinity functions, not normalized probability densities. The blue tail decays much more slowly.

This particular rational function is proportional to the density of a Student-t distribution with one degree of freedom, which in one dimension is also called the Cauchy distribution:

\[f(t)=\frac{1}{\pi(1+t^2)}.\]

The constant \(1/\pi\) ensures that the one-dimensional density integrates to one. For t-SNE we do not need this constant; we use the shape as a pairwise affinity kernel and normalize the resulting pair weights later. In particular, the kernel \((1+\|\mathbf y_i-\mathbf y_j\|^2)^{-1}\) should not be confused with the full normalized density of a multivariate Student-t distribution.

We can now define the similarity of a pair in our two- or three-dimensional map:

\[v_{ij}=\frac{1}{1+\|\mathbf y_i-\mathbf y_j\|^2}.\]

The same logic as before applies. If we want a probability distribution, divide by the sum of all eligible weights. Because we now work with symmetric pairs rather than separately normalized neighborhoods, the denominator covers all ordered pairs of different points:

\[\boxed{q_{ij}=\frac{v_{ij}}{Z},\qquad Z=\sum_{k\ne l}v_{kl}=\sum_{k\ne l}\frac{1}{1+\|\mathbf y_k-\mathbf y_l\|^2}.}\]

Thus \(q_{ii}=0\), \(q_{ij}=q_{ji}\), and \(\sum_{i\ne j}q_{ij}=1\). We finally have the two objects we wanted: \(P=(p_{ij})\) records the neighbor relationships of the original observations, while \(Q=(q_{ij})\) records the neighbor relationships of the current drawing.

The heavy tail gives the map more freedom to place unrelated regions apart while still distinguishing neighbors. It does not guarantee the recovery of exact global distances or the original topology. The map remains a compromise.

Interactive experiment: What the tail changes

Move two map points farther apart. The two values below are unnormalized weights for the same distance, so their different rates of decay are directly comparable.

AB
Gaussian: 0.1353Student-t: 0.2000

What Should Count as an Error?

The simplest way to compare \(P\) and \(Q\) might be to subtract their entries. But this would hide an important asymmetry: losing a genuine neighbor is not the same mistake as accidentally placing two unrelated observations rather close together.

We can derive an alternative by asking how surprising a map would be when the original data generates a random pair \((i,j)\) according to probabilities \(p_{ij}\). Our low-dimensional map assigns that pair probability \(q_{ij}\). If the map gives a tiny probability to a pair that the original data chooses frequently, it makes a poor prediction.

Information theory measures the surprise of an event assigned probability \(q\) by \(-\log q\). A common event assigned very low probability is surprising; a common event assigned high probability is not. Since the original distribution \(P\) determines how often events occur, the average surprise under the map is

\[H(P,Q)=-\sum_{i\ne j}p_{ij}\log q_{ij}.\]

This is the cross-entropy. Why don't we simply minimize it? In fact we can, because the entropy of the original distribution,

\[H(P)=-\sum_{i\ne j}p_{ij}\log p_{ij},\]

is fixed. The original similarities do not change when we move points on the map. Subtracting this fixed baseline from cross-entropy gives

\[H(P,Q)-H(P)=\sum_{i\ne j}p_{ij}\log\frac{p_{ij}}{q_{ij}}.\]

The result is the Kullback–Leibler divergence, and we have now derived why it is the natural cost function here:

\[\boxed{C(\mathcal Y)=D_{\mathrm{KL}}(P\|Q)=\sum_{i\ne j}p_{ij}\log\frac{p_{ij}}{q_{ij}}.}\]

The divergence is nonnegative and equals zero only when the two probability distributions agree. It is not a distance in the usual geometric sense: generally \(D_{\mathrm{KL}}(P\|Q)\ne D_{\mathrm{KL}}(Q\|P)\). Individual summands can even be negative; the complete sum cannot. The logarithm in this formula is the natural logarithm. In the perplexity definition we used base two to express entropy in bits; these conventions are compatible.

The direction \(P\|Q\) is deliberate. When an original pair has appreciable \(p_{ij}\) but the map gives it a tiny \(q_{ij}\), the logarithmic mismatch receives substantial weight. When \(p_{ij}\) is tiny, even a considerable relative disagreement contributes little by itself. Of course, the \(q_{ij}\) must still sum to one, so moving points can change many terms at once. The overall effect is a preference for preserving strong local relationships rather than every long-range separation. 2

The original problem has become an optimization problem. We no longer have to invent the coordinates of the map directly. We need to find the positions \(\mathbf y_i\) that make this cost function small.

If the Map Is Wrong, Which Way Should a Point Move?

This is the point where an equation often appears without explanation: the gradient of the t-SNE objective. But it too follows from a few simple observations.

First, only \(Q\) depends on the map positions. The values in \(P\) are fixed. Recall that \(q_{ij}=v_{ij}/Z\). Taking the logarithm gives \(\log q_{ij}=\log v_{ij}-\log Z\). Inserting this into the cost function, and using \(\sum_{i\ne j}p_{ij}=1\), produces

\[C(\mathcal Y)=\underbrace{\sum_{i\ne j}p_{ij}\log p_{ij}}_{\text{constant}}-\sum_{i\ne j}p_{ij}\log v_{ij}+\log Z.\]

This rearrangement is useful because it separates the objective into two parts that exert opposing influences. The term involving \(-\log v_{ij}\) penalizes pairs that should be close but are far apart. The \(\log Z\) term arises only because all map affinities must sum to one; it prevents every point from being moved together without considering the rest.

Now let \(r_{ij}^2=\|\mathbf y_i-\mathbf y_j\|^2\). Because the squared distance is a dot product, its derivative with respect to \(\mathbf y_i\) is

\[\frac{\partial r_{ij}^2}{\partial\mathbf y_i}=2(\mathbf y_i-\mathbf y_j).\]

We also know that \(v_{ij}=1/(1+r_{ij}^2)\), so \(-\log v_{ij}=\log(1+r_{ij}^2)\). By the chain rule,

\[\frac{\partial[-\log v_{ij}]}{\partial\mathbf y_i}=\frac{2(\mathbf y_i-\mathbf y_j)}{1+r_{ij}^2}=2v_{ij}(\mathbf y_i-\mathbf y_j).\]

There is a factor of two in this derivative. A second factor of two appears because each symmetric pair is included twice in our sum, once as \((i,j)\) and once as \((j,i)\). Consequently, the first variable term contributes

\[4\sum_{j\ne i}p_{ij}v_{ij}(\mathbf y_i-\mathbf y_j).\]

Next comes the denominator \(Z\). Its presence might look like an unpleasant technical detail, but it is the source of the repulsion that makes the algorithm work. Differentiating \(Z=\sum_{k\ne l}v_{kl}\) gives

\[\frac{\partial Z}{\partial\mathbf y_i}=-4\sum_{j\ne i}v_{ij}^2(\mathbf y_i-\mathbf y_j).\]

And because the derivative of \(\log Z\) is \((1/Z)\,\partial Z\), while \(v_{ij}/Z=q_{ij}\), we obtain

\[\frac{\partial\log Z}{\partial\mathbf y_i}=-4\sum_{j\ne i}q_{ij}v_{ij}(\mathbf y_i-\mathbf y_j).\]

Put the two contributions together and the pieces cancel into a compact expression:

\[\boxed{\frac{\partial C}{\partial\mathbf y_i}=4\sum_{j\ne i}(p_{ij}-q_{ij})\frac{\mathbf y_i-\mathbf y_j}{1+\|\mathbf y_i-\mathbf y_j\|^2}.}\]

Every symbol now has a reason to be there. The vector \(\mathbf y_i-\mathbf y_j\) points along the line between the map points. The difference \(p_{ij}-q_{ij}\) describes whether their relationship is underrepresented or overrepresented. The denominator \(1+\|\mathbf y_i-\mathbf y_j\|^2\) is the Student-t kernel we chose to solve the crowding problem. The leading four comes from differentiating a squared distance and counting two directions of a symmetric pair.

Remember that a gradient points uphill. To reduce the error, we must move against it. If \(p_{ij}>q_{ij}\), a pair is less similar on the map than in the original data; descending the gradient pulls the points toward one another. If \(p_{ij}<q_{ij}\), the pair is more similar on the map than in the original data, so gradient descent pushes the points apart. This repulsion is not a physical rule appended to t-SNE; it is what the probability mismatch produces when differentiated.

p > q : attractionp < q : repulsion Too little similarity in the mapToo much similarity in the map
The arrows indicate the direction of the coordinate update, which is the negative gradient, not the gradient itself.

From One Small Movement to an Entire Map

We begin with some initial positions \(\mathbf y_i^{(0)}\), often very close to the origin or obtained from an informative initialization such as PCA. We compute the gradient at these positions, move slightly downhill, and repeat.

With a learning rate \(\eta>0\), the basic update is

\[\mathbf y_i^{(t+1)}=\mathbf y_i^{(t)}-\eta\,\nabla_{\mathbf y_i}C(\mathcal Y^{(t)}).\]

Too small a learning rate can make progress slow; too large a rate can make the optimization unstable. The objective is nonlinear and nonconvex, so there can be different useful solutions depending on initialization and settings. No theorem promises that one run will discover a unique best-looking map.

Many implementations add momentum. If \(\boldsymbol\Delta_i^{(t)}\) is the previous displacement of a point, the next can be chosen as

\[\begin{aligned}\boldsymbol\Delta_i^{(t+1)}&=\mu\boldsymbol\Delta_i^{(t)}-\eta\,\nabla_{\mathbf y_i}C,\\\mathbf y_i^{(t+1)}&=\mathbf y_i^{(t)}+\boldsymbol\Delta_i^{(t+1)}.\end{aligned}\]

The coefficient \(\mu\) retains some motion from the preceding step. The intuition is the same as pushing a heavy object down an uneven slope: persistence can reduce unproductive zigzagging. The original t-SNE publication used momentum and adaptive learning rates as practical optimization devices. 2

Nothing in this derivation required \(\mathbf y_i\) to have exactly two coordinates. The gradient is a vector in whatever output dimension we choose. In 2D it has two components; in 3D it has three. The objective, pair similarities, and optimization principle are otherwise the same. The third coordinate can provide room for a configuration that cannot be drawn without overlap in a plane, although it also makes visual inspection dependent on viewpoint.

Why the First Steps Are Treated Differently

An additional problem arises at the beginning of optimization. When all map points start in a tiny region, the final separation of meaningful groups has not yet formed. Nearby observations should first recognize and gather around one another, but repulsion is already acting on the entire map.

Early exaggeration temporarily makes the attractive part stronger. We replace \(p_{ij}\) in the force expression by \(\alpha p_{ij}\), where \(\alpha>1\), for an initial phase:

\[\mathbf g_i^{(\alpha)}=4\sum_{j\ne i}(\alpha p_{ij}-q_{ij})\frac{\mathbf y_i-\mathbf y_j}{1+\|\mathbf y_i-\mathbf y_j\|^2}.\]

This change lets strongly connected observations contract into compact neighborhoods before we restore the normal balance of attraction and repulsion. It is helpful to think of the first phase as roughly establishing local groups, and the second as allowing those groups to expand and arrange themselves. That is an intuition about the optimization, not a promise of perfect clustering.

There is an important technical distinction. Since \(\sum p_{ij}=1\), multiplying every \(p_{ij}\) by \(\alpha\) makes the sum equal \(\alpha\), not one. During early exaggeration, these scaled numbers are no longer a probability distribution. We should therefore treat the expression as a modified gradient or force rule, not casually call it the gradient of the same normalized KL objective.

The 2008 paper used a factor of four during an initial phase of 50 iterations in its reported experiments; the 2014 tree-based implementation described a factor of twelve for 250 iterations. Those values came from particular implementations and experiments. They are not constants that follow from the definition of t-SNE. 2, 3

Later, Linderman and Steinerberger studied the exaggeration phase mathematically. Under explicit assumptions about well-separated input clusters and parameter choices, they proved rapid contraction of embedded clusters and identified a relationship to spectral clustering. The relevant regime includes a parameter order \(\alpha\sim n/10\) together with a suitable step size, not the inverse \(10/n\). This is a conditional theoretical result, not a universal prescription for all datasets or implementations. 4

Interactive experiment: Let the points find their places

The same 36 synthetic observations have four input coordinates. Choose a 2D or 3D output, then start the exact pairwise t-SNE update. In 3D, rotate the projected view; the underlying optimization really has three coordinates. Colors only reveal which input groups were generated; they are not given to the optimizer.

Iteration: 0
36 points · exact pairwise t-SNE
The map begins with overlapping positions. Run the optimizer to see nearby observations attract and non-neighbors repel.

What Two Dimensions Show, and What Three Can Add

The usual t-SNE plot has two axes because it is immediately readable. You can inspect a two-dimensional scatter plot without manipulating it, print it, or place it next to other figures. This simplicity is a substantial advantage. Yet the fact that we can draw a point cloud in two dimensions does not mean its structure is inherently two-dimensional.

A three-dimensional embedding assigns each observation a third coordinate and optimizes all three simultaneously. This can sometimes accommodate relationships that crowd together in 2D. However, it also introduces an observational problem: when a 3D cloud is rendered on a flat screen, two distant points may appear to overlap because one lies behind the other. Rotation, depth cues, and an interactive viewpoint are therefore essential.

Neither choice changes the conceptual objective: preserve pairwise neighbor probabilities as well as the available space allows. The extra dimension offers more room; it does not restore the full original geometry. If the original data has dozens or hundreds of meaningful dimensions, even 3D is a drastic simplification.

The interactive optimizer above allows the same synthetic four-dimensional input observations to be embedded in either 2D or 3D. The 3D visualization can be rotated, but the underlying map is actually optimized in three coordinates; it is not merely a 2D result drawn with perspective.

A Real Dataset in the Browser: MNIST

The handwritten digits of MNIST provide a useful counterpoint to our imaginary monsters. Here no learned representation is required to begin: each 28-by-28 image already gives us exactly 784 numerical pixel values. The dataset has 60,000 training images and 10,000 test images. We can treat the pixel intensities themselves as features, calculate pairwise distances, and ask t-SNE to preserve the resulting neighborhoods. Digit labels can be shown as colors afterward, but they must not enter the similarity calculation if we want an unsupervised experiment. [10]

Why not calculate all 70,000 images in a beginner's browser demo? The raw pixel bytes are actually modest—70,000 times 784 equals 54,880,000 bytes, before compression—but the exact t-SNE optimizer compares pairs. For \(n\) points, this introduces approximately \(n^2\) interactions. Seventy thousand points produce almost five billion ordered pairs. The difficulty is not primarily downloading the digits; it is optimizing their arrangement repeatedly.

For teaching, a sample of 200–700 real MNIST digits is more practical. Their actual neighbor probabilities can be computed and the map watched during gradient descent. Showing the original image on hover would make relationships inspectable, while digit labels should remain outside the similarity calculation in an unsupervised experiment. The optimizer in this chapter instead uses a small synthetic dataset so its exact pairwise computation stays responsive.

Preprocessing also changes the input geometry. Comparing original 28-by-28 pixel features with averaged 14-by-14 features, for example, would test how much neighborhood structure survives that representation change. Such choices are part of the experiment, not merely implementation details.

A different design is appropriate when the goal is simply to browse tens of thousands of digits. Then the best experience is often to calculate a large embedding beforehand and send only two or three coordinates per sample, plus labels and thumbnail images, to the browser. A 70,000-point 2D map needs only 560,000 bytes for two 32-bit floating-point coordinates per point, before compression. Rendering a prepared result is cheap; optimizing it is the expensive part. The distinction between these two experiences—watching an algorithm think and exploring its finished map—is worth keeping explicit.

The distinction is useful when choosing an experiment: a prepared large map is inexpensive to browse, while recomputing an exact embedding is the costly operation.

Beyond a Static Map: Where Does a New Point Go?

Imagine that our MNIST map already contains 10,000 handwritten digits. A new image arrives. The new image is again a vector of 784 pixel intensities, but classical t-SNE has not learned a function that takes this vector and returns a location. It has learned a finite list of positions for the original observations. Formally, the result is \(\{\mathbf y_1,\ldots,\mathbf y_n\}\), not a general rule \(f:\mathbb R^d\to\mathbb R^s\). This is why an ordinary t-SNE plot is not automatically a streaming visualization.

One solution is to start over with all \(n+1\) observations. That changes the neighborhood probabilities and requires another optimization; it may also rearrange existing points. If we were using the map as an interface, the entire landscape could shift whenever one observation arrived. Even if the new result were mathematically reasonable, it would make it harder for a human observer to maintain a sense of location.

A simpler alternative is to freeze the old map and solve only for the location of the newcomer. First, calculate its input-space distance to every existing observation. Use the same idea as before to obtain a probability distribution over candidate neighbors:

\[p_{i\mid *}=\frac{\exp(-\|\mathbf x_* -\mathbf x_i\|^2/(2\sigma_*^2))}{\sum_{k=1}^{n}\exp(-\|\mathbf x_* -\mathbf x_k\|^2/(2\sigma_*^2))}.\]

The asterisk marks the new observation. The bandwidth \(\sigma_*\) can be selected by a target perplexity, just as for any other point. Next, propose a position \(\mathbf y_*\) on the old map. The browser demo stores these positions in SVG screen units, so it uses a scale factor \(a^2=3900\) in the output-space kernel. This is a coordinate-unit choice. Its similarities are normalized only over candidates for this new observation:

\[v_i(\mathbf y_*)=\left(1+\frac{\|\mathbf y_*-\mathbf y_i\|^2}{a^2}\right)^{-1},\qquad q_{i\mid *}(\mathbf y_*)=\frac{v_i(\mathbf y_*)}{\sum_{k=1}^{n}v_k(\mathbf y_*)}.\]

We want a position at which the map's proposed neighbor probabilities resemble the input probabilities. For this small, conditional out-of-sample problem we can minimize

\[C_*(\mathbf y_*)=\sum_{i=1}^{n}p_{i\mid *}\log\frac{p_{i\mid *}}{q_{i\mid *}(\mathbf y_*)}.\]

Notice the difference from classical, joint t-SNE: the denominator covers only the candidate neighbors of the new point, and the existing coordinates are held fixed. We have derived a useful placement procedure, not a claim that we have exactly updated the original global t-SNE solution. Reoptimizing all positions and pair probabilities would generally produce another answer.

Even this smaller problem has a gradient. Differentiating the scaled kernel and its normalizing sum gives

\[\nabla_{\mathbf y_*}C_* = \frac{2}{a^2}\sum_{i=1}^{n}(p_{i\mid *}-q_{i\mid *})\,v_i(\mathbf y_*)\,(\mathbf y_*-\mathbf y_i).\]

The factor is two, not four, because the conditional objective does not count each symmetric pair twice. We can optimize only \(\mathbf y_*\) while all \(\mathbf y_i\) remain fixed. This makes a practical, explainable insertion demo possible, even though it is not equivalent to retraining t-SNE.

The map at time tA new observation at time t + 1 Old coordinates already existOnly the new green point is placed
A new observation can be fitted into an existing map without moving the old points.

An out-of-sample insertion keeps the previous layout stationary. It is a useful approximation, not an exact global re-embedding.

Interactive experiment — place a new observation. Choose a family of synthetic observations, insert a new point, and watch its conditional neighbor distribution determine where it moves. The existing map stays fixed. This experiment is a schematic frozen-map insertion rather than a full incremental t-SNE implementation.

Interactive experiment — where does a new point go?

Each family is a toy set of observations with five input features. The older points have fixed locations on a schematic map. Insert a new observation; only its coordinates are optimized to match conditional neighbor probabilities. This is not a full t-SNE recomputation.

8
Select a family and click Insert new point.

There is another possibility: instead of storing only a list of coordinates, we can learn a mapping \(f_\theta(\mathbf x)\). A neural network, for example, can be trained so that the neighbor probabilities produced by its outputs have low t-SNE-style divergence. Once trained, a fresh observation gets an output position simply by evaluating the network. Van der Maaten introduced this approach as parametric t-SNE in 2009. The additional machinery buys something the original embedding does not provide: an explicit, reusable function. 12

We should distinguish speed from faithfulness. A fixed learned function makes new points cheap to project, but it need not reproduce the result of rerunning the nonparametric optimization when the underlying data distribution changes. In a live system, the representation, neighborhood index, and projection model may all need periodic revision.

When the Data Itself Moves

Inserting a previously unseen observation is not the same problem as following the same observation through time. Imagine a sensor attached to a motor. At each moment it reports a high-dimensional spectrum. We want to understand whether its operating state drifts gradually or changes abruptly. If we calculate a completely independent t-SNE map at every moment, points may rotate, flip, expand, or rearrange for reasons that have nothing to do with changes in the motor. After all, translating or rotating an entire t-SNE map leaves all its pairwise distances unchanged. The objective itself does not give the map an absolute orientation.

This creates two competing wishes. At each time \(t\), we want to represent the current data faithfully. At the same time, a reading that has barely changed should not jump across the screen. An intuitive way to express both demands is to start with the ordinary t-SNE objective at each time and penalize movement between successive maps:

\[C_{\mathrm{time}}=\sum_{t=1}^{T}D_{\mathrm{KL}}(P^{(t)}\|Q^{(t)})+\lambda\sum_{t=2}^{T}\sum_i\|\mathbf y_i^{(t)}-\mathbf y_i^{(t-1)}\|^2.\]

This is a simplified illustrative temporal regularization objective in the spirit of dynamic t-SNE, not a statement that every dynamic t-SNE implementation uses this exact expression. The first term asks whether the individual maps represent their current neighborhoods. The second says that large unexplained jumps should cost something. The nonnegative parameter \(\lambda\) controls the compromise. With \(\lambda=0\), the frames are free to differ. If \(\lambda\) is excessively large, a genuinely changing structure may appear almost frozen. Rauber, Falcão, and Telea investigated this trade-off in their original work on dynamic t-SNE. 13

There are necessary assumptions. We must know which observation at time \(t\) corresponds to which observation at time \(t+1\), or have an explicit rule for matching them. If observations enter and leave the system, a simple one-to-one temporal penalty is not enough. And if our goal is to study motion, we must take care not to mistake differences of alignment or optimization convergence for physical change.

This distinction has practical consequences. A static map is useful for exploring a finite collection. Out-of-sample placement adds observations without destabilizing the old layout. A parametric map produces reusable coordinates for new inputs. A dynamic map explicitly represents time and balances visual stability against current geometric fidelity. The four are related, but they solve different problems. Thinking of them as one algorithm would obscure the very design choices we are trying to understand.

What We Can—and Cannot—Learn From the Picture

Once the algorithm finishes, we may see islands, bands, arcs, or loose clouds. These shapes are useful clues, but they must not be interpreted more literally than the optimization warrants.

First, local neighborhoods are the strongest evidence. Two observations that remain close over several reasonable parameter settings are candidates for genuine similarity in the original feature space. To establish that they are in fact similar, inspect their original vectors or reconstruct the nearest-neighbor relationships there. A dramatic separation between two distant clusters is much less directly meaningful.

Second, the size of an island is not the original variance of its points. Because every point uses a locally adapted input bandwidth and the map includes global repulsion, dense and sparse regions can be expanded or compressed differently. A wide cluster on the plot is not necessarily a more variable population; a compact cluster is not necessarily a denser one.

Third, islands are not automatically natural categories. A continuous structure can be broken into apparent groups by a particular set of parameters; true subgroups can also be merged. The number of colored clusters in an image is not something that t-SNE itself proves. It is especially easy to be misled when labels or colors are added only after looking at a compelling plot.

Fourth, global directions are arbitrary. Translating, rotating, or reflecting the whole embedding does not alter its pairwise distances or objective value. Unlike PCA components, the first and second t-SNE coordinates are not named underlying features. The same applies to the orientation of a 3D embedding.

Finally, one plot is not an experiment. Different random initializations, perplexities, feature preprocessing choices, and optimization settings can lead to different pictures. Research on single-cell datasets has shown how strongly global arrangement can depend on parameter choices and initialization. A scientific interpretation should therefore compare plausible settings and validate findings in the original data. 5, 6

A Sensible Procedure Before Trusting a Map

Start with the original measurements. Decide what the features mean, whether they should be scaled, and whether Euclidean distance is a reasonable proxy for meaningful similarity. There is no mathematical miracle at the end of the pipeline that repairs a poorly chosen notion of distance at the beginning.

For very high-dimensional or noisy observations, a moderate preliminary PCA reduction can be useful. This is not the same as doing PCA instead of t-SNE. PCA can first remove some noise and reduce computational burden; t-SNE then attempts to preserve neighborhoods in the resulting representation. The suitability of this sequence depends on the data.

Next choose the visualization dimension: 2D for compact plots, or 3D when an interactive view and the additional space are valuable. Try several reasonable perplexities and repeat the optimization with different seeds or initializations. Record the preprocessing, metric, output dimension, perplexity, learning rate, exaggeration schedule, and iteration count so that the result can be reproduced.

In the simplest implementation, every iteration must consider interactions among approximately \(n(n-1)\) ordered pairs, giving quadratic cost. Barnes–Hut t-SNE uses a tree-based approximation for distant repulsive interactions, reducing much of the gradient calculation to about \(O(n\log n)\). Later interpolation-based methods further improve large-scale computation. These accelerations make the algorithm useful on datasets that would be expensive to process with the exact pairwise force calculation. 3, 7

Finally, ask a specific scientific question of the plot. Are nearest neighbors of a particular observation preserved? Do known but unused labels agree with emergent neighborhoods? Does a transition appear across multiple settings? Can a proposed group be distinguished in the original measurements? These questions convert a picture into evidence.

The Whole Construction, in One Thought

We began with high-dimensional points and admitted that a two- or three-dimensional picture could not preserve their entire geometry. We therefore decided to prioritize who is near whom. A smoothly decreasing Gaussian turned distances into preferences, and normalization made those preferences probabilities. Entropy gave us a way to specify how many effective neighbors each point should consider; symmetrization gave us one distribution over pairs.

The crowding problem then explained why we should not use an equally short-tailed Gaussian on the map. A heavy-tailed Student-t kernel produced a second pair distribution that permits more separation. We derived the KL divergence by asking how surprising the original neighbor pairs would be under the map, then differentiated that objective and discovered forces that pull underrepresented neighbors together and push overrepresented pairs apart.

This is the central idea of t-SNE: do not ask a low-dimensional map to preserve everything. First decide what constitutes a meaningful relationship, express that relationship as a probability, and let optimization construct a picture that preserves those probabilities as closely as possible. The result is not an exact miniature of the original space. It is a carefully negotiated view of its neighborhoods. And if we later want to insert fresh observations or follow a changing process, we must make another choice: should the old map remain fixed, should we learn a reusable projection function, or should the map itself evolve through time? These are not implementation footnotes. They are new mathematical objectives, arising from new questions.

References