t-SNE (t-Distributed Stochastic Neighbor Embedding) › Building the Mathematical Foundation

To understand how t-SNE works, we need to follow its core logic: the algorithm tries to preserve relationships between nearby points when moving from high dimensions to low dimensions. Instead of trying to preserve exact distances, t-SNE preserves similarities. Think of it as keeping friends close together rather than measuring exact distances between people in a room.

The algorithm operates in three main steps. First, it measures how similar points are in the original high-dimensional space. Second, it creates a low-dimensional representation and measures similarities there. Third, it adjusts the low-dimensional positions until the two sets of similarities match as closely as possible.

We start with our original data points $x_i$ and $x_j$ in high-dimensional space. The first step is straightforward: compute how far apart these points are using squared Euclidean distance.

$$d_{ij} = ||x_i - x_j||^2$$

where:

• $x_i$: the $i$-th data point in the original high-dimensional space
• $x_j$: the $j$-th data point in the original high-dimensional space
• $d_{ij}$: squared Euclidean distance between points $i$ and $j$
• $||\cdot||^2$: squared L2 norm (sum of squared differences across all dimensions)

This gives us raw distances, but t-SNE doesn't work directly with distances. Instead, it converts them into probabilities that express similarity. Why probabilities? Because they naturally normalize the scale and make similarities relative rather than absolute. A distance of 5 units means something different if most points are separated by 100 units versus if most are separated by 1 unit.

To convert distances into similarities, we use a Gaussian (normal) distribution centered at each point. The conditional probability $p_{j|i}$ represents the probability that point $x_i$ would pick $x_j$ as its neighbor if neighbors were selected in proportion to their probability density under this Gaussian.

$$p_{j|i} = \frac{\exp(-d_{ij}/2\sigma_i^2)}{\sum_{k \neq i} \exp(-d_{ik}/2\sigma_i^2)}$$

where:

• $p_{j|i}$: conditional probability that point $i$ picks point $j$ as its neighbor
• $\sigma_i$: bandwidth of the Gaussian kernel centered at point $i$ (controls the effective neighborhood size)
• $\exp(\cdot)$: exponential function (converts distances to similarities)
• $\sum_{k \neq i}$: sum over all points $k$ except point $i$ (normalization)

Let's break down what each part does:

• The numerator $\exp(-d_{ij}/2\sigma_i^2)$ gives higher values to closer points and lower values to distant points. The exponential ensures all values are positive.
• The variance $\sigma_i^2$ controls the width of the Gaussian. Larger $\sigma_i$ means we consider more distant neighbors.
• The denominator normalizes these values so they sum to 1, making them proper probabilities.

The variance $\sigma_i$ is different for each point, which is unusual but important. In dense regions of the data, we use smaller $\sigma_i$ to focus on very close neighbors. In sparse regions, we use larger $\sigma_i$ to capture the local structure. The perplexity parameter controls these variances indirectly. Higher perplexity means the Gaussian is wider, effectively considering more neighbors.

Now we have conditional probabilities $p_{j|i}$, but there's a problem: $p_{j|i}$ might not equal $p_{i|j}$. To fix this asymmetry, we symmetrize by averaging:

$$p_{ij} = \frac{p_{j|i} + p_{i|j}}{2n}$$

where:

• $p_{ij}$: symmetric joint probability of similarity between points $i$ and $j$
• $n$: total number of data points in the dataset
• The factor $2n$ ensures that $\sum_{i}\sum_{j} p_{ij} = 1$

This gives us a symmetric joint probability distribution where $p_{ij} = p_{ji}$. The division by $2n$ ensures all probabilities sum to 1 over all pairs. These symmetrized probabilities $p_{ij}$ represent our target: they encode which points should be close together in the visualization.

Next, we need to measure similarities in the low-dimensional space where we want to visualize our data. Let $y_i$ and $y_j$ be the low-dimensional coordinates we're trying to find. We could use the same Gaussian distribution here, but that creates a problem called crowding.

The crowding problem happens because high-dimensional spaces have much more "room" than low-dimensional spaces. Imagine 100 points evenly distributed in 50 dimensions. When we try to place them in 2 dimensions while keeping all pairwise distances the same, they'll be forced too close together in the center. Moderate distances in high dimensions become impossible to represent in low dimensions.

To solve this, t-SNE uses a heavy-tailed distribution in low dimensions: the Student t-distribution with one degree of freedom (also called the Cauchy distribution). The similarity between points $y_i$ and $y_j$ is:

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

where:

• $q_{ij}$: similarity between points $i$ and $j$ in the low-dimensional space
• $y_i$: low-dimensional representation of point $i$ (coordinates we're optimizing)
• $y_j$: low-dimensional representation of point $j$
• $||y_i - y_j||^2$: squared Euclidean distance between points in low-dimensional space
• $(1 + ||y_i - y_j||^2)^{-1}$: Student t-distribution kernel with one degree of freedom
• $\sum_{k} \sum_{l \neq k}$: sum over all distinct pairs of points for normalization

The key difference from the Gaussian is $(1 + ||y_i - y_j||^2)^{-1}$ instead of $\exp(-||y_i - y_j||^2)$. The t-distribution has heavier tails, meaning it decreases more slowly as distance increases. This allows moderately distant points in low dimensions to still have reasonable similarity values, giving them room to spread out instead of crowding in the center.

The denominator again normalizes to make these proper probabilities that sum to 1.

Now we have two probability distributions: $p_{ij}$ representing similarities we want to preserve from the high-dimensional space, and $q_{ij}$ representing similarities in our current low-dimensional representation. We need to adjust the positions $y_i$ until these distributions match.

To measure how well they match, t-SNE uses the Kullback-Leibler (KL) divergence, a standard measure of how different two probability distributions are:

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

where:

• $C$: Kullback-Leibler divergence (cost function to minimize)
• $p_{ij}$: target similarity between points $i$ and $j$ in high-dimensional space
• $q_{ij}$: current similarity between points $i$ and $j$ in low-dimensional space
• $\log$: natural logarithm
• The sum is over all pairs of points

The KL divergence has useful properties for this task. When $p_{ij}$ is large (points should be close) but $q_{ij}$ is small (they're far apart), the log ratio is large and positive, heavily penalizing the mismatch. When $p_{ij}$ and $q_{ij}$ are both small (points should be far apart and are), the contribution is minimal. The cost function is asymmetric: it cares more about keeping close points together than about separating distant points, which aligns with t-SNE's goal of preserving local structure.

To minimize this cost, we use gradient descent. We need the derivative of the cost with respect to each low-dimensional position $y_i$:

$$\frac{\partial C}{\partial y_i} = 4 \sum_j (p_{ij} - q_{ij})(y_i - y_j)(1 + ||y_i - y_j||^2)^{-1}$$

where:

• $\frac{\partial C}{\partial y_i}$: gradient of the cost function with respect to position $y_i$
• The sum is over all other points $j$
• $(p_{ij} - q_{ij})$: difference between target and current similarities
• $(y_i - y_j)$: vector pointing from $y_j$ to $y_i$ in low-dimensional space
• $(1 + ||y_i - y_j||^2)^{-1}$: weighting factor from the t-distribution
• The factor of 4 comes from the derivative of the symmetrized KL divergence

This gradient has an elegant interpretation. For each pair of points:

• If $p_{ij} > q_{ij}$ (points should be closer than they are), the term $(p_{ij} - q_{ij})$ is positive. The gradient pulls $y_i$ toward $y_j$.
• If $p_{ij} < q_{ij}$ (points should be farther apart), the term is negative. The gradient pushes $y_i$ away from $y_j$.
• The magnitude of the force depends on how different the probabilities are.
• The factor $(1 + ||y_i - y_j||^2)^{-1}$ comes from the t-distribution and ensures the gradient doesn't explode for distant points.

The algorithm iteratively updates each $y_i$ by taking small steps in the direction opposite to this gradient, gradually moving points until the low-dimensional similarities $q_{ij}$ match the high-dimensional similarities $p_{ij}$ as closely as possible.

Understanding a few mathematical properties helps explain t-SNE's behavior in practice.

Non-convexity and local minima: The cost function has multiple local minima because it's non-convex. Different random initializations can lead to different final configurations. This is why running t-SNE multiple times with different random seeds is important to ensure your visualization is stable.

Scale invariance: Multiplying all distances by a constant factor doesn't change the relative similarities. This means t-SNE results don't depend on the absolute scale of your data, only on the relative distances. However, the method is not rotation-invariant, so the orientation of your final visualization is arbitrary.

The role of perplexity: The perplexity parameter is a smoothness control. Lower values (5-15) create more localized neighborhoods. Higher values (30-50) create broader neighborhoods. The choice of perplexity significantly affects the result. A common default is 30, which typically balances local and global structure well for most datasets. The perplexity should generally be between 5 and 50, and you should try several values to see which gives the most interpretable visualization.