t-SNE-tutorial

An illustrated introduction to the t-SNE algorithm

In the Big Data era, data is not only becoming bigger and bigger; it is also becoming more and more complex. This translates into a spectacular increase of the dimensionality of the data. For example, the dimensionality of a set of images is the number of pixels in any image, which ranges from thousands to millions.

Computers have no problem processing that many dimensions. However, we humans are limited to three dimensions. Computers still need us (thankfully), so we often need ways to effectively visualize high-dimensional data before handing it over to the computer.

How can we possibly reduce the dimensionality of a dataset from an arbitrary number to two or three, which is what we're doing when we visualize data on a screen?

The answer lies in the observation that many real-world datasets have a low intrinsic dimensionality, even though they're embedded in a high-dimensional space. Imagine that you're shooting a panoramic landscape with your camera, while rotating around yourself. We can consider every picture as a point in a 16,000,000-dimensional space (assuming a 16 megapixels camera). Yet, the set of pictures approximately lie in a three-dimensional space (yaw, pitch, roll). This low-dimensional space is embedded within the high-dimensional space in a complex, nonlinear way. Hidden in the data, this structure can only be recovered via specific mathematical methods.

This is the topic of manifold learning, also called nonlinear dimensionality reduction, a branch of machine learning (more specifically, unsupervised learning). It is still an active area of research today to develop algorithms that can automatically recover a hidden structure in a high-dimensional dataset.

This post is an introduction to a popular dimensonality reduction algorithm: t-distributed stochastic neighbor embedding (t-SNE). Developed by Laurens van der Maaten and Geoffrey Hinton (see the original paper here), this algorithm has been successfully applied to many real-world datasets. Here, we'll follow the original paper and describe the key mathematical concepts of the method, when applied to a toy dataset (handwritten digits). We'll use Python and the scikit-learn library.

Visualizing handwritten digits

Let's first import a few libraries.

That's an impressive list of imports.

import numpy as np from numpy import linalg from numpy.linalg import norm from scipy.spatial.distance import squareform, pdist

We import sklearn.

import sklearn from sklearn.manifold import TSNE from sklearn.datasets import load_digits from sklearn.preprocessing import scale

We'll hack a bit with the t-SNE code in sklearn 0.15.2.

from sklearn.metrics.pairwise import pairwise_distances from sklearn.manifold.t_sne import (_joint_probabilities, _kl_divergence) from sklearn.utils.extmath import _ravel

Random state.

RS = 20150101

We'll use matplotlib for graphics.

import matplotlib.pyplot as plt import matplotlib.patheffects as PathEffects import matplotlib %matplotlib inline

We import seaborn to make nice plots.

import seaborn as sns sns.set_style('darkgrid') sns.set_palette('muted') sns.set_context("notebook", font_scale=1.5, rc={"lines.linewidth": 2.5})

We'll generate an animation with matplotlib and moviepy.

from moviepy.video.io.bindings import mplfig_to_npimage import moviepy.editor as mpy

Now we load the classic handwritten digits datasets. It contains 1797 images with \(8*8=64\) pixels each.

digits = load_digits() digits.data.shape

print(digits['DESCR'])

Here are the images:

nrows, ncols = 2, 5 plt.figure(figsize=(6,3)) plt.gray() for i in range(ncols * nrows): ax = plt.subplot(nrows, ncols, i + 1) ax.matshow(digits.images[i,...]) plt.xticks([]); plt.yticks([]) plt.title(digits.target[i]) plt.savefig('images/digits-generated.png', dpi=150)

Now let's run the t-SNE algorithm on the dataset. It just takes one line with scikit-learn.

We first reorder the data points according to the handwritten numbers.

X = np.vstack([digits.data[digits.target==i] for i in range(10)]) y = np.hstack([digits.target[digits.target==i] for i in range(10)])

digits_proj = TSNE(random_state=RS).fit_transform(X)

Here is a utility function used to display the transformed dataset. The color of each point refers to the actual digit (of course, this information was not used by the dimensionality reduction algorithm).

def scatter(x, colors): # We choose a color palette with seaborn. palette = np.array(sns.color_palette("hls", 10))

# We create a scatter plot.
f = plt.figure(figsize=(8, 8))
ax = plt.subplot(aspect='equal')
sc = ax.scatter(x[:,0], x[:,1], lw=0, s=40,
                c=palette[colors.astype(np.int)])
plt.xlim(-25, 25)
plt.ylim(-25, 25)
ax.axis('off')
ax.axis('tight')

# We add the labels for each digit.
txts = []
for i in range(10):
    # Position of each label.
    xtext, ytext = np.median(x[colors == i, :], axis=0)
    txt = ax.text(xtext, ytext, str(i), fontsize=24)
    txt.set_path_effects([
        PathEffects.Stroke(linewidth=5, foreground="w"),
        PathEffects.Normal()])
    txts.append(txt)

return f, ax, sc, txts

Here is the result.

scatter(digits_proj, y) plt.savefig('images/digits_tsne-generated.png', dpi=120)

We observe that the images corresponding to the different digits are clearly separated into different clusters of points.

Mathematical framework

Let's explain how the algorithm works. First, a few definitions.

A data point is a point \(x_i\) in the original data space \(\mathbf{R}^D\), where \(D=64\) is the dimensionality of the data space. Every point is an image of a handwritten digit here. There are \(N=1797\) points.

A map point is a point \(y_i\) in the map space \(\mathbf{R}^2\). This space will contain our final representation of the dataset. There is a bijection between the data points and the map points: every map point represents one of the original images.

How do we choose the positions of the map points? We want to conserve the structure of the data. More specifically, if two data points are close together, we want the two corresponding map points to be close too. Let's \(\left| x_i - x_j \right|\) be the Euclidean distance between two data points, and \(\left| y_i - y_j \right|\) the distance between the map points. We first define a conditional similarity between the two data points:

\(p_{j|i} = \frac{\exp\left(-\left| x_i - x_j\right|^2 \big/ 2\sigma_i^2\right)}{\displaystyle\sum_{k \neq i} \exp\left(-\left| x_i - x_k\right|^2 \big/ 2\sigma_i^2\right)}\)

This measures how close \(x_j\) is from \(x_i\), considering a Gaussian distribution around \(x_i\) with a given variance \(\sigma_i^2\). This variance is different for every point; it is chosen such that points in dense areas are given a smaller variance than points in sparse areas. The original paper details how this variance is computed exactly.

Now, we define the similarity as a symmetrized version of the conditional similarity:

\(p_{ij} = \frac{p_{j|i} + p_{i|j}}{2N}\)

We obtain a similarity matrix for our original dataset. What does this matrix look like?

Similarity matrix

The following function computes the similarity with a constant \(\sigma\).

def _joint_probabilities_constant_sigma(D, sigma): P = np.exp(-D2/2 * sigma2) P /= np.sum(P, axis=1) return P

We now compute the similarity with a \(\sigma_i\) depending on the data point (found via a binary search, according to the original t-SNE paper). This algorithm is implemented in the _joint_probabilities private function in scikit-learn's code.

Pairwise distances between all data points.

D = pairwise_distances(X, squared=True)

Similarity with constant sigma.

P_constant = _joint_probabilities_constant_sigma(D, .002)

Similarity with variable sigma.

P_binary = _joint_probabilities(D, 30., False)

The output of this function needs to be reshaped to a square matrix.

P_binary_s = squareform(P_binary)

We can now display the distance matrix of the data points, and the similarity matrix with both a constant and variable sigma.

plt.figure(figsize=(12, 4)) pal = sns.light_palette("blue", as_cmap=True)

plt.subplot(131) plt.imshow(D[::10, ::10], interpolation='none', cmap=pal) plt.axis('off') plt.title("Distance matrix", fontdict={'fontsize': 16})

plt.subplot(132) plt.imshow(P_constant[::10, ::10], interpolation='none', cmap=pal) plt.axis('off') plt.title("$p_{j|i}$ (constant $\sigma$)", fontdict={'fontsize': 16})

plt.subplot(133) plt.imshow(P_binary_s[::10, ::10], interpolation='none', cmap=pal) plt.axis('off') plt.title("$p_{j|i}$ (variable $\sigma$)", fontdict={'fontsize': 16}) plt.savefig('images/similarity-generated.png', dpi=120)

We can already observe the 10 groups in the data, corresponding to the 10 numbers.

Let's also define a similarity matrix for our map points.

\(q_{ij} = \frac{f(\left| x_i - x_j\right|)}{\displaystyle\sum_{k \neq i} f(\left| x_i - x_k\right|)} \quad \textrm{with} \quad f(z) = \frac{1}{1+z^2}\)

This is the same idea as for the data points, but with a different distribution (t-Student with one degree of freedom, or Cauchy distribution, instead of a Gaussian distribution). We'll elaborate on this choice later.

Whereas the data similarity matrix \(\big(p_{ij}\big)\) is fixed, the map similarity matrix \(\big(q_{ij}\big)\) depends on the map points. What we want is for these two matrices to be as close as possible. This would mean that similar data points yield similar map points.

A physical analogy

Let's assume that our map points are all connected with springs. The stiffness of a spring connecting points \(i\) and \(j\) depends on the mismatch between the similarity of the two data points and the similarity of the two map points, that is, \(p_{ij} - q_{ij}\). Now, we let the system evolve according to the laws of physics. If two map points are far apart while the data points are close, they are attracted together. If they are nearby while the data points are dissimilar, they are repelled.

The final mapping is obtained when the equilibrium is reached.

Here is an illustration of a dynamic graph layout based on a similar idea. Nodes are connected via springs and the system evolves according to law of physics (example by Mike Bostock).

Algorithm

Remarkably, this physical analogy stems naturally from the mathematical algorithm. It corresponds to minimizing the Kullback-Leiber divergence between the two distributions \(\big(p_{ij}\big)\) and \(\big(q_{ij}\big)\):

\(KL(P||Q) = \sum_{i, j} p_{ij} , \log \frac{p_{ij}}{q_{ij}}.\)

This measures the distance between our two similarity matrices.

To minimize this score, we perform a gradient descent. The gradient can be computed analytically:

\(\frac{\partial , KL(P || Q)}{\partial y_i} = 4 \sum_j (p_{ij} - q_{ij}) g\left( \left| x_i - x_j\right| \right) u_{ij} \quad \textrm{where} , g(z) = \frac{z}{1+z^2}.\)

Here, \(u_{ij}\) is a unit vector going from \(y_j\) to \(y_i\). This gradient expresses the sum of all spring for