Barnes-Hut-SNE
Authors: Laurens van der Maaten
ICLR 2013

Abstract: The paper presents an O(N log N)-implementation of t-SNE -- an embedding technique that is commonly used for the visualization of high-dimensional data in scatter plots and that normally runs in O(N^2). The new implementation uses vantage-point trees to compute sparse pairwise similarities between the input data objects, and it uses a variant of the Barnes-Hut algorithm - an algorithm used by astronomers to perform N-body simulations - to approximate the forces between the corresponding points in the embedding. Our experiments show that the new algorithm, called Barnes-Hut-SNE, leads to substantial computational advantages over standard t-SNE, and that it makes it possible to learn embeddings of data sets with millions of objects.

Barnes-Hut-SNE
Laurens van der Maaten
Pattern Recognition and Bioinformatics Group, Delft University of Technology
Mekelweg 4, 2628 CD Delft, The Netherlands
lvdmaaten@gmail.com
Abstract
The paper presents an O(Nlog N)-implementation of t-SNE — an embedding
technique that is commonly used for the visualization of high-dimensional data
in scatter plots and that normally runs in O(N2). The new implementation uses
vantage-point trees to compute sparse pairwise similarities between the input data
objects, and it uses a variant of the Barnes-Hut algorithm to approximate the forces
between the corresponding points in the embedding. Our experiments show that
the new algorithm, called Barnes-Hut-SNE, leads to substantial computational ad-
vantages over standard t-SNE, and that it makes it possible to learn embeddings
of data sets with millions of objects.
1 Introduction
Data-visualization techniques are an essential tool for any data analyst, as they allow the analyst to
visually explore the data and generate hypotheses. One of the key limitations of traditional visual-
ization techniques such as histograms, scatter plots, and parallel coordinate plots (see, e.g., [10] for
an overview) is that they only facilitate the visualization of one or a few data variables at a time.
In order to get an idea of the structure of all variables in the data, it is therefore necessary to per-
form an automatic analysis of the data before making the visualization, for instance, by learning a
low-dimensional embedding of the data. In such an embedding, each data object is represented by a
low-dimensional point in such a way, that nearby points correspond to similar data objects and that
distant points correspond to dissimilar data objects. The low-dimensional embedding can readily be
visualized in, e.g., a scatter plot or a parallel coordinate plot.
A plethora of embedding techniques have been proposed over the last decade, e.g., [5, 15, 20, 23,
25, 26]. For creating two- or three-dimensional embeddings that can be readily visualized in a
scatter plot, a family of techniques based onstochastic neighbor embedding(SNE; [11]) has recently
become very popular. These techniques compute anN×Nsimilarity matrix in both the original data
space and in the low-dimensional embedding space; the similarities take the form of a probability
distribution over pairs of points in which high probabilities correspond to similar objects or points.
The probabilities are generally deﬁned as normalized Gaussian or Student-t kernels, which makes
that SNE focuses on preserving local data structure. The embedding is learned by minimizing the
Kullback-Leibler divergence between the probability distributions in the original data space and the
embedding space with respect to the locations of the points in the embedding. As the resulting
cost function is non-convex, this minimization is typically performed using ﬁrst- or second-order
gradient-descent techniques [5, 11, 27]. The gradient of the Kullback-Leibler divergence may be
interpreted as an N-body system in which all of the N points exert forces on each other.
One of the key limitations of SNE (and of its variants) is that its computational and memory com-
plexity scales quadratically in the number of data objectsN. In practice, this limits the applicability
of SNE to data sets with only a few thousand points. To visualize larger data sets, landmark imple-
mentations of SNE may be used [25], but this is hardly a satisfactory solution.
1
arXiv:1301.3342v2  [cs.LG]  8 Mar 2013
In this paper, we develop a new algorithm for (t-)SNE that requires only O(Nlog N) computation
and O(N) memory. Our new algorithm computes a sparse approximation of the similarities between
the original data objects using vantage-point trees [31], and subsequently, approximates the forces
between the points in the embedding using a Barnes-Hut algorithm [1] — an algorithm commonly
used by astronomers to performN-body simulations. The Barnes-Hut algorithm reduces the number
of pairwise forces that needs to be computed by exploiting the fact that the forces exerted by a group
of points on a point that is relatively far away are all very similar.
2 Related work
A large body of previous work has focused on decreasing the computational complexity of algo-
rithms that scale quadratically in the amount of data when implemented naively. Most of these
studies focus on speeding up nearest-neighbor searches using space-partitioning (metric) trees (e.g.,
B-trees [2], cover trees [3], and vantage-point trees [31]) or using locality sensitive hashing ap-
proaches (e.g., [12, 29]). Motivated by their strong performance reported in earlier work in [17], we
opt to use metric trees to approximate the similarities of the input objects in our algorithm.
Several prior studies have also developed algorithms to speed up N-body computations. Most
prominently, [7, 8] developed a dual-tree algorithm that is similar in spirit to the Barnes-Hut al-
gorithm we use in this work. The dual-tree algorithm does not consider interactions between single
points and groups of points like the Barnes-Hut algorithm, but it only considers interactions between
groups of points. In preliminary experiments (see appendix), we found the dual-tree and Barnes-
Hut algorithms to perform on par when used in the context of t-SNE — we opt for the Barnes-Hut
algorithm here because it is conceptually simpler. Prior work [6] has also used the fast Gaussian
transform [9, 30] (a special case of a fast multipole method [19]) to speed up the computation of
Gaussian N-body interactions. Since in t-SNE, the forces exerted on the bodies are non-Gaussian,
such an approach cannot readily be applied here.
3 t-Distributed Stochastic Neighbor Embedding
t-Distributed Stochastic Neighbor Embedding (t-SNE) minimizes the divergence between two dis-
tributions: a distribution that measures pairwise similarities between the original data objects and a
distribution that measures pairwise similarities between the corresponding points in the embedding.
Suppose we are given a data set of objects D= {x1,x2,..., xN}and a function d(xi,xj) that
computes a distance between a pair of objects, e.g., their Euclidean distance. Our aim is to learn
an s-dimensional embedding in which each object is represented by a point, E= {y1,y2,..., yN}
with yi ∈Rs. To this end, t-SNE deﬁnes joint probabilities pij that measure the pairwise similarity
between objects xi and xj by symmetrizing two conditional probabilities as follows:
pj|i = exp(−d(xi,xj)2/2σ2
i)∑
k̸=iexp(−d(xi,xk)2/2σ2
i), p i|i = 0 (1)
pij = pj|i + pi|j
2N . (2)
Herein, the bandwidth of the Gaussian kernels σi is set such that the perplexity of the conditional
distribution Pi equals a predeﬁned perplexity u. The optimal value of σi varies per object, and
is found using a simple binary search; see [11] for details. A heavy-tailed distribution is used to
measure the similarity qij between the two corresponding points yi and yj in the embedding:
qij = (1 +∥yi −yj∥2)−1
∑
k̸=l(1 +∥yk −yl)∥2)−1 , q ii = 0. (3)
In the embedding, a normalized Student-t kernel is used to measure similarities rather than a nor-
malized Gaussian kernel to account for the difference in volume between high- and low-dimensional
spaces [25]. The locations of the embedding points yi are learned by minimizing the Kullback-
Leibler divergence between the joint distributions P and Q:
C(E) =KL(P||Q) =
∑
i̸=j
pijlog pij
qij
. (4)
2
This cost function is non-convex; it is typically minimized by descending along the gradient:
∂C
∂yi
= 4
∑
j̸=i
(pij −qij)qijZ(yi −yj), (5)
where we deﬁned the normalization term Z = ∑
k̸=l(1 +∥yk −yl)∥2)−1. The evaluation of both
joint distributions P and Q is O(N2), because their respective normalization terms sum over all
N2 pairs of points. Since t-SNE scales quadratically in the number of objects N, its applicability is
limited to data sets with only a few thousand data objects; beyond that, learning becomes very slow.
4 Barnes-Hut-SNE
Barnes-Hut-SNE uses metric trees to approximate P by a sparse distribution in which only O(uN)
values are non-zero, and approximates the gradients ∂C
∂yi
using a Barnes-Hut algorithm.
4.1 Approximating Input Similarities
As the input similarities are computed using a (normalized) Gaussian kernel, probabilitiespij corre-
sponding to dissimilar input objectsiand jare (nearly) inﬁnitesimal. Therefore, we can use a sparse
approximation to the probabilities pij without a substantial negative effect on the quality of the ﬁ-
nal embeddings. In particular, we compute the sparse approximation by ﬁnding the ⌊3u⌋nearest
neighbors of each of the N data objects, and redeﬁning the pairwise similarities pij as:
pj|i =
{
exp(−d(xi,xj)2/2σ2
i )∑
k∈Ni exp(−d(xi,xk)2/2σ2
i ) , if j ∈Ni
0, otherwise
(6)
pij = pj|i + pi|j
2N . (7)
Herein, Nirepresents the set of the⌊3u⌋nearest neighbors ofxi, and σiis set such that the perplexity
of the conditional distribution equals u. The nearest neighbor sets are found in O(uNlog N) time
by building a vantage-point tree on the data set.
Vantage-point tree. In a vantage-point tree, each node stores a data object and the radius of a
(hyper)ball that is centered on this object [31]. All non-leaf nodes have two children: data objects
that are located inside the ball are stored under the left child of the node, whereas data objects that
are located outside the ball are stored under the right child. The tree is constructed by presenting
the data objects one-by-one, traversing the tree based on whether the current data object lies inside
or outside a ball, and creating a new leaf node in which the object is stored. The radius of the new
leaf node is set to the median distance between its object and all other objects that lie inside the ball
represented by its parent node. To construct a vantage-point tree, the objects need not necessarily
be points in a high-dimensional feature space; the availability of a metric d(xi,xj) sufﬁces. (In our
experiments, however, we use xi ∈RD and d(xi,xj) =∥xi −xj∥.)
A nearest-neighbor search is performed using a depth-ﬁrst search on the tree that computes the dis-
tance of the objects stored in the nodes to the target object, whilst maintaining i) a list of the current
nearest neighbors and ii) the distance τ to the furthest nearest neighbor in the current neighbor list.
The value of τ determines whether or not a node should be explored: if there can still be objects
inside the ball whose distance to the target object is smaller than τ, the left node is searched, and if
there can still be objects outside the ball whose distance to the target object is smaller than τ, the
right node is searched. The order in which children are searched depends on whether or not the tar-
get object lies inside or outside the current node ball: the left child is examined ﬁrst if the object lies
inside the ball, because the odds are that the nearest neighbors of the target object are also located
inside the ball. The right child is searched ﬁrst whenever the target object lies outside of the ball.
4.2 Approximating t-SNE Gradients
To approximate the t-SNE gradient, we start by splitting the gradient into two parts as follows:
∂C
∂yi
= 4(Fattr −Frep) = 4

∑
j̸=i
pijqijZ(yi −yj) −
∑
j̸=i
q2
ijZ(yi −yj)

, (8)
3
where Fattr denotes the sum of all attractive forces (the left sum), whereas Frep denotes the sum of
all repulsive forces (the right sum). Computing the sum of all attractive forces, Fattr, is computa-
tionally efﬁcient; it can be done by summing over all non-zero elements of the sparse distributionP
in O(uN). (Note that the term qijZ = (1 +∥yi −yj∥2)−1 can be computed in O(1).) However, a
naive computation of the sum of all repulsive forces,Frep, is O(N2). We now develop a Barnes-Hut
algorithm to approximate Frep efﬁciently in O(Nlog N).
Consider three points yi, yj, and yk with ∥yi−yj∥≈∥yi−yk∥≫∥yj−yk∥. In this situation, the
contributions of yj and yk to Frep will be roughly equal. The Barnes-Hut algorithm [1] exploits this
by i) constructing a quadtree on the current embedding, ii) traversing the quadtree using a depth-ﬁrst
search, and iii) at every node in the quadtree, deciding whether the corresponding cell can be used
as a “summary” for the gradient contributions of all points in that cell.
Figure 1: Quadtree constructed on a two-
dimensional t-SNE embedding of 500 MNIST
digits (the colors of the points correspond to the
digit classes). Note how the quadtree adapts to
the local point density in the embedding.
Quadtree. A quadtree is a tree in which each
node represents a rectangular cell with a partic-
ular center, width, and height. Non-leaf nodes
have four children that split up the cell into four
smaller cells (quadrants) that lie “northwest”,
“northeast”, “southwest”, and “southeast” of
the center of the parent node (see Figure 1 for
an illustration). Leaf nodes represent cells that
contain at most one point of the embedding; the
root node represents the cell that contains the
complete embedding. In each node, we store
the center-of-mass of the embedding points that
are located inside the corresponding cell, ycell,
and the total number of points that lie inside the
cell, Ncell. A quadtree hasO(N) nodes and can
be constructed in O(N) time by inserting the
points one-by-one, splitting a leaf node when-
ever a second point is inserted in its cell, and
updating ycell and Ncell of all visited nodes.
Approximating the gradient.To approximate
the repulsive part of the gradient, Frep, we
note that if a cell is sufﬁciently small and suf-
ﬁciently far away from point yi, the contri-
butions q2
ijZ(yi −yj) to Frep will be roughly similar for all points yj inside that cell. We
can, therefore, approximate these contributions by Ncellq2
i,cellZ(yi −ycell), where we deﬁne
qi,cellZ = (1 +∥yi −ycell∥2)−1. We ﬁrst approximate FrepZ = q2
ijZ2(yi −yj) by perform-
ing a depth-ﬁrst search on the quadtree, assessing at each node whether or not that node may be
used as a “summary” for all the embedding points that are located in the corresponding cell. During
this search, we construct an estimate of Z = ∑
i̸=j(1 +∥yi −yj∥2)−1 in the same way. The two
approximations thus obtained are then used to compute Frep via Frep = FrepZ
Z .
We use the condition proposed by [1] to decide whether a cell may be used as a “summary” for all
points in that cell. The condition compares the distance of the cell to the target point with its size:
∥yi −ycell∥2/rcell <θ, (9)
where rcell represents the length of the diagonal of the cell under consideration and θis a threshold
that trades off speed and accuracy (higher values of θ lead to poorer approximations). In prelimi-
nary experiments, we also explored various other conditions that take into account the rapid decay
of the Student-t tail, but we did not ﬁnd to lead these alternative conditions to lead to a better
accuracy-speed trade-off. (The problem of more complex conditions is that they require expensive
computations at each cell. By contrast, the condition in Equation 9 can be evaluated very rapidly.)
Dual-tree algorithms. Whilst the Barnes-Hut algorithm considers point-cell interactions, further
speed-ups may be obtained by computing only cell-cell interactions. This can be done using a
dual-tree algorithm [7] that simultaneously traverses the quadtree twice, and for every pair of nodes
decides whether the interaction between the corresponding cells can be used as “summary” for the
interactions between all points inside these two cells. Perhaps surprisingly, we did ﬁnd such an
4
Computation time Nearest neighbor error 
Figure 2: Computation time (in seconds) required to embed70,000 MNIST digits using Barnes-Hut-
SNE (left) and the 1-nearest neighbor errors of the corresponding embeddings (right) as a function
of the trade-off parameter θ.
approach to perform on par with the Barnes-Hut algorithm in preliminary experiments. The com-
putational advantages of the dual-tree algorithm evaporate because after computing an interaction
between two cells, one still needs to determine to which set of points the interaction applies. This can
be done by searching the cell or by storing a list of children in each node during tree construction.
Both these approaches are computationally costly. (It should be noted that the dual-tree algorithm
is, however, much faster in approximating the value of the t-SNE cost function.) The results of our
experiments with dual-tree algorithms are presented in the appendix.
5 Experiments
We performed experiments on four large data sets to evaluate the performance of Barnes-Hut-SNE.
Code for our algorithm is available from http://homepage.tudelft.nl/19j49/tsne.
Data sets.We performed experiments on four data sets: i) the MNIST data set, ii) the CIFAR-10 data
set, iii) the NORB data set, and iv) the TIMIT data set. The MNIST data set contains N= 70,000
grayscale handwritten digit images of size D= 28×28 = 784pixels, each of which corresponds
to one of ten classes. The CIFAR-10 data set [14] is an annotated subset of the 80 million tiny
images data set [24] that contains N = 70,000 RGB images of size 32×32 pixels, leading to a
D=32×32×3=3 ,072-dimensional input objects; each image corresponds to one of ten classes. The
(small) NORB data set [16] contains grayscale images of toys from ﬁve different classes, rendered on
a uniform background under 6 lighting conditions, 9 elevations (30 to 70 degrees every 5 degrees),
and 18 azimuths (0 to 340 every 20 degrees). All images contain D = 96×96 = 9,216 pixels.
The TIMIT data set contains speech data from which MFCC, delta, and delta-delta features were
extracted, leading to D= 39-dimensional features [22]; each frame in the data has one of 39 phone
labels. We used the TIMIT training set of N=1,105,455 frames in our experiments.
Experimental setup. In all our experiments, we follow the experimental setup of [25] as closely
as possible. In particular, we initialize the embedding points by sampling from a Gaussian with a
variance of 10−4. We run a gradient-descent optimizer for 1,000 iterations, setting the initial step
size to 200; the step size is updated during the optimization use the scheme of [13]. We use an
additional momentum term that has weight 0.5 during the ﬁrst 250 iterations, and 0.8 afterwards.
The perplexity uis ﬁxed to 30. Following [25], all data sets with a dimensionality Dlarger than 50
were preprocessed using PCA to reduce their dimensionality to 50.
During the ﬁrst 250 learning iterations, we multiplied allpij-values by a user-deﬁned constantα>1.
As explained in [25], this trick enables t-SNE to ﬁnd a better global structure in the early stages of the
optimization. In preliminary experiments, we found that this trick becomes increasingly important to
obtain good embeddings when the data set size increases, as it becomes harder for the optimization
5
Computation time Nearest neighbor error 
Figure 3: Compution time (in seconds) required to embed MNIST digits (left) and the 1-nearest
neighbor errors of the corresponding embeddings (right) as a function of data set size N for both
standard t-SNE and Barnes-Hut-SNE. Note that the required computation time, which is shown on
the y-axis of the left ﬁgure, is plotted on a logarithmic scale.
to ﬁnd a good global structure when there are more points in the embedding because there is less
space for clusters to move around. In our experiments, we ﬁx α=12 (by contrast, [25] used α=4).
We present the results of three sets of experiments. In the ﬁrst experiment, we investigate the effect
of the trade-off parameter θon the speed and the quality of embeddings produced by Barnes-Hut-
SNE on the MNIST data set. In the second experiment, we investigate the computation time required
to run Barnes-Hut-SNE as a function of the number of data objectsN (also on the MNIST data set).
In the third experiment, we construct and visualize embeddings of all four data sets.
Results. Figure 2 presents the results of an experiment in which we varied the speed-accuracy
trade-off parameter θused to learn the embedding. The ﬁgure shows the computation time required
to construct embeddings of all 70,000 MNIST digit images, as well as the 1-nearest neighbor error
(computed based on the digit labels) of the corresponding embeddings. The results presented in the
ﬁgure show that the trade-off parameter θmay be increased to a value of approximately 0.5 without
negatively affecting the quality of the embedding. At the same time, increasing the value ofθto 0.5
leads to very substantial improvements in terms of the amount of computation required: the time
required to embed all 70,000 MNIST digits is reduced to just 645 seconds when θ=0.5. (Note that
the special case θ= 0corresponds to standard t-SNE [25]; we did not run an experiment with θ= 0
because standard t-SNE would take days to complete on the full MNIST data set.)
In Figure 3, we compare standard t-SNE and Barnes-Hut-SNE in terms of i) the computation time
required for the embedding of MNIST digit images as a function of the data set size N and ii)
the 1-nearest neighbor errors of the corresponding embeddings. (Note that the y-axis of the left
ﬁgure, which represents the required computation time in seconds, uses a logarithmic scale.) In the
experiments, we ﬁxed the parameterθthat trades off speed and accuracy to0.5. The results presented
in the ﬁgure show that Barnes-Hut-SNE is orders of magnitude faster than standard t-SNE, whilst
the difference in quality of the constructed embeddings (which is measured by the nearest-neighbor
errors) is negligible. Most prominently, the computational advantages of Barnes-Hut-SNE rapidly
increase as the number of objects in the data set N increases.
Figure 4 presents embeddings of all four data sets constructed using Barnes-Hut-SNE. The colors
of the points indicate the classes of the corresponding objects; the titles of the plots indicate the
computation time that was used to construct the corresponding embeddings. As before, we ﬁxedθ=
0.5 in all four experiments. The results in the ﬁgure shows that Barnes-Hut-SNE can construct high-
quality embeddings of, e.g., the 70,000 MNIST handwritten digit images in just over 10 minutes.
(Although our MNIST embedding contains many more points, it may be compared with that in [25].
Visually, the structure of the two embeddings is very similar.) The results also show that Barnes-
6
Figure 4: Barnes-Hut-SNE visualizations of four data sets: MNIST handwritten digits (top-left),
CIFAR-10 tiny images (top-right), NORB object images (bottom-left), and TIMIT speech frames
(bottom-right). The colors of the point indicate the classes of the corresponding objects. The titles of
the ﬁgures indicate the computation time that was used to construct the corresponding embeddings.
Figure best viewed in color.
Hut-SNE makes it practical to embed data sets with more than a million data points: the TIMIT
embedding shows all 1,105,455 data points, and was constructed in less than four hours.
A version of the MNIST embedding in which the original digit images are shown is presented in Fig-
ure 5. The results show that, like standard t-SNE, Barnes-Hut-SNE is very good at preserving local
structure of the data in the embedding: for instance, the visualization clearly shows that orientation
is one of the main sources of variation within the cluster of ones.
6 Conclusion and Future Work
We presented a new t-SNE algorithm [25], called Barnes-Hut-SNE, that i) constructs a sparse ap-
proximation to the similarities between input objects using vantage-point trees, and ii) approxi-
mates the t-SNE gradient using a variant of the Barnes-Hut algorithm. The new algorithm runs in
O(Nlog N) rather than O(N2), and requires only O(N) memory. Our experimental evaluation of
Barnes-Hut-SNE shows that it is substantially faster than standard t-SNE, and that it facilitates the
visualization of data sets with millions of data objects in scatter plots.
A drawback of Barnes-Hut-SNE is that it does not provide any error bounds [21]. Indeed, there
exist alternative algorithms that do provide such error bounds ( e.g., [28]); we aim to explore these
alternatives in future work to see whether they can be used to bound the error made in our t-SNE
7
Figure 5: Barnes-Hut-SNE visualization of all 70,000 MNIST handwritten digit images (con-
structed in 10 minutes and 45 seconds). Zoom in on the visualization for more detailed views.
gradient computations, and to bound the error in the ﬁnal embedding. Another limitation of Barnes-
Hut-SNE is that it can only be used to embed data in two or three dimensions. Generalizations to
higher dimensions are infeasible because the size of the tree grows exponentially in the dimensional-
ity of the embedding space. Having said that, this limitation is not very severe since t-SNE is mainly
used for visualization ( i.e. for embedding in two or three dimensions). Moreover, it is relatively
straightforward to replace the quadtree by metric trees that scale better to high-dimensional spaces.
In future work, we plan to further scale up our algorithm by developing parallelized implementations
that can run on data sets that are too large to be fully stored in memory. We also aim to investigate
the effect of varying the value of θduring the optimization. In addition, we plan to explore to what
extent adapted versions of our algorithm (that use metric trees instead of quadtrees) can be used to
speed up techniques for relational embedding (e.g., [4, 18]).
Acknowledgments
The author is supported by EU-FP7 Social Signal Processing (SSPNet) and by the Netherlands In-
stitue for Advanced Study (NIAS). The author thanks Geoffrey Hinton for many helpful discussions,
and two anonymous reviewers for their helpful comments.
8
References
[1] J. Barnes and P. Hut. A hierarchical O(N log N) force-calculation algorithm. Nature, 324(4):446–449,
1986.
[2] R. Bayer and E. McCreight. Organization and maintenance of large ordered indexes. Acta Informatica,
1(3):173–189, 1972.
[3] A. Beygelzimer, S. Kakade, and J. Langford. Cover trees for nearest neighbor. In Proceedings of the
International Conference on Machine Learning, pages 97–104, 2006.
[4] A. Bordes, J. Weston, R. Collobert, and Y . Bengio. Learning structured embeddings of knowledge bases.
In Proceedings of the 25th Conference on Artiﬁcial Intelligence (AAAI), 2011.
[5] M. ´A. Carreira-Perpi˜n´an. The elastic embedding algorithm for dimensionality reduction. In Proceedings
of the 27th International Conference on Machine Learning, pages 167–174, 2010.
[6] N. DeFreitas, Y . Wang, M. Mahdaviani, and D. Lang. Fast Krylov methods for N-body learning. In
Advances in Neural Information Processing Systems, volume 18, pages 251–258, 2006.
[7] A.G. Gray and A.W. Moore. N-body problems in statistical learning. In Advances in Neural Information
Processing Systems, pages 521–527, 2001.
[8] A.G. Gray and A.W. Moore. Rapid evaluation of multiple density models. In Proceedings of the Interna-
tional Conference on Artiﬁcial Intelligence and Statistics, 2003.
[9] L. Greengard and V . Rokhlin. A fast algorithm for particle simulations.Journal of Computational Physics,
73:325–348, 1987.
[10] J. Heer, M. Bostock, and V . Ogievetsky. A tour through the visualization zoo. Communications of the
ACM, 53:59–67, 2010.
[11] G.E. Hinton and S.T. Roweis. Stochastic Neighbor Embedding. In Advances in Neural Information
Processing Systems, volume 15, pages 833–840, 2003.
[12] P. Indyk and R. Motwani. Approximate nearest neighbors: Towards removing the curse of dimensionality.
In Proceedings of 30th Symposium on Theory of Computing, 1998.
[13] R.A. Jacobs. Increased rates of convergence through learning rate adaptation. Neural Networks, 1:295–
307, 1988.
[14] A. Krizhevsky. Learning multiple layers of features from tiny images. Technical report, University of
Toronto, 2009.
[15] N.D. Lawrence. Spectral dimensionality reduction via maximum entropy. In Proceedings of the Inter-
national Conference on Artiﬁcial Intelligence and Statistics, JMLR W&CP , volume 15, pages 51–59,
2011.
[16] Y . LeCun, F.J. Huang, and L. Bottou. Learning methods for generic object recognition with invariance to
pose and lighting. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition,
pages 97–104, 2004.
[17] T. Liu, A.W. Moore, A. Gray, and K. Yang. An investigation of practical approximate nearest neighbor
algorithms. In Advances in Neural Information Processing Systems, volume 17, pages 825–832, 2004.
[18] A. Paccanaro and G.E. Hinton. Learning distributed representations of concepts using linear relational
embedding. IEEE Transactions on Knowledge and Data Engineering, 13(2):232–244, 2001.
[19] V . Rokhlin. Rapid solution of integral equations of classic potential theory. Journal of Computational
Physics, 60:187–207, 1985.
[20] S.T. Roweis and L.K. Saul. Nonlinear dimensionality reduction by Locally Linear Embedding. Science,
290(5500):2323–2326, 2000.
[21] J.K. Salmon and M.S. Warren. Skeletons from the treecode closet. Journal of Computational Physics ,
111(1):136–155, 1994.
[22] F. Sha and L.K. Saul. Large margin hidden Markov models for automatic speech recognition. InAdvances
in Neural Information Processing Systems, volume 19, pages 1249–1456, 2007.
[23] J.B. Tenenbaum, V . de Silva, and J.C. Langford. A global geometric framework for nonlinear dimension-
ality reduction. Science, 290(5500):2319–2323, 2000.
[24] A. Torralba, R. Fergus, and W.T. Freeman. 80 million tiny images: A large dataset for non-parametric
object and scene recognition. IEEE Transactions on Pattern Analysis and Machine Intelligence ,
30(11):1958–1970, 2008.
[25] L.J.P. van der Maaten and G.E. Hinton. Visualizing data using t-SNE. Journal of Machine Learning
Research, 9(Nov):2431–2456, 2008.
9
[26] J. Venna, J. Peltonen, K. Nybo, H. Aidos, and S. Kaski. Information retrieval perspective to nonlinear
dimensionality reduction for data visualization. Journal of Machine Learning Research , 11(Feb):451–
490, 2010.
[27] M. Vladymyrov and M. ´A. Carreira-Perpi ˜n´an. Partial-Hessian strategies for fast learning of nonlinear
embeddings. In Proceedings of the International Conference on Machine Learning , pages 345–352,
2012.
[28] M.S. Warren and J.K. Salmon. A parallel hashed octtree N-body algorithm. In Proceedings of the
ACM/IEEE Conference on Supercomputing, pages 12–21, 1993.
[29] Y . Weiss, A. Torralba, and R. Fergus. Spectral hashing. In Advances in Neural Information Processing
Systems, pages 1753–1760, 2008.
[30] C. Yang, R. Duraiswami, N.A. Gumerov, and L. Davis. Improved fast Gauss transform and efﬁcient
kernel density estimation. In Proceedings of the IEEE International Conference on Computer Vision ,
pages 664–671, 2003.
[31] P.N. Yianilos. Data structures and algorithms for nearest neighbor search in general metric spaces. In
Proceedings of the ACM-SIAM Symposium on Discrete Algorithms, pages 311–321, 1993.
10
A Experiments with Dual-Tree t-SNE
We also performed experiments with a dual-tree implementation [7] of t-SNE. Dual-tree t-SNE differs from
Barnes-Hut-SNE in that it considers only cell-cell instead of point-cell interactions. It simultaneously traverses
the quadtree twice, and decides for each pair of nodes whether the interaction between these nodes can be used
as a “summary” for all points in the cells corresponding to these two nodes. We use the following condition to
check whether the interaction between a pair of nodes may be used as a “summary” interaction:
∥ycell1 −ycell2∥2/max(rcell1,rcell2) <ρ, (10)
where ycell1 and ycell2 represent the center-of-mass of the two cells,rcell1 and rcell2 represent the diameter of
the two cells, and ρis a speed-accuracy trade-off parameter (similar to θin Barnes-Hut-SNE).
Figure 6 presents the results of an experiment in which we investigate the inﬂuence of the trade-off parameter
ρon the learning time and the quality of the embedding on the MNIST data set. The results in the ﬁgure may
be readily compared to those in Figure 2. The results in the ﬁgure show that, whilst the dual-tree algorithm
provides additional speed-ups compared to the Barnes-Hut algorithm, the quality of the embedding also de-
teriorates much faster as the trade-off parameter ρincreases. The quality of the embedding obtained with a
dual-tree algorithm with ρ= 0.25 roughly equals that of a Barnes-Hut embedding with θ= 0.5, and these two
embeddings are constructed in roughly the same time (viz. in approximately 650–700 seconds). Figure 7 shows
the performance of dual-tree t-SNE with ρ=0.25 as a function of the number of MNIST digits N. The results
in Figure 7 can be readily compared to those in Figure 3. Again, the results show that dual-tree t-SNE performs
roughly on par with Barnes-Hut-SNE, irrespective of the size of the data set N.
Computation time Nearest neighbor error
Figure 6: Computation time (in seconds) required to embed 70,000 MNIST digits using dual-tree
t-SNE (left) and the 1-nearest neighbor errors of the corresponding embeddings (right) as a function
of the trade-off parameter ρ. The results may be compared to those in Figure 2.
Computation time Nearest neighbor error
Figure 7: Compution time (in seconds) required to embed MNIST digits (left) and the 1-nearest
neighbor errors of the corresponding embeddings (right) as a function of data set size N for both
standard t-SNE and dual-tree t-SNE. The results may be compared to those in Figure 3.
11


Peer review:
[Review]:
Thanks a bunch for these insightful reviews and for the useful pointers to related work (some of which I was not aware of)!

In preliminary experiments, I compared locality-sensitive hashing and vantage-point trees in the initial nearest-neighbor (in the high-dimensional space). I found vantage-point trees to perform considerably better, which is why I used them in the final implementation. The strong performance I obtained when using metric trees appears to be in line with the results presented by Liu, Moore, Gray & Yang (2004). I agree with the first reviewer that there are many other (approximate) nearest-neighbor algorithms that could be used here instead. I will clarify this in the paper, and include references to relevant related work.

The work by Nando de Freitas's lab on n-body simulations is very interesting indeed. I don't think it can readily be applied to t-SNE though, as it appears to heavily rely on the (improved) fast Gauss transform, i.e. on the assumption that Gaussian kernels are used. To the best of my knowledge, there is no existing work that uses fast multipole methods to evaluate Student-t kernels (the fast Gauss transform is an example of a fast multipole method), so extension of this work to t-SNE appears non-trivial. It is also unclear whether fast multipole methods would actually outperform Barnes-Hut in practice, because multipole methods tend to have constants that are much worse. Having said that, this is indeed a very interesting direction for future work! I will clarify this in the paper, and make sure to include the relevant references.

I was not aware of the work by Alex Gray's lab on dual-tree algorithms for n-body simulations; indeed, this work seems readily applicable to t-SNE. I'm presently coding up a dual-tree version of my algorithm, and will try to include empirical evaluations with the dual-tree approach in the final version of the paper. I hope to post an updated version of the paper with these results on Arxiv in a week or two.

I agree with the first reviewer that it is interesting to study if the accuracy-speed trade-off can be adapted during the optimization, but I am not sure that I agree that looser bounds should be used in the beginning of the optimization. In fact, the first 100 or so iterations are essential in identifying the global structure of the data --- doing a poor job in those iterations often implies getting stuck in poor local optima. (I guess one can think of it as errors propagating over time in the optimization.) So an optimal strategy may actually be the opposite of what the reviewer suggests: use tight bounds in the early stages of the optimization and looser bounds later on. It's certainly an interesting direction for future work!