Heavy-tailed kernels reveal a finer cluster structure in t-SNE visualisations

T-distributed stochastic neighbour embedding (t-SNE) is a widely used data visualisation technique. It differs from its predecessor SNE by the low-dimensional similarity kernel: the Gaussian kernel was replaced by the heavy-tailed Cauchy kernel, solving the"crowding problem"of SNE. Here, we develop an efficient implementation of t-SNE for a $t$-distribution kernel with an arbitrary degree of freedom $\nu$, with $\nu\to\infty$ corresponding to SNE and $\nu=1$ corresponding to the standard t-SNE. Using theoretical analysis and toy examples, we show that $\nu<1$ can further reduce the crowding problem and reveal finer cluster structure that is invisible in standard t-SNE. We further demonstrate the striking effect of heavier-tailed kernels on large real-life data sets such as MNIST, single-cell RNA-sequencing data, and the HathiTrust library. We use domain knowledge to confirm that the revealed clusters are meaningful. Overall, we argue that modifying the tail heaviness of the t-SNE kernel can yield additional insight into the cluster structure of the data.


Introduction
T-distributed stochastic neighbour embedding (t-SNE) [12] and related methods [15,13] are used for data visualisation in many scientific fields dealing with thousands or even millions of high-dimensional samples. They range from single-cell cytometry [1] and transcriptomics [16,19], where samples are cells and features are proteins or genes, to population genetics [4], where samples are people and features are single-nucleotide polymorphisms, to humanities [14], where samples are books and features are words.
T-SNE was developed from an earlier method called SNE [5]. The central idea of SNE was to describe pairwise relationships between high-dimensional points in terms of normalised affinities: close neighbours have high affinity whereas distant samples have near-zero affinity. SNE then positions the points in two dimensions such that the Kullback-Leibler divergence between the high-and lowdimensional affinities is minimised. This worked to some degree but suffered from what was later called the "crowding problem". The idea of t-SNE was to adjust the kernel transforming pairwise low-dimensional distances into affinities: the Gaussian kernel was replaced by the heavy-tailed Cauchy kernel (t-distribution with one degree of freedom ν), ameliorating the crowding problem.
The choice of the specific heavy-tailed kernel was mostly motivated by mathematical and computational simplicity: a t-distribution with ν = 1 has a density proportional to 1/(1+x 2 ) which is mathematically compact and fast to compute. However, a t-distribution with any finite ν has heavier tails than the Gaussian distribution (which corresponds to ν → ∞). It is therefore reasonable to explore the whole spectrum of the values of ν from ∞ to 0. Given that t-SNE (ν = 1) outperforms SNE (ν = ∞), it might be that for some data sets ν < 1 would offer additional insights into the structure of the data.
While this seems like a straightforward extension and has already been discussed in the literature [10,18], no efficient implementation of this idea has been available until now. T-SNE is usually optimised via adaptive gradient descent. While it is easy to write down the gradient for an arbitrary value of ν, the exact t-SNE from the original paper requires O(n 2 ) time and memory, and cannot be run for sample sizes much larger than n = 10 000. Efficient approximations have been developed allowing to run approximate t-SNE for much larger sample sizes [11,9], but until now have only been implemented for ν = 1. As a result, the effect of ν = 1 on large real-life datasets has remained unknown.
Here we show that the recent FIt-SNE approximation [9] can be modified to use an arbitrary value of ν and demonstrate that ν < 1 can reveal "hidden" structure, invisible with standard t-SNE.

t-SNE with arbitrary degree of freedom
SNE defines directional affinity of point x j to point x i as .
For each i, this forms a probability distribution over all points j = i (all p i|i are set to zero). The variance of the Gaussian kernel σ 2 i is chosen such that the perplexity of this probability distribution exp − ln(2) · j =i p j|i log 2 p j|i has some pre-specified value. In symmetric SNE (SSNE) 5 and t-SNE the affinities are symmetrised and normalised to form a probability distribution on the set of all pairs (i, j).
The points are then arranged in a low-dimensional space to minimise the Kullback-Leibler (KL) divergence between p ij and the affinities in the lowdimensional space, q ij : Here k(d) is a kernel that transforms Euclidean distance d between any two points into affinities, and y i are low-dimensional coordinates (all q ii are set to 0). SNE uses the Gaussian kernel k(d) = exp(−d 2 ). T-SNE uses the t-distribution with one degree of freedom (also known as Cauchy distribution): k(d) = 1/(1 + d 2 ). Here we consider a general t-distribution kernel ( ) As in [18], we use a simplified version defined as ( ) This kernel corresponds to the scaled t-distribution with ν = 2α − 1. This means that using ( ) instead of ( ) in t-SNE produces an identical output apart from the global scaling by 2ν/(ν + 1). At the same time, ( ) allows to use any α > 0, including α ∈ (0, 1/2] corresponding to negative ν, i.e. it allows kernels with tails heavier than any possible t-distribution. 6 Yang et al. [18] use the same kernel but with α replaced by 1/α, and call it "heavy-tailed SNE" (HSSNE). The gradient of the loss function (see Appendix or [18]) is Any implementation of exact t-SNE can be easily modified to use this expression instead of the α = 1 gradient. 5 In the following text we will not make a distinction between the symmetric SNE (SSNE) and the original, asymmetric, SNE. 6 Equivalently, we could use an even simpler kernel k(d) = (1 + d 2 ) −α that differs from ( ) only by scaling. We prefer ( ) because of the explicit Gaussian limit at α → ∞. Modern t-SNE implementations make two approximations. First, they set most p ij to zero, apart from only a small number of close neighbours [11,9], accelerating the attractive force computations (that can be very efficiently parallelised). This carries over to the α = 1 case. The repulsive forces are approximated in FIt-SNE by interpolation on a grid, further accelerated with the Fourier transform [9]. This interpolation can be carried out for the α = 1 case in full analogy to the α = 1 case (see Appendix).
Importantly, the runtime of FIt-SNE with α = 1 is practically the same as with α = 1. For example, embedding MNIST (n = 70 000) with perplexity 50 as described below took 90 seconds with α = 1 and 97 seconds with α = 0.5 on a computer with 4 double-threaded cores, 3.4 GHz each. 7

Toy examples
We first applied exact t-SNE with various values of α to a simple toy data set consisting of several well-separated clusters. Specifically, we generated a 10dimensional data set with 100 data points in each of the 10 classes (1000 points overall). The points in class i were sampled from a Gaussian distribution with covariance I 10 and mean µ i = 4e i where e i is the i-th basis vector. We used perplexity 50, and default optimisation parameters (1000 iterations, learning rate 200, early exaggeration 12, length of early exaggeration 250, initial momentum 0.5, switching to 0.8 after 250 iterations).  Figure 1 shows the t-SNE results for α = 100, α = 1, and α = 0.1. A tdistribution with ν = 2α − 1 = 199 degrees of freedom is very close to the Gaussian distribution, so here and below we will refer to the α = 100 result as SNE. We see that class separation monotonically increases with decreasing α: t-SNE ( Figure 1B) separates the classes much better than SNE ( Figure 1A), but t-SNE with α = 0.5 separates them much better still ( Figure 1C).
In the above toy example, the choice between different values of α is mostly aesthetic. This is not the case in the next toy example. Here we change the dimensionality to 20 and shift 50 points in each class by 2e 10+i and the remaining 50 points by −2e 10+i (where i is the class number). The intuition is that now each of the 10 classes has a "dumbbell" shape. This shape is invisible in SNE ( Figure 2A) and hardly visible in standard t-SNE ( Figure 2B), but becomes apparent with α = 0.5 ( Figure 2C). In this case, decreasing α below 1 is necessary to bring out the fine structure of the data.

Mathematical analysis
We showed that decreasing α increases cluster separation (Figures 1, 2). Why does this happen? An informal argument is that in order to match the betweencluster affinities p ij , the distance between clusters in the t-SNE embedding needs to grow when the kernel becomes progressively more heavy-tailed [12].
To quantify this effect, we consider an example of two standard Gaussian clusters in 10 dimensions (n = 100 in each) with the between-centroid distance set to 5 √ 2; these clusters can be unambiguously separated. We use exact t-SNE (perplexity 50) with various values of α from 0.2 to 3.0 and measure the cluster separation in the embedding. As a scale-invariant measure of separation we used between-centroids distance divided by the root-mean-square within- Between/within distance ratio Fig. 3. Separation in the t-SNE visualisation between the two well-separated clusters as a function of α. Separation was measured as the between-centroids distance divided by the root-mean-square within-cluster distance.
cluster distance. Indeed, we observed a monotonic decrease of this measure with growing α (Figure 3).
The informal argument mentioned above can be replaced by the following formal one. Consider two high-dimensional clusters (n points in each) with all pairwise within-cluster distances equal to D w and all pairwise between-cluster distances equal to D b D w (this can be achieved in the space of 2n dimensions). In this case, the p ij matrix has only two unique non-zero values: all within-cluster affinities are given by p w and all between-cluster affinities by p b , where K(D) is the Gaussian kernel corresponding to the chosen perplexity value. Consider an exact t-SNE mapping to the space of the same dimensionality. In this idealised case, t-SNE can achieve zero loss by setting within-and betweencluster distances d w and d b in the embedding such that q w = p w and Plugging in the expression for k(d) and denoting the constant right-hand side by c < 1, we obtain The left-hand side can be seen as a measure of class separation close to the one used in Figure 3, and the right-hand side monotonically decreases with increasing α.
In the simulation shown in Figure 3, the p ij matrix does not have only two unique elements, the target dimensionality is two, and the t-SNE cannot possibly achieve zero loss. Still, qualitatively we observe the same behaviour: approximately power-law decrease of separation with increasing α.

Real-life data sets
We now demonstrate that these theoretical insights are relevant to practical use cases on large-scale data sets. Here we use approximate t-SNE (FIt-SNE).
MNIST We applied t-SNE with various values of α to the MNIST data set (Figure 4), comprising n = 70 000 grayscale 28 × 28 images of handwritten digits. As a pre-processing step, we used principal component analysis (PCA) to reduce the dimensionality from 784 to 50. We used perplexity 50 and default optimisation parameters apart from learning rate that we increased to η = 1000. 8 For easier reproducibility, we initialised the t-SNE embedding with the first two PCs (scaled such that PC1 had standard deviation 0.0001).
To the best of our knowledge, Figure 4A is the first existing SNE (α = 100) visualisation of the whole MNIST: we are not aware of any SNE implementation 8 To get a good t-SNE visualisation of MNIST, it is helpful to increase either the learning rate or the length of the early exaggeration phase. Default optimisation parameters often lead to some of the digits being split into two clusters. In the cytometric context, this phenomenon was described in detail by [2].
that can handle a dataset of this size. It produces a surprisingly good visualisation but is nevertheless clearly outperformed by standard t-SNE (α = 1, Figure 4B): many digits coalesce together in SNE but get separated into clearly distinct clusters in t-SNE. Remarkably, reducing α to 0.5 makes each digit further split into multiple separate sub-clusters ( Figure 4C), revealing a fine structure within each of the digits.
To demonstrate that these sub-clusters are meaningful, we computed the average MNIST image for some of the sub-clusters ( Figure 4D). In each case, the shapes appear to be meaningfully distinct: e.g. for the digit "4", the handwriting is more italic in one sub-cluster, more wide in another, and features a non-trivial homotopy group (i.e. has a loop) in yet another one. Similarly, digit "2" is separated into three sub-clusters, with the most abundant one showing a loop in the bottom-left and the next abundant one having a sharp angle instead. Digit "1" is split according to the stroke angle. Re-running t-SNE using random initialisation with different seeds yielded consistent results. Points that appear as outliers in Figure 4C mostly correspond to confusingly written digits.
MNIST has been a standard example for t-SNE starting from the original t-SNE paper [12], and it has been often observed that t-SNE preserves meaningful within-digit structure. Indeed, the sub-clusters that we identified in Figure 4C are usually close together in Figure 4B. 9 However, standard t-SNE does not separate them into visually isolated sub-clusters, and so does not make this internal structure obvious.
Single-cell RNA-sequencing data For the second example, we took the transcriptomic dataset from [16], comprising n = 23 822 cells from adult mouse cortex (sequenced with Smart-seq2 protocol). Dimensions are genes, and the data are the integer counts of RNA transcripts of each gene in each cell. Using a custom expert-validated clustering procedure, the authors divided these cells into 133 clusters. In Figure 5, we used the cluster ids and cluster colours from the original publication. Figure 5A shows the standard t-SNE (α = 1) of this data set, following common transcriptomic pre-processing steps as described in [7]. Briefly, we rownormalised and log-transformed the data, selected 3000 most variable genes and used PCA to further reduce dimensionality to 50. We used perplexity 50 and PCA initialisation. The resulting t-SNE visualisation is in a reasonable agreement with the clustering results, however it lumps many clusters together into contiguous "islands" or "continents" and overall suggests many fewer than 133 distinct clusters.
Reducing the number of degrees of freedom to α = 0.6 splits many of the contiguous islands into "archipelagos" of smaller disjoint areas ( Figure 5B). In many cases, this roughly agrees with the clustering results of [16]. Figure 5C shows a zoom-in into the Vip clusters (west-southwest part of panel B) that provide one such example: isolated islands correspond well to the individual clusters (or sometimes pairs of clusters). Importantly, the cluster labels in this data set are not ground truth; nevertheless the agreement between cluster labels and t-SNE with α = 0.6 provides additional evidence that this data categorisation is meaningful.
HathiTrust library For the final example, we used the HathiTrust library data set [14]. The full data set comprises 13.6 million books and can be described with several million features that represent word counts of each word in each book. We used the pre-processed data from [14]: briefly, the word counts were row-normalised, log-transformed, projected to 1280 dimensions using random linear projection with coefficients ±1, and then reduced to 100 PCs. 10 The available meta-data include author name, book title, publication year, language, and Library of Congress classification (LCC) code. For simplicity, we took a n = 408 291 subset consisting of all books in Russian language. We used perplexity 50 and learning rate η = 10 000. Figure 6A shows the standard t-SNE visualisation (α = 1) coloured by the publication year. The most salient feature is that pre-1917 books cluster together (orange/red colours): this is due to the major reform of Russian orthography implemented in 1917, leading to most words changing their spelling. However, not much of a substructure can be seen among the books published after (or before) 1917. In contrast, t-SNE visualisation with α = 0.5 fragments the corpus into a large number of islands ( Figure 6B). We can identify some of the islands by inspecting the available meta-data. For example, mathematical literature (LCC code QA, n = 6490 books) is not separated from the rest in standard t-SNE, but occupies the leftmost island in t-SNE with α = 0.5 (blue contour lines in both panels). Several neighbouring islands correspond to the physics literature (LCC code QC, n = 5104 books; not shown). In an attempt to capture something radically different from mathematics, we selected all books authored by several famous Russian poets 11 (n = 1369 in total). This is not a curated list: there are non-poetry books authored by these authors, while many other poets were not included (the list of poets was not cherry-picked; we made the list before looking at the data). Nevertheless, when using α = 0.5, the poetry books printed after 1917 seemed to occupy two neighbouring islands, and the ones printed before 1917 were reasonably isolated as well ( Figure 6B, black contour lines). In the standard t-SNE visualisation poetry was not at all separated from the surrounding population of books.

Related work
Yang et al. [18] introduced symmetric SNE with the kernel family k(d) = 1 (1 + αd 2 ) 1/α , calling it "heavy-tailed symmetric SNE" (HSSNE). This is exactly the same kernel family as ( ), but with α replaced by 1/α. However, Yang et al. did not show any examples of heavier-tailed kernels revealing additional structure compared to α = 1 and did not provide an implementation suitable for large sample sizes. Interestingly, Yang et al. argued that gradient descent is not suitable for HSSNE and suggested an alternative optimisation algorithm; here we demonstrated that the standard t-SNE optimisation works reasonably well in a wide range of α values (but see Discussion).
Van der Maaten [10] discussed the choice of the degree of freedom in the t-distribution kernel in the context of parametric t-SNE. He argued that ν > 1 might be warranted when embedding the data in more than two dimensions. He also implemented a version of parametric t-SNE that optimises over ν. However, similar to [18], [10] did not contain any examples of ν < 1 being actually useful in practice.
UMAP [13] is a promising recent algorithm closely related to an earlier largeVis [15]; both are similar to t-SNE but modify the repulsive forces to make them amenable for a sampling-based stochastic optimisation. UMAP uses the following family of similarity kernels: which reduces to Cauchy when a = b = 1 and is more heavy-tailed when 0 < b < 1. UMAP default is a ≈ 1.6 and b ≈ 0.9 with both parameters adjusted via the min dist input parameter (default value 0.1). Decreasing min dist all the way to zero corresponds to decreasing b to 0.79. In our experiments, we observed that modifying min dist (or b directly) led to an effect qualitatively similar to modifying α in t-SNE. For some data sets this required manually decreasing b below 0.79. In case of MNIST, b = 0.3, but not b = 0.79, revealed sub-digit structure ( Figure S1) -an effect that has not been described before (cf. [13] where McInnes et al. state that min dist is "an essentially aesthetic parameter"). In other words, the same conclusion seems to apply to UMAP: heavy-tailed kernels reveal a finer cluster structure. A more in-depth study of the relationships between the two algorithms is beyond the scope of this paper. The horizontal axis is on the log scale. The α values were sampled on a grid with step 0.01 for α < 1, 0.25 for 1 ≤ α ≤ 5 and 1 for α > 5. The black line shows KL divergence (left axis) with minimum at α = 1.5. Running gradient descent with α = 0.5 for 10 000 iterations ( Figure S3) lowered KL divergence down to 3.6, which was still above the minimum value. Blue lines show neighbourhood preservation (the fraction of k nearest neighbours of each point that remain within k nearest neighbours in the embedding, averaged over all n = 70 000 points) for k = 10, k = 50, and k = 100. and we are not aware of any a priori argument in favour of α = 1. As α = 1 still yields a t-distribution kernel (scaled t-distribution to be precise), we prefer not to use a separate acronym (HSSNE [18]). If needed, one can refer to t-SNE with α < 1 as "heavy-tailed" t-SNE.
We found that lowering α below 1 makes progressively finer structure apparent in the visualisation and brings out smaller clusters, which -at least in the data sets studied here -are often meaningful. In a way, α < 1 can be thought of as a "magnifying glass" for the standard t-SNE representation. We do not think that there is one ideal value of α suitable for all data sets and all situations; instead we consider it a useful adjustable parameter of t-SNE, complementary to the perplexity. We observed a non-trivial interaction between α and perplexity: Small vs. large perplexity makes the affinity matrix p ij represent the local vs. global structure of the data [7]. Small vs. large α makes the embedding represent the finer vs. coarser structure of the affinity matrix. In practice, it can make sense to treat it as a two-dimensional parameter space to explore. However, for large data sets (n 10 6 ), it is computationally unfeasible to substantially increase the perplexity from its standard range of 30-100 (as it would prohibitively increase the runtime), and so α becomes the only available parameter to adjust.
One important caveat is to be kept in mind. It is well-known that t-SNE, especially with low perplexity, can find "clusters" in pure noise, picking up random fluctuations in the density [17]. This can happen with α = 1 but gets exacerbated with lower values of α. A related point concerns clustered real-life data where separate clusters (local density peaks) can sometimes be connected by an area of lower but non-zero density: for example, [16] argued that many pairs of their 133 clusters have intermediate cells. Our experiments demonstrate that lowering α can make such clusters more and more isolated in the embedding, creating a potentially misleading appearance of perfect separation (see e.g. Figure 1). In other words, there is a trade-off between bringing out finer cluster structure and preserving continuities between clusters.
Choosing a value of α that yields the most faithful representation of a given data set is challenging because it is difficult to quantify "faithfulness" of any given embedding [8]. For example, for MNIST, KL divergence is minimised at α ≈ 1.5 (Figure 7), but it may not be the ideal metric to quantify the embedding quality [6]. Indeed, we found that k-nearest neighbour (KNN) preservation [8] peaked elsewhere: the peak for k = 10 was at α ≈ 1.0, for k = 50 at α ≈ 0.9, and for k = 100 at α ≈ 0.8 (Figure 7). We stress that we do not think that KNN preservation is the most appropriate metric here; our point is that different metrics can easily disagree with each other. In general, there may not be a single "best" embedding of high-dimensional data in a two-dimensional space. Rather, by varying α, one can obtain different complementary "views" of the data.
Very low values of α correspond to kernels with very wide and very flat tails, leading to vanishing gradients and difficult convergence. We found that α = 0.5 was about the smallest value that could be safely used ( Figure S2). In fact, it may take more iterations to reach convergence for 0.5 < α < 1 compared to α = 1. As an example, running t-SNE on MNIST with α = 0.5 for ten times longer than we did for Figure 4C, led to the embedding expanding much further (which leads to a slow-down of FIt-SNE interpolation) and, as a result, resolving additional subclusters ( Figure S3). On a related note, when using only one single MNIST digit as an input for t-SNE with α = 0.5, the embedding also fragments into many more clusters ( Figure S4), which we hypothesise is due to the points rapidly expanding to occupy a much larger area compared to what happens in the full MNIST embedding ( Figure S4). This can be counterbalanced by increasing the strength of the attractive forces ( Figure S4). Overall, the effect of the embedding scale on the cluster resolution remains an open research question.
In conclusion, we have shown that adjusting the heaviness of the kernel tails in t-SNE can be a valuable tool for data exploration and visualisation. As a practical recommendation, we suggest to embed any given data set using various values of α, each inducing a different level of clustering, and hence providing insight that cannot be obtained from the standard α = 1 choice alone. 12

Appendix
The loss function, up to a constant term p ij log p ij , can be rewritten as follows: where we took into account that p ij = 1. The first term in Eq. (1) contributes attractive forces to the gradient while the second term yields repulsive forces. The gradient is The first expression is more convenient for numeric optimisation while the second one can be more convenient for mathematical analysis. For the kernel Plugging Eq. 4 into Eq. 3, we obtain the expression for the gradient [18] 13 For numeric optimisation it is convenient to split this expression into the attractive and the repulsive terms. Plugging Eq. 4 into Eq. 2, we obtain It is noteworthy that the expression for F attr has w ij raised to the 1/α power, which cancels out the fractional power in k(d). This makes the runtime of F attr computation unaffected by the value of α. In FIt-SNE, the sum over j in F attr is approximated by the sum over 3Π approximate nearest neighbours of point i obtained using Annoy [3], where Π is the provided perplexity value. The 3Π heuristic comes from [11]. The remaining p ij values are set to zero. The F rep can be approximated using the interpolation scheme from [9]. It allows fast approximate computation of the sums of the form j K( y i − y j ) and j K( y i − y j )y j , where K(·) is any smooth kernel, by using polynomial interpolation of K on a fine grid. 14 All kernels appearing in F rep are smooth.   . Toy example with ten "dumbbell"-shaped clusters from Figure 2, here embedded with α = 0.1. Top-left plot shows the result after 1000 gradient descent iterations (default). Note that the dumbbell shape is lost: whereas the number of visible clusters increased as α was lowered from 100 to 0.5 (Figure 2), it decreased when it was further lowered to 0.1. We believe the reason for this is that the strong repulsion between dumbbells "squashes" them in the beginning of optimisation into very compact blobs. It is likely that longer optimisation would resolve the dumbbell shapes. This is difficult to test because the kernel with α = 0.1 is extremely wide and flat, leading to slow convergence. Top-right plot shows the result after 5000 iterations. Here a few outlying points get pushed to the periphery. Zooming-in to the main 10 clusters (bottom-left) still does not resolve the dumbbell shapes. Further zooming in on one of the dumbbells (bottom-right) shows that the points are squashed into 1D which may be a sign of poor convergence.
In a separate sets of experiments, we observed the similar phenomenon with MNIST: α = 0.2 after 1000 iterations yielded fewer clusters than α = 0.5. Our conclusion is that smaller values of α should be used with caution.  Figure 4C; it was obtained with 1000 gradient descent iterations (the default value). The top-right panel corresponds to 10 000 iterations and has many more isolated sub-clusters. This can also be seen in the bottom row showing the respective zoom-ins into the digit "4". At the same time, the embedding after 1000 iterations is not misleading and is simply a coarser-grained version of the embedding after 10 000 iterations.
Using 10 000 iterations is impractical: whereas 1000 iterations were finished in 1.5 minutes, 10 000 iterations took 4 hours 30 minutes. This is because FIt-SNE interpolation scheme uses regular interpolation grid with the number of nodes growing quadratically with the embedding size. While the left embedding is contained within [−50, 50] 2 , the right one expands to [−200, 200] 2 . In principle, an implementation based on the fast multipole method (FMM) could be developed to dramatically accelerate the gradient computation in this setting where most of the embedding space is "empty space", but current FIt-SNE implementation does not support it.
Note that the standard t-SNE embedding with α = 1 also expands much further after 10 000 iterations, compared to the 1000 iterations. However, with α = 1 it does not resolve additional sub-clusters, at least in MNIST.  Figure 4C. Note the large number of isolated clusters. We believe this happens because the embedding rapidly expands to a larger area, compared to Figure S3 (bottom-left). One evidence for that is that re-running t-SNE after adding several random Gaussian clusters with n = 7 000 each, roughly recovers the shape of the digit "4" archipelago from the full MNIST embedding (middle). Right: α = 0.5 and exaggeration factor 1.5 [7], i.e. all attractive forces are multiplied by 1.5 after the end of the early exaggeration phase (during the early exaggeration they are multiplied by 12). This roughly recovers the sub-clusters from the full MNIST embedding ( Figure S3). The relationship between α and exaggeration remains for future work.