I have change the header to [WIP] (work in progress). Please don't forget to change it to [MRG] (meaning ready for merge review) to attract the attention of more reviewers once your are done with your to-do list.

The narrative documentation will also need an update: http://scikit-learn.org/dev/modules/manifold.html#t-distributed-stochastic-neighbor-embedding-t-sne . The source file is in doc/modules/manifold.rst. I will add a new todo item.

The openmp flag cannot be enabled by default as it will not be build-able with clang (which is the default compiler under OSX) and probably not under old versions of Microsoft Visual C++ required to build Python extensions for Python 2.7 under Windows. So if we really want to use openmp we should at least have a something in setup.py to check that the compiler supports it and turn it of if it is not supported.

Also the openmp implementation in GCC (libgomp) will crash whenever it's use in conjunction with Python multiprocessing (that does fork without exec under posix by default). There is a patch to make libgomp robust to forks but it has not been reviewed yet. So unless openmp brings a significant speedup vs the serial execution, I would just remove parallel support altogether. If we decide to keep it, it should be turned off by default and only enabled with a specific parameters (such as n_jobs commonly used in the rest of the library, despite the bad name).

Do the tests pass on your box? On travis there is a bunch of direct errors (e.g. missing attributes) that make me thing that you have not pushed all the code your run locally (maybe the generated C files are not up to date?).

I would also rename bhtsne.pyx to _barnes_hut_tsne.py for discoverability.

Could you add some test on some toy data that compares the perplexity of the 2 methods and check that they are approximately similar?

Hey Christopher,

Thanks a lot for contributing this. I haven't had time to review it, but
it is exciting! I think that the Barnes-Hut t-SNE algorithm is very
important to us to scale up.

Thanks @GaelVaroquaux and @ogrisel for your comments! @nicktrav and I are slowly working through our TODO list and all of @ogrisel's wonderful sugestions (thank you!) -- we're hoping to be ready for a fuller review in the coming week or two.

Hi @GaelVaroquaux, @ogrisel & @AlexanderFabisch thanks for the comments! 11 hours ago I thought this was going to take a week or two, but I believe we're now ready for review :) I'll try to respond to everyone's comments all at once.

• We've removed all OpenMP requirements -- the overhead of launching threads was bigger than the tradeoff, even for large datasets with 10k elements. Once we've figured out building scikit-learn with OpenMP we can revisit this question, but let's table the issue for the time being.
• We've modified the narrative docs to reflect the new performance and a extra parameters. The change is minimal since the algorithm performs more or less the same, just faster.
• Added toy data testing & comparing the perplexity of the 2 methods.

The rest of the changes we've immediately incorporated. Tests should now pass too, which I think that means this is ready for merge review, so I've changed the title to [MRG].

cc @nicktrav

Could you add the performance comparison to examples/manifold and make sure that it does not take much longer than 10^1 seconds (n_samples=1000-2000)?

I can reproduce your compile error with gcc 4.9.2-10ubuntu13 from ubuntu 15.04. I looks like a Cython bug: when calling printf("some string\n") in Cython, it creates a local char* typed variable that is not a const and that causes the -Werror=security-format to cause a compiler error with gcc 4.9.

When calling printf("some string with a formatted value: %d\n", i) then the generated C code has a const char* as first argument and everything is fine.

The following hack seems to work:

diff --git a/sklearn/manifold/_barnes_hut_tsne.pyx b/sklearn/manifold/_barnes_hut_tsne.pyx
index a5579b5..3978c80 100644
--- a/sklearn/manifold/_barnes_hut_tsne.pyx
+++ b/sklearn/manifold/_barnes_hut_tsne.pyx
@@ -14,6 +14,8 @@ from libc.math cimport sqrt, log
 cimport numpy as np
 import numpy as np

+cdef char* EMPTY_STRING = ""
+
 cdef extern from "math.h":
     float fabsf(float x) nogil

@@ -751,13 +753,19 @@ def gradient(float[:,:] pij_input,
     err = insert_many(qt, pos_output)
     assert err == 0, "[t-SNE] Insertion failed"
     if verbose > 10:
-        printf("[t-SNE] Computing gradient\n")
+        # XXX: format hack to workaround lack of `const char *` type
+        # in the generated C code that triggers error with gcc 4.9
+        # and -Werror=format-security
+        printf("[t-SNE] Computing gradient\n%s", EMPTY_STRING)
     sum_Q = compute_gradient(pij_input, pos_output, neighbors, forces,
                              qt.root_node, theta, dof, skip_num_points, -1)
     C = compute_error(pij_input, pos_output, neighbors, sum_Q, n_dimensions,
                       verbose)
     if verbose > 10:
-        printf("[t-SNE] Checking tree consistency \n")
+        # XXX: format hack to workaround lack of `const char *` type
+        # in the generated C code
+        # and -Werror=format-security
+        printf("[t-SNE] Checking tree consistency\n%s", EMPTY_STRING)
     cdef long count = count_points(qt.root_node, 0)
     m = ("Tree consistency failed: unexpected number of points=%i "
          "at root node=%i" % (count, qt.root_node.cumulative_size))

There are other occurrences of the same issue in the file. The above patch does not fix them all.

Alas, I was overly optimistic in how much time I'd have to work on this. Can we make a fresh TODO list of what this PR needs?

@cemoody welcome to the club ;)
Maybe reset it to #4887 and fix the errors.

@ogrisel is the cython bug fixed in newer cython? Should we report?

I think having a sparse matrix variant would be nice but I would be happy to merge without it.

I added @ogrisel's fix to #4887. However, the transform method can't work the way it is implemented now. I haven't really looked at it, but it currently assumes that the input dimension is the same as the output dimensions somehow. It concatenates the embedding with the new data. I think it should just concatenate the training data with the new data, right? But then the tests should be tighter.

Fun api issue: calling transform will change _X_fit, so calling transform changes again afterwards will change the results...

Also, calling transform changes the embedding of the training data. I thought that was supposed to be avoided using skip_num_points.
transform also returns the transform of the training data plus the additional data, which breaks api.
Maybe we should try to merge this without the transform? Or does someone have time to look into it?

Is this project still in progress? I'd really love to be able to use this new algorithm.

@zachmayer the algorithm is in master and will be in the imminent release. It uses dense data structures, though, so the memory requirements will be the same as before. There is no-one working on using sparse data structures atm afaik.

I can live with that! Awesome!

Oh, I didn't realize it was in master! That's super awesome & I'm incredibly excited!!! Thanks @amueller for all of your hard work :)

@amueller Can you link to the PR/commit that puts it in master? Thanks!

@zachmayer I think it's at this commit: fe09a07

This one: #4887

I left your PR open, as I removed the "transform" method. It changed the training data, which is not really allowed. I think that might not have been what you intended, was it?

I assume this should be closed?

It can be reopened if that's appropriate.

@cemoody @amueller is there still ongoing work to optimize memory consumption occurring in calculating the pairwise distances? I'm getting a MemoryError on a 75,000x50 dataset.

@jamiis AFAIK nobody is working on it at the moment. Do you want to give it a try?

In principle, shouldn't there be a way to avoid calculating the full distance matrix? The original author of t-SNE says "Barnes-Hut algorithms are intended to operate on large datasets for which you cannot keep such distance matrices in memory" (see lvdmaaten/bhtsne#7), so I'm wondering if the current approach in the sklearn version (starting by calculating the full distance matrix) is inherently flawed. Though honestly I'm not sure what the correct alternative is....

I think the best way to implement it is a sparse distance matrix where only the relevant entries are included. However, it is a little bit complicated to handle all special cases correctly, e.g. precomputed distance matrices.

Any progress? I'm unable to find an implementation that lets you both work on large datasets and vary the number of iterations.

Moving the discussion to #7089. PR welcome ;)