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 '17 at 9:25
  • $\begingroup$ @swish how come? That is the formula provided... $\endgroup$ – SumNeuron Feb 17 '17 at 10:01
  • $\begingroup$ Variance is already a $\sigma^2$ afaik $\endgroup$ – swish Feb 17 '17 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$ – Alexey Golyshev Feb 17 '17 at 11:12
  • $\begingroup$ @SumNeuron gist.githubusercontent.com/sw1sh/… $\endgroup$ – swish Feb 17 '17 at 11:50

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 '17 at 20:48
  • $\begingroup$ @SumNeuron Good question. But I don't know. $\endgroup$ – Alexey Golyshev Mar 11 '17 at 5:10
  • 2
    $\begingroup$ @sumneuron You can set the perplexity like this Method -> {"TSNE", "Perplexity" -> 100} $\endgroup$ – user5601 Mar 16 '17 at 21:49
  • $\begingroup$ @user5601 thanks! Just f.y.i. the original publication Perplexity was best set between 5-50 :) $\endgroup$ – SumNeuron Mar 17 '17 at 8:26
  • 6
    $\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 '17 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.