ActUp: Analyzing and Consolidating tSNE and UMAP

1 Introduction

Dimensionality Reduction (DR) algorithms are invaluable for qualitatively inspecting high-dimensional data and are widely used across scientific disciplines. Broadly speaking, these algorithms transform a high-dimensional input into a faithful lower-dimensional embedding. This embedding aims to preserve similarities among the points, where similarity is often measured by distances in the corresponding spaces.

tSNE [Van der Maaten and Hinton, 2008] [Van Der Maaten, 2014] and UMAP [McInnes et al., 2018] are two widely popular DR algorithms due to their efficiency and interpretable embeddings. Both algorithms establish analogous similarity measures, share comparable loss functions, and find an embedding through gradient descent. Despite these similarities, tSNE and UMAP have several key differences. First, although both methods obtain similar results, UMAP prefers large inter-cluster distances while tSNE leans towards large intra-cluster distances. Second, UMAP runs significantly faster as it performs efficient sampling during gradient descent. While attempts have been made to study the gaps between the algorithms [Kobak and Linderman, 2021], [Bohm et al., 2020], [Damrich and Hamprecht, 2021], there has not yet been a comprehensive analysis of their methodologies nor a method that can obtain both tSNE and UMAP embeddings at UMAP speeds.

We believe that this is partly due to their radically different presentations. While tSNE takes a computational angle, UMAP originates from category theory and topology. Despite this, many algorithmic choices in UMAP and tSNE are presented without theoretical justification, making it difficult to know which algorithmic components are necessary.

In this paper, we make the surprising discovery that the differences in both the embedding structure and computational complexity between tSNE and UMAP can be resolved via a single algorithmic choice - the normalization factor. We come to this conclusion by deriving both algorithms from first principles and theoretically showing the effect that the normalization has on the gradient structure. We supplement this by identifying every implementation and hyperparameter difference between the two methods and implementing tSNE and UMAP in a common library. Thus, we study the effect that each choice has on the embeddings and show both quantitatively and qualitatively that, other than the normalization of the pairwise similarity matrices, none of these parameters significantly affect the outputs.

Based on this analysis, we introduce the necessary changes to the UMAP algorithm such that it can produce tSNE embeddings as well. We refer to this algorithm as Gradient Dimensionality Reduction (GDR) to emphasize that it is consistent with the presentations of both tSNE and UMAP. We experimentally validate that GDR can simulate both methods through a thorough quantitative and qualitative evaluation across many datasets and settings. Lastly, our analysis provides insights for further speed improvements and allows GDR to perform gradient descent faster than the standard implementation of UMAP.

In summary, our contributions are as follows:

1. We perform the first comprehensive analysis of the differences between tSNE and UMAP, showing the effect of each algorithmic choice on the embeddings.
2. We theoretically and experimentally show that changing the normalization is a sufficient condition for switching between the two methods.
3. We release simple, plug-and-play implementations of GDR, tSNE and UMAP that can toggle all of the identified hyperparameters. Furthermore, GDR obtains embeddings for both algorithms faster than UMAP.

2 RelatedWork

When discussing tSNE we are referring to [Van Der Maaten, 2014] which established the nearest neighbor and sampling improvements and is generally accepted as the standard tSNE method. A popular subsequent development was presented in [Linderman et al., 2019], wherein Fast Fourier Transforms were used to accelerate the comparisons between points. Another approach is LargeVis [Tang et al., 2016], which modifies the embedding functions to satisfy a graph-based Bernoulli probabilistic model of the low-dimensional dataset. As the more recent algorithm, UMAP has not had as many variations yet. One promising direction, however, has extended UMAP's second step as a parametric optimization on neural network weights [Sainburg et al., 2020].

Many of these approaches utilize the same optimization structure where they iteratively attract and repel points. While most perform their attractions along nearest neighbors in the high-dimensional space, the repulsions are the slowest operation and each method approaches them differently. tSNE samples repulsions by utilizing Barnes-Hut (BH) trees to sum the forces over distant points. The work in [Linderman et al., 2019] instead calculates repulsive forces with respect to specifically chosen interpolation points, cutting down on the O(n log n) BH tree computations. UMAP and LargeVis, on the other hand, simplify the repulsion sampling by only calculating the gradient with respect to a constant number of points. These repulsion techniques are, on their face, incompatible with one another, i.e., several modifications have to be made to each algorithm before one can interchange the repulsive force calculations.

There is a growing amount of work that compares tSNE and UMAP through a more theoretical analysis [Damrich and Hamprecht, 2021][Bohm et al., 2020][Damrich et al., 2022][Kobak and Linderman, 2021]. [Damrich and Hamprecht, 2021] find that UMAP's algorithm does not optimize the presented loss and provide its effective loss function. Similarly [Bohm et al., 2020] analyze tSNE and UMAP through their attractive and repulsive forces, discovering that UMAP diverges when using O(n) repulsions per epoch. We expand on the aforementioned findings by showing that the forces are solely determined by the choice of normalization, giving a practical treatment to the proposed ideas. Lastly, [Damrich et al., 2022] provides the interesting realization that tSNE and UMAP can both be described through contrastive learning approaches. Our work differs from theirs in that we analyze the full space of parameters in the algorithms and distill the difference to a single factor, allowing us to connect the algorithms without the added layers of contrastive learning theory. Lastly, the authors in [Kobak and Linderman, 2021] make the argument that tSNE can perform UMAP's manifold learning if given UMAP's initialization. Namely, tSNE randomly initializes the low dimensional embedding whereas UMAP starts from a Laplacian Eigenmap [Belkin and Niyogi, 2003] projection. While this may help tSNE preserve the local kNN structure of the manifold, it is not true of the macrolevel distribution of the embeddings. Lastly, [Wang et al., 2021] discusses the role that the loss function has on the resulting embedding structure. This is in line with our results, as we show that the normalization's effect on the loss function is fundamental in the output differences between tSNE and UMAP.

3 Comparison of tSNE and UMAP

We begin by formally introducing the tSNE and UMAP algorithms. Let X ∈ Rn×D be a high dimensional dataset of n points and let Y ∈ Rn×d be a previously initialized set of n points in lower-dimensional space such that d < D. Our aim is to define similarity measures between the points in each space and then find the embedding Y such that the pairwise similarities in Y match those in X.

To do this, both algorithms define high- and lowdimensional non-linear functions p : X × X → [0, 1] and q : Y × Y → [0, 1]. These form pairwise similarity matrices P(X),Q(Y ) ∈ Rn×n , where the i, j-th matrix entry represents the similarity between points i and j. Formally,

where d(xi, xj) is the high-dimensional distance function, σ and τ are point-specific variance scalars 1 , ρi = minl̸=i d(xi, xl), and a and b are constants. Note that the tSNE denominators in Equation 1 are the sums of all the numerators. We thus refer to tSNE's similarity functions as being normalized while UMAP's are unnormalized.

The high-dimensional p values are defined with respect to the point in question and are subsequently symmetrized. WLOG, let pij = S(pj|i, pi|j) for some symmetrization function S. Going forward, we write pij and qij without the superscripts when the normalization setting is clear from the context.

Given these pairwise similarities in the high- and lowdimensional spaces, tSNE and UMAP attempt to find the embedding Y such that Q(Y ) is closest to P(X). Since both similarity measures carry a probabilistic interpretation, we find an embedding by minimizing the KL divergence KL(P∥Q). This gives us:

In essence, tSNE minimizes the KL divergence of the entire pairwise similarity matrix since its P and Q matrices sum to 1. UMAP instead defines Bernoulli probability distributions {pij, 1 - pij}, {qij, 1 - qij} and sums the KL divergences between the n2 pairwise probability distributions 2 .

3.1 Gradient Calculations

We now describe and analyze the gradient descent approaches in tSNE and UMAP. First, notice that the gradients of each algorithm change substantially due to the differing normalizations. In tSNE, the gradient can be written as an attractive AtSNE i and a repulsive Rtsne i force acting on point yi with

where Z is the normalization term in qtSNE ij . On the other hand, UMAP's attractions and repulsions 3 are presented as [McInnes et al., 2018]

1 In practice, we can assume that 2σ2 i is functionally equivalent to τi, as they are both chosen such that the entropy of the resulting distribution is equivalent.

2 Both tSNE and UMAP set the diagonals of P and Q to 0

3 The ε value is only inserted for numerical stability

In the setting where a = b = 1 and ε = 0, Equations 6, 7 can be written as 4

We remind the reader that we are overloading notation - p and q are normalized when they are in the tSNE setting and are unnormalized in the UMAP setting.

In practice, tSNE and UMAP optimize their loss functions by iteratively applying these attractive and repulsive forces. It is unnecessary to calculate each such force to effectively estimate the gradient, however, as the pij term in both the tSNE and UMAP attractive forces decays exponentially. Based on this observation, both methods establish a nearest neighbor graph in the high-dimensional space, where the edges represent nearest neighbor relationships between xi and xj. It then suffices to only perform attractions between points yi and yj if their corresponding xi and xj are nearest neighbors.

This logic does not transfer to the repulsions, however, as the Student-t distribution has a heavier tail so repulsions must be calculated evenly across the rest of the points. tSNE does this by fitting a Barnes-Hut tree across Y during every epoch. If yk and yl are both in the same tree leaf then we assume qik = qil, allowing us to only calculate O(log(n)) similarities. Thus, tSNE estimates all n - 1 repulsions by performing one such estimate for each cell in Y 's Barnes-Hut tree. UMAP, on the other hand, simply obtains repulsions by sampling a constant number of points uniformly and only applying those repulsions. These repulsion schemas are depicted in Figure 3. Note, tSNE collects all of the gradients before a full momentum gradient descent step whereas UMAP moves each point immediately upon calculating a force.

There are a few differences between the two algorithms' gradient descent loops. First, the tSNE learning rate stays constant over training while UMAP's linearly decreases. Second, tSNE's gradients are strengthened by adding a "gains" term which scales gradients based on whether they point in the same direction from epoch to epoch 5 . We refer to these two elements as gradient amplification.

4 We derive this in section A.2 in the supplementary material

Note that UMAP's repulsive force has a 1 - pik term that is unavailable at runtime, as xi and xk may not have been nearest neighbors. In practice, UMAP estimates these 1 - pik terms by using the available pij values 6 . We also note that UMAP does not explicitly multiply by pij and 1 - pik. Instead, it samples the forces proportionally to these scalars. For example, if pij = 0.1 then we apply that force without the pij multiplier once every ten epochs. We refer to this as scalar sampling.

3.2 The Choice of Normalization

We now present a summary of our theoretical results before providing their formal statements. As was shown in [Bohm et al., 2020], the ratio between attractive and repulsive magnitudes determines the structure of the resulting embedding. Given this context, Theorem 1 shows that the normalization directly changes the ratio of attraction/repulsion magnitudes, inducing the difference between tSNE and UMAP embeddings. Thus, we can toggle the normalization to alternate between their outputs. Furthermore, Theorem 2 shows that the attraction/repulsion ratio in the normalized setting is independent of the number of repulsive samples collected. This second point allows us to accelerate tSNE to UMAP speeds without impacting embedding quality by simply removing the dependency on Barnes-Hut trees and calculating 1 per-point repulsion as in UMAP. We now provide the necessary definitions for the theorems.

Assume that the pij terms are given. We now consider the dataset Y probabilistically by defining a set of random variables vij = yi - yj and assume that all O(n2 ) vij vectors are i.i.d. around a non-zero mean. Let rij = (1 + |vij|2 )-1 and define Z = ∑n2 i,j rij as the sum over n2 pairs of points and ˜ Z = ∑ni,j rij as the sum over n pairs of points. Then applying n per-point repulsions gives us the force acting on point yi of E[|Rtsne |] = E[∑nj ||(r2 ij /Z 2 ) · vij||]. We now define an equivalent force term in the setting where we have 1 per-point repulsion: E[| ˜Rtsne |] = E[||(r2 ij / ˜ Z 2 )·vij||]. Note that we have a constant number c of attractive forces acting on each point, giving E[|Atsne |] = c · ptsne ij E[||(rij/Z) · vij||] and E[| ˜ Atsne |] = c · ptsne ij E[||(rij/ ˜ Z) · vij||].

Thus, |Atsne | and |Rtsne | represent the magnitudes of the forces when we calculate tSNE's O(n) per-point repulsions while | ˜ Atsne | and | ˜Rtsne | represent the forces when we have UMAP's O(1) per-point repulsions. Given this, we have the following theorems:

5 This term has not been mentioned in the literature but is present in common tSNE implementations.

6 When possible, we use index k to represent repulsions and j to represent attractions to highlight that pik is never calculated in UMAP. See Section A.6 in the supplementary material for details.

The proofs are given in Sections A.3 and A.4 of the supplementary material. We point out that ptsne ij is normalized over the sum of all cn attractions that are sampled, giving us the estimate ptsne ij ∼ 1/(cn). Theorem1's result is visualized in the gradient plots in Figure 2. There we see that, for nonnegligible values of d(xi, xj), the UMAP repulsions can be orders of magnitude larger than the corresponding tSNE ones, even when accounting for the magnitude of the attractions. Furthermore, Section 5 evidences that toggling the normalization is sufficient to switch between the algorithms' embeddings and that no other hyperparameter accounts for the difference in inter-cluster distances between tSNE and UMAP.

4 Unifying tSNE and UMAP

This leads us to GDR- a modification to UMAP that can recreate both tSNE and UMAP embeddings at UMAP speeds. We choose the general name Gradient Dimensionality Reduction to imply that it is both UMAP and tSNE.

Our algorithm follows the UMAP optimization procedure except that we (1) replace the scalar sampling by iteratively processing attractions/repulsions and (2) apply the gradients after having collected all of them, rather than immediately upon processing each one. The first change accommodates the gradients under normalization since the normalized repulsive forces do not have the 1 - pik term to which UMAP samples proportionally. The second change allows for performing momentum gradient descent for faster convergence in the normalized setting.

Since we follow the UMAP optimization procedure, GDR defaults to producing UMAP embeddings. In the case of replicating tSNE, we simply normalize the P and Q matrices and scale the learning rate. Although we only collect O(1) attractions and repulsions for each point, their magnitudes are balanced due to Theorems 1 and 2. We refer to GDR as GDRumap if it is in the unnormalized setting and as GDRtsne if it is in the normalized setting. We note that changing the normalization necessitates gradient amplification.

By allowing GDR to toggle the normalization, we are free to choose the simplest options across the other parameters. GDR therefore defaults to tSNE's asymmetric attraction and a and b scalars along with UMAP's distance-metric, initialization, nearest neighbors, and pij symmetrization.

The supplementary material provides some further information on the flexibility of GDR (A.1), such as an accelerated version of the algorithm where we modify the gradient formulation such that it is quicker to optimize. This change induces a consistent 2× speedup of GDR over UMAP. Despite differing from the true KL divergence gradient, we find that the resulting embeddings are comparable. Our repository also provides a CUDA kernel that calculates GDRumap and GDRtsne embeddings in a distributed manner on a GPU.

4.1 Theoretical Considerations

UMAP's theoretical framework identifies the existence of a locally-connected manifold in the high-dimensional space under the UMAP pseudo-distance metric ˜ d. This pseudodistance metric is defined such that the distance from point xj to xi is equal to ˜ d(xi, xj)=d(xi, xj)-minl̸=i d(xi, xl). Despite this being a key element of the UMAP foundation, we find that substituting the Euclidean distance for the pseudodistance metric seems to have no effect on the embeddings, as seen in Tables 3 and 4. It is possible that the algorithm's reliance on highly non-convex gradient descent deviates enough from the theoretical discussion that the pseudodistance metric loses its applicability. It may also be the case that this pseudo-distance metric, while insightful from a theoretical perspective, is not a necessary calculation in order to achieve the final embeddings.

Furthermore, many of the other differences between tSNE and UMAP are not motivated by the theoretical foundation of either algorithm. The gradient descent methodology is entirely heuristic, so any differences therein do not impact the theory. This applies to the repulsion and attraction sampling and gradient descent methods. Moreover, the highdimensional symmetrization function, embedding initialization, symmetric attraction, and a, b scalars can all be switched to their alternative options without impacting either method's consistency within its theoretical presentation. Thus, each of these heuristics can be toggled without impacting the embedding's interpretation, as most of them do not interfere with the theory and none affect the output.

We also question whether the choice of normalization is necessitated by either algorithm's presentation. tSNE, for example, treats the normalization of P and Q as an assumption and provides no further justification. In the case of UMAP, it appears that the normalization does not break the assumptions of the original paper [McInnes et al., 2018, Sec. 2,3]. We therefore posit that the interpretation of UMAP as finding the best fit to the high-dimensional data manifold extends to tSNE as well, as long as tSNE's gradients are calculated under the pseudo-distance metric in the high-dimensional space. We additionally theorize that each method can be paired with either normalization without contradicting the foundations laid out in its paper.

We evidence the fact that tSNE can preserve manifold structure at least as well as UMAP in Table 2, where BarnesHut tSNE without normalization cleanly maintains the structure of the Swiss Roll dataset. We further discuss these manifold learning claims in the supplementary material (A.5).

For all of these reasons, we make the claim that tSNE and UMAP are computationally consistent with one another. That is, we conjecture that, up to minor changes, one could have presented UMAP's theoretical foundation and implemented it with the tSNE algorithm or vice-versa.

4.2 Frobenius Norm for UMAP

Finally, even some of the standard algorithmic choices can be modified without significantly impacting the embeddings. For example, UMAP and tSNE both optimize the KL divergence, but we see no reason that the Frobenius norm cannot be substituted in its place. Interestingly, the embeddings in Figure 8 in the supplementary material show that optimizing the Frobenius norm in the unnormalized setting produces outputs that are indistinguishable from the ones obtained by minimizing the KL-divergence. To provide a possible indication as to why this occurs, Figure 2 shows that the zerogradient areas between the KL divergence and the Frobenius norm strongly overlap, implying that a local minimum under one objective satisfies the other one as well.

We bring this up for two reasons. First, the Frobenius norm is a significantly simpler loss function to optimize than the KL divergence due to its convexity. We hypothesize that there must be simple algorithmic improvements that can exploit this property. Further detail is given in Section A.7 in the supplementary material. Second, it is interesting to consider that even fundamental assumptions such as the objective function can be changed without significantly affecting the embeddings across datasets.

5 Results

Metrics. There is no optimal way to compare embeddings - an analysis at the point-level loses global information while studying macro-structures loses local information. To account for this, we employ separate metrics to study the embeddings at the micro- and macro-scales. Specifically, we use the kNN accuracy to analyze preservation of local neighborhoods as established in [Van Der Maaten et al., 2009] and the V-measure [Rosenberg and Hirschberg, 2007] to study the embedding's global cluster structures 7 .

5.1 Hyperparameter Effects

We first show that a majority of the differences between tSNE and UMAP do not significantly affect the embeddings. Specifically, Table 4 shows that we can vary the hyperparameters in Table 1 with negligible change to the embeddings of any discussed algorithm. Equivalent results on other datasets can be found in Tables 8 and 9 in the supplementary material. Furthermore, Table 3 provides quantitative evidence that the hyperparameters do not affect the embeddings across datasets; similarly, Table 9 in the supplementary material confirms this finding across algorithms.

Looking at Table 4, the initialization and the symmetric attraction induce the largest variation in the embeddings. For the initialization, the relative positions of clusters change but the relevant inter-cluster relationships remain consistent 8 . Enabling symmetric attraction attracts yj to yi when we attract yi to yj . Thus, switching from asymmetric to symmetric attraction functionally scales the attractive force by 2. This leads to tighter tSNE clusters that would otherwise be evenly spread out across the embedding, but does not affect UMAP significantly. We thus choose asymmetric attraction for GDR as it better recreates tSNE embeddings.

We show the effect of single hyperparameter changes for combinatorial reasons. However, we see no significant difference between changing one hyperparameter or any number of them. We also eschew including hyperparameters that have no effect on the embeddings and are the least interesting. These include the exact vs. approximate nearest neighbors, gradient clipping, and the number of epochs.

5.2 Effect of Normalization

Although Theorem 2 shows that we can take fewer repulsive samples without affecting the repulsion's magnitude, we must also verify that the angle of the repulsive force is preserved as well. Towards this end, we plot the average angle between the tSNE Barnes-Hut repulsions and the UMAP sampled repulsions in Figure 4. We see that, across datasets, the direction of the repulsion remains consistent throughout the optimization process. Thus, since both the magnitude and the direction are robust to the number of samples taken, we conclude that one can obtain tSNE embeddings with O(1) per-point repulsions.

We now show that toggling the normalization allows tSNE to simulate UMAP embeddings and vice versa. Table 2 shows exactly this. First note that tSNE in the unnormalized setting has significantly more separation between clusters in a manner similar to UMAP. The representations are fuzzier than the UMAP ones as we are still estimating O(n) repulsions, causing the embedding to fall closer to the mean of the multimodal datasets. To account for the n× more repulsions, we scale each repulsion by 1/n for the sake of convergence. This is a different effect than normalizing by ∑ pij as we are not affecting the attraction/repulsion ratio in Theorem 1.

7 We provide formalization of these metrics in the supplementary material (A.8)

8 As such, we employ the Laplacian Eigenmap initialization on small datasets (<100K) due to its predictable output and the random initialization on large datasets (>100K) to avoid slowdowns.

The analysis is slightly more involved in the case of UMAP. Recall that the UMAP algorithm approximates the pij and 1 - pik gradient scalars by sampling the attractions and repulsions proportionally to pij and 1 - pik, which we referred to as scalar sampling. However, the gradients in the normalized setting (Equation 5) lose the 1 - pik scalar on repulsions. The UMAP optimization schema, then, imposes an unnecessary weight on the repulsions in the normalized setting as the repulsions are still sampled according to the nolonger-necessary 1 - pik scalar. Accounting for this requires dividing the repulsive forces by 1 - pik, but this (with the momentum gradient descent and stronger learning rate) leads to a highly unstable training regime. We refer the reader to Figure 7 in the supplementary material for details.

This implies that stabilizing UMAP in the normalized setting requires removing the sampling and instead directly multiplying by pij and 1 - pik. Indeed, this is exactly what we do in GDR. Under this change, GDRumap and GDRtsne obtain effectively identical embeddings to the default UMAP and tSNE ones. This is confirmed in the kNN accuracy and K-means V-score metrics in Table 3.

5.3 Time Efficiency

We lastly discuss the speeds of UMAP, tSNE, GDR, and our accelerated version of GDR in section A.1 of the supplementary material due to space concerns. Our implementations of UMAP and GDR perform gradient descent an order of magnitude faster than the standard UMAP library, implying a corresponding speedup over tSNE. We also provide an acceleration by doing GDR with scalar sampling that provides a further 2× speedup. Despite the fact that this imposes a slight modification onto the effective gradients, we show that this is qualitatively insignificant in the resulting embeddings.

6 Conclusion &FutureWork

We discussed the set of differences between tSNE and UMAP and identified that only the normalization significantly impacts the outputs. This provides a clear unification of tSNE and UMAP that is both theoretically simple and easy to implement. Beyond this, our analysis has uncovered multiple misunderstandings regarding UMAP and tSNE while hopefully also clarifying how these methods work.

We raised several questions regarding the theory of gradient-based DR algorithms. Is there a setting in which the UMAP pseudo-distance changes the embeddings? Does the KL divergence induce a better optimization criterium than the Frobenius norm? Is it true that UMAP's framework can accommodate tSNE's normalization? We hope that we have facilitated future research into the essence of these algorithms through identifying all of their algorithmic components and consolidating them in a simple-to-use codebase.