Six months ago @M.R. asked about an implementation of the t-Distributed Stochastic Neighbor Embedding (t-SNE) algorithm by van der Maaten and Hinton (2008). (@M.R.'s question)

@Alexey Golyshev gave a solid answer utilizing RLink. However, I thought it would be more interesting* to try and implement t-SNE in Mathematica (v11+). However I have run into a bit of difficulty and would appreciate the community's help.

*it may also have a minimal relation to MacOS related issues with RLink.

Note: I like to code using Mathematica's symbolic notation (and hence it looks like crude even with the Mathematica Stack Exchange Plugin). A notebook with the code can be found here (rather than a mess of a \[VariousLongSymbolNames] for a post).


Page 2587 gives a the pseudo-code for the t-SNE algorithm:

Algorithm 1

Data: data set $\chi=\{x_1, x_2, \dots, x_n\}$,
cost function parameters: perplexity Perp,
optimization parameters: number of iterations $T$, learning rate $\eta$, momentum $\alpha(t)$.
Result: low-dimensional data representation $Y^{(T)}=\{y_1, y_2, \dotsm y_n\}.$

----compute pariwase affinities $p_{j|i}$ with perplexity Perp (using Equation 1)
----set $p_{ij}=\frac{p_{j|i}+p_{i|j}}{2n}$
----sample initial solution $Y^{(0)}=\{y_1,y_2,\dots,y_n\}$ from $N(0,10^{-4}I)$
----for $t=1$ to $T$ do
--------computer low-dimensional affinities $q_{ij}$ (using Equation 4)
--------compute gradient $\frac{\delta C}{\delta Y}$ (using Equation 5)
--------set $Y^{(t)}=Y^{(t-1)}+\eta \frac{\delta C}{\delta Y} + \alpha(t)(Y^{(t-1)}-Y^{(t-2)})$


1 $$ p_{j\vert i}=\frac{\exp(-\| x_i-x_j\|^2/2\sigma^2)}{\sum_{k\neq i}\exp(-\| x_i-x_k\|^2/2\sigma^2)} $$

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

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


So here is what I have been able to implement so far... (see code)


  • how to better handle functions requiring access to elements of the Dataset. I like assigning $x_i$ values inside Table (which does work), however I fear that this is against convention?
  • how can I write this following Mathematica conventions?
  • how to update the learning term?
  • any obvious mistakes?




One may notice that the pseudo-code from the original paper has a typo. It has a recurrence relation for the gradient, which depends on the previous two gradients; only one gradient is initialized, however. Ergo one needs to define $T^{-1}$ to be something (perhaps) all zeros.

In addition, as pointed out in the comments, although variance is continually referred to as $\sigma_i$ in the paper, that is a repeated typo and should be $\sigma^2_i$.

  • $\begingroup$ I is an IdentityMatrix. Also you don't need to square Variance and there exists SquaredEuclideanDistance. $\endgroup$
    – swish
    Feb 17, 2017 at 9:25
  • $\begingroup$ @swish how come? That is the formula provided... $\endgroup$
    – SumNeuron
    Feb 17, 2017 at 10:01
  • $\begingroup$ Variance is already a $\sigma^2$ afaik $\endgroup$
    – swish
    Feb 17, 2017 at 11:11
  • $\begingroup$ @SumNeuron You can take the Python code here and make the naive translation of it to the Wolfram language + Compile. $\endgroup$ Feb 17, 2017 at 11:12
  • $\begingroup$ @SumNeuron gist.githubusercontent.com/sw1sh/… $\endgroup$
    – swish
    Feb 17, 2017 at 11:50

1 Answer 1


In Mathematica 11.1 DimensionReduce supports t-SNE.

enter image description here

  • 3
    $\begingroup$ Which implementation? The Barnes Hut? Or the traditional? Can one set the Perplexity? This would be important given that Perplexity basically defined the scope on which t-SNE searches for local structure in high dimensions. $\endgroup$
    – SumNeuron
    Mar 10, 2017 at 20:48
  • $\begingroup$ @SumNeuron Good question. But I don't know. $\endgroup$ Mar 11, 2017 at 5:10
  • 2
    $\begingroup$ @sumneuron You can set the perplexity like this Method -> {"TSNE", "Perplexity" -> 100} $\endgroup$
    – user5601
    Mar 16, 2017 at 21:49
  • $\begingroup$ @user5601 thanks! Just f.y.i. the original publication Perplexity was best set between 5-50 :) $\endgroup$
    – SumNeuron
    Mar 17, 2017 at 8:26
  • 7
    $\begingroup$ @SumNeuron: we use the original exact-gradient method for very small datasets, and a parallelized C++ implementation of the Barnes-Hut approximation method for larger datasets. $\endgroup$
    – Sebastian
    Mar 27, 2017 at 20:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.