MRG: Speed up euclidean_distances with Cython

[vene]

euclidean_distances usually pops up near the top as cumulative time in every profiler output that uses it. (I'm talking about our new toy in http://jenkins-scikit-learn.github.com/scikit-learn-speed/) so optimizing it sounds like a good idea.

Here is my first go. I'm new at writing cython so feedback is highly appreciated. I have a sinking feeling that I'm copying around memory when I don't want that.

Here is a hackish benchmark script: https://gist.github.com/3305769

I have intentionally not included the .c file in this PR, so the commits will be lighter until this is done.

TODO:

• fast CSR dot product (turns out almost impossible?)

P.S. be gentle, I'm typing with a broken finger :)

[mblondel]

@dwf did something similar during the Granada sprint but @GaelVaroquaux and him decided to not merge the code (don't remember why though).

[GaelVaroquaux]

[1]@dwf did something similar during the Granada sprint but
[2]@GaelVaroquaux and him decided to not merge the code (don't remember
why though).

Me neither. I was thinking about this recently and wanted to look at the
code to understand why it wasn't merged, but couldn't find it. @dwf, do
you remember?

[dwf]

Basically it's because in a lot of cases, as @jakevdp pointed out to me, BLAS is faster. Past a certain number of features/number of examples, you are fighting a land war in Asia, so to speak. Thus, better to use np.dot() (or, if you want to call into dgemm/sgemm from Cython, that could be an option too).

One thing that we still need is an out= argument, though, for algorithms where repeated distance computation are the bottleneck, to prevent multiple allocations.

[vene]

On Aug 9, 2012, at 19:13 , David Warde-Farley notifications@github.com wrote:

Basically it's because in a lot of cases, as @jakevdp pointed out to me, BLAS is faster. Past a certain number of features/number of examples, you are fighting a land war in Asia, so to speak. Thus, better to use np.dot() (or, if you want to call into dgemm/sgemm from Cython, that could be an option too).

I am not multiplying arrays by hand, but calling gemm. The speedup comes from only looping once over the array to add the norm of X, the norm of Y, and the clipping. For the Y=None case even more can be gained, I guess. I didn't find any cases where it would be slower. No WOW factor, though.

As for the sparse case, I'm calling safe_sparse_dot outside of cython for now, and passing in the result. The speedup is insignificant, but the norm calculation by itself is faster. (the dot dominates by a factor, of course). Did you also implement the sparse calculation, and could there be cases where scipy.sparse * wins?

And either way, could I please see your code? I want to gulp up all the cython tricks I can.

One thing that we still need is an out= argument, though, for algorithms where repeated distance computation are the bottleneck, to prevent multiple allocations.

Partly implemented. Should I keep going with this, guys?

—
Reply to this email directly or view it on GitHub.

---

Vlad N.
http://vene.ro

[mblondel]

Besides speed up, another benefit is that your implementation avoids intermediary memory allocations (like X * X, for example).

[dwf]

On Thu, Aug 9, 2012 at 1:19 PM, Vlad Niculae notifications@github.comwrote:

I am not multiplying arrays by hand, but calling gemm. The speedup comes
from only looping once over the array to add the norm of X, the norm of Y,
and the clipping. For the Y=None case even more can be gained, I guess. I
didn't find any cases where it would be slower. No WOW factor, though.

Ahh ok. Sorry I didn't look over your code, I'm in the middle of writing my
thesis proposal :S

And either way, could I please see your code? I want to gulp up all the
cython tricks I can.

#483

One thing that we still need is an out= argument, though, for algorithms
where repeated distance computation are the bottleneck, to prevent multiple
allocations.

Partly implemented. Should I keep going with this, guys?

Sounds good to me.

[GaelVaroquaux]

On Thu, Aug 09, 2012 at 10:19:34AM -0700, Vlad Niculae wrote:

I am not multiplying arrays by hand, but calling gemm. The speedup comes
from only looping once over the array to add the norm of X, the norm of Y,
and the clipping. For the Y=None case even more can be gained, I guess.
I didn't find any cases where it would be slower. No WOW factor, though.

I don't have time to look at the code (currently in panic mode), but this
sounds really a good approach, I should say. I am hopefully.

[jakevdp]

Looks good - could you put up some comparison benchmarks?

There's been a lot of buzz on the cython list lately about the best way to pass pointers from numpy to C. (see https://groups.google.com/forum/?fromgroups#!topic/cython-users/8uuxjB_wbBQ[1-25]). It seems the <double*>arr.data pattern may break with newer versions of the numpy C API. According to cython folks, we should be using typed memoryviews for this sort of application, though that requires cython 0.17+ and python 2.5+.

[jakevdp]

Defining arrays like this does not allow fast indexing (this was true in cython <0.16; I'm not sure if it's been fixed that since then).

[jakevdp]

I should specify: I believe defining numpy arrays as a function argument leads to fast indexing, but defining arrays as variables does not.

[dwf]

Really? I've never noticed that, and I'm almost sure I have seen some old generated C that was doing indexing properly on such local variables.

It's definitely the case if you do "cdef np.ndarray foo" that you get Python overhead, but I think if a C type and ndim are there you're alright.

[amueller]

@vene could you compare if indexing via .data is faster?

[vene]

EDIT: Cancel, that, I managed.

The speed is about the same. The generated C code is slightly different. Compare

[vene]

XX[ii] = 0

->

__pyx_t_10 = __pyx_v_ii;
*__Pyx_BufPtrStrided1d(__pyx_t_7sklearn_7metrics_14euclidean_fast_DOUBLE *, __pyx_pybuffernd_XX.rcbuffer->pybuffer.buf, __pyx_t_10, __pyx_pybuffernd_XX.diminfo[0].strides) = 0.0;

versus

(<DOUBLE *>(XX.data))[ii] = 0

->

(((__pyx_t_7sklearn_7metrics_14euclidean_fast_DOUBLE *)__pyx_v_XX->data)[__pyx_v_ii]) = 0.0;

[amueller]

Thanks for giving that a shot. Good to know.

[agramfort]

cython -a gives me some yellow lines that hurts my eyes :)

[dwf]

@jakevdp I recall a discussion in Austin with you and Kurt Smith that there were some memoryview performance issues. Has that been resolved?

[vene]

On Aug 9, 2012, at 23:35 , Alexandre Gramfort notifications@github.com wrote:

cython -a gives me some yellow lines that hurts my eyes :)

I just learned about that moments before submitting the PR, sorry. It's still not obvious to me how to fix some of them.

I can conclude that it's not cool to return an array, but I was gonna use an out= parameter anyway, so that's that. But what's with the RefNannyDeclarations that make all the lines where I declare a function go yellow?

Should I aim for no yellow at all?

—
Reply to this email directly or view it on GitHub.

---

Vlad N.
http://vene.ro

[GaelVaroquaux]

On Thu, Aug 09, 2012 at 03:48:12PM -0700, Vlad Niculae wrote:

cython -a gives me some yellow lines that hurts my eyes :)

I just learned about that moments before submitting the PR, sorry. It's
still not obvious to me how to fix some of them.

I can conclude that it's not cool to return an array,

You are sure? Yellow lines are important only if they occur in a for
loop. If they occur at the beginning or the end of a function, they are
not more (or less costly) than a Python call in the enclosing scope. Now
if you want to call that function many, many times, its a different
matter.

[dwf]

On Thu, Aug 9, 2012 at 6:51 PM, Gael Varoquaux notifications@github.comwrote:

You are sure? Yellow lines are important only if they occur in a for
loop. If they occur at the beginning or the end of a function, they are
not more (or less costly) than a Python call in the enclosing scope. Now
if you want to call that function many, many times, its a different
matter.

Indeed, you should not pay attention to yellow lines at the beginning/end
of a def or cpdef function. You'll see that even if you don't return an
array. It's a bit of a nuisance of the annotator.

I will add to what @GaelVaroquaux said that yellow lines inside a for loop
may indicate that there are other yellow lines you should pay attention to
outside the for loop (i.e. you forgot to cdef something).

[amueller]

I can't believe I didn't know about cython -a ...

[jakevdp]

@dwf - the performance issues come up when passing around many small slices: essentially, the overhead of creating a memoryview slice, though small, can be significant in some circumstances. i.e. if you're slicing an array in order to multiply two 3-dimensional vectors, it's not so good. If you're slicing in order to multiply two 1000-dimensional vectors, then you're probably OK.

The way this PR is written, typed memory views will be as fast or faster than using the cython numpy interface, and will lead to code that's more stable in numpy 1.7+ (with respect to the arr.data problem mentioned above). The only question is whether we want to be dependent on bleeding edge cython (0.16+, or perhaps 0.17), and whether we're OK losing compatibility with python 2.4. Perhaps @GaelVaroquaux can weigh-in on that last bit.

[GaelVaroquaux]

On Fri, Aug 10, 2012 at 06:32:10AM -0700, Jacob Vanderplas wrote:

The only question is whether we want to be dependent on bleeding
edge cython (0.16+, or perhaps 0.17),

That worries me a bit. It forces developpers to instal there own version
of Cython rather than relying on the packaged version for many systems.
That said, we can say that we have a dependency on a recent Cython only
for this file. We would have to be pretty explicit about it in the
comments at the top of the file.

and whether we're OK losing compatibility with python 2.4.

We do not care about Python 2.4.

I see two options:

• either making this specific file depend on a recent
  version of cython. I am worried that the dependency may spread as more
  and more developers want to use fancy cython features.
• or having what works now written as short-lived (and with a comment
  saying so), and what will work in the future as commented out source lines
  of code. This means that we have a 'future deprecation' problem with
  the '.data' access, but I think that we can live with it for a release
  or two.

[vene]

@jakevdp I wrote a cython csr_matmat_t function to compute X * Y.T, imagining it would be faster for this specific setting. My implementation computes every element by finding the common nonzero indices in the rows. I get the good result, but the default, 2-pass implementation kicks my ass. I thought X * Y.T would be a simpler case. Is there any way to improve this or should I keep using scipy's dot?

[jakevdp]

@vene - I think two passes is required for this. Take a look at the scipy implementation:
https://github.com/scipy/scipy/blob/master/scipy/sparse/sparsetools/csr.h#L476
You may be able to improve on this if you're tricky, but the scipy version of the algorithm is very optimized.

[jakevdp]

@GaelVaroquaux good points. I think, though, that because cython is only required for developers and not users, we can safely depend on the latest stable version. If I'm not mistaken, in the past we've settled on everyone using the latest stable version in order to prevent minor changes causing huge diffs in the generated C code.

Also, there is a way to maintain future compatibility without using typed memory views: access the pointer with <double *> np.PyArray_DATA(arr) rather than <double *> arr.data. This has a bit of overhead (it does some reference counting in the background), so it probably shouldn't be used inside loops, but the advantage is it will not break with future numpy versions.

[vene]

On Aug 10, 2012, at 17:48 , Jacob Vanderplas notifications@github.com wrote:

@GaelVaroquaux good points. I think, though, that because cython is only required for developers and not users, we can safely depend on the latest stable version. If I'm not mistaken, in the past we've settled on everyone using the latest stable version in order to prevent minor changes causing huge diffs in the generated C code.

Actually, I toyed with the cd_fast.pyx to see what its cython -a output looks like, and it turns out the C file is generated with an older version than the one I have (pip-installed). So this is not a rule.

---

Vlad N.
http://vene.ro

[mblondel]

You can do the dot product in one pass like this:
https://github.com/mblondel/dot-bench/blob/master/dot.cpp#L200

This assumes that the indices are sorted. Indices can be sorted by csr_matrix.sort_indices(). Internally the object knows if the indices are sorted or not so if the indices are already sorted (which is usually the case), csr_matrix.sort_indices() won't sort the indices again.

[vene]

On Aug 10, 2012, at 18:12 , Mathieu Blondel notifications@github.com wrote:

You can do the dot product in one pass like this:
https://github.com/mblondel/dot-bench/blob/master/dot.cpp#L200

Mathieu, you just blew my mind! :)

[amueller]

This is the same in the svm implementation we are using. Did you get it from there @mblondel?
Btw, did any one ever notice that the ddot.c we are shipping claims to be from '78 ?

[amueller]

For the cython issue: we are using array.data in all the cython implementations, right?

[mblondel]

@amueller Yes that's the same as in libsvm.

[GaelVaroquaux]

Re the fact that it's only needed for devs: I don't buy this argument.
We need to lower the bar for contribution. I believe that easy
contribution is a cornerstone of a dynamic project.

On 8/10/12, Mathieu Blondel notifications@github.com wrote:

@amueller Yes that's the same as in libsvm.

---

Reply to this email directly or view it on GitHub:
#1006 (comment)

[vene]

OK, new observations.
0. I have some yellow np.empty() memory allocations that I'm not particularily fond of.

1. The out= argument cannot be used in the sparse case because there is no (trivial) way to compute XY^T in preallocated memory. @mblondel's implementation should save us here.
2. I made a call to dsyrc, if you think this is too much, we can drop back to dgemm instead. The speedup is not so relevant, we still have to walk the whole output array to fill the lower triangle. Maybe for the symmetrical case, we can do away with half of the matrix, in case all calling points agree? It'd be an API change, though. Maybe uncalled for.
3. TODO: optimize the X sparse - Y dense case, and the precomputed norms case. (this is a note to self).

[GaelVaroquaux]

I don't think that the np.empty is a problem in itself: the yellow
lines correspond to things that would happen anyhow. what could help
would to preallocate memory when calling this function many times.

On 8/10/12, Vlad Niculae notifications@github.com wrote:

OK, new observations.
0. I have some yellow np.empty() memory allocations that I'm not
particularily fond of.

1. The out= argument cannot be used in the sparse case because there is no
   (trivial) way to compute XY^T in preallocated memory. @mblondel's
   implementation should save us here.
2. I made a call to dsyrc, if you think this is too much, we can drop back
   to dgemm instead. The speedup is not so relevant, we still have to walk
   the whole output array to fill the lower triangle. Maybe for the symmetrical
   case, we can do away with half of the matrix, in case all calling points
   agree? It'd be an API change, though. Maybe uncalled for.
3. TODO: optimize the X sparse - Y dense case, and the precomputed norms
   case. (this is a note to self).

---

Reply to this email directly or view it on GitHub:
#1006 (comment)

[vene]

Fixed, the branch was indeed unreachable because csr_matrix always creates a new copy.

How does this look? Should I rename euclidean_fast.pyx to _euclidean_fast.pyx? Should I rebase into a nicer history?

[GaelVaroquaux]

How does this look?

I'll check ASAP.;

Should I rename euclidean_fast.pyx to _euclidean_fast.pyx?

Yes please,

Should I rebase into a nicer history?

If it's not too much of a pain. But only after review, as it will screw
up the github history.

[GaelVaroquaux]

How does this look?

I'll check ASAP.;

This PR looks good. I am 👍 for merge. Before you merge, @vene: did you
fix the seeds in the RNGs of the scikit-learn-speed project? We want to
be able to see the impact of this merge :)

[vene]

This PR looks good. I am 👍 for merge. Before you merge, @vene: did you
fix the seeds in the RNGs of the scikit-learn-speed project? We want to
be able to see the impact of this merge :)

Yes, that is supposed to be fixed now. I hope the impact will be visible.

[GaelVaroquaux]

Regarding the scikit-learn speed page: the readme on github specifies that the URL is at vene.github.com, but I have the impression that it hasn't been updating lately. Is this the right URL?

[vene]

Regarding the scikit-learn speed page: the readme on github specifies that the URL is at vene.github.com, but I have the impression that it hasn't been updating lately. Is this the right URL?

You are correct in that it has moved to http://scikit-learn.github.com/scikit-learn-speed

Thanks for pointing that out. The idea was to link it from the main website and to alias it to speed.scikit-learn.org.

[amueller]

this bit here was refactored in to _build_utils iirc.

[amueller]

Can this be merged? Is another review needed? what is the test coverage?

[GaelVaroquaux]

Can this be merged?

I think that it is good.

what is the test coverage?

Dunno

[amueller]

@vene can you change the setup.py I pointed to and run a coverage report, then I'll merge?

[vene]

Coverage dropped from 90% to 89% in the pairwise module. Line numbers don't match so I have to look at it more closely, I think it should be even higher than 90%. I don't think I can tackle this today though.

[vene]

Coverage report now:

Name                                   Stmts   Miss  Cover   Missing
--------------------------------------------------------------------
sklearn.metrics                           13      0   100%   
sklearn.metrics.cluster                   13      0   100%   
sklearn.metrics.cluster.supervised       151      2    99%   81, 271
sklearn.metrics.cluster.unsupervised      27      1    96%   80
sklearn.metrics.metrics                  249     18    93%   68-70, 285, 290, 608, 655, 677, 720, 842-843, 845, 915, 950-955, 980-981
sklearn.metrics.pairwise                 162     13    92%   94, 96, 270, 275, 416, 520, 544, 583, 662, 671, 681, 686, 689
--------------------------------------------------------------------
TOTAL                                    615     34    94%   
----------------------------------------------------------------------

Coverage report before:

Name                                   Stmts   Miss  Cover   Missing
--------------------------------------------------------------------
sklearn.metrics                           13      0   100%   
sklearn.metrics.cluster                   13      0   100%   
sklearn.metrics.cluster.supervised       138      3    98%   83, 253, 781
sklearn.metrics.cluster.unsupervised      27      1    96%   80
sklearn.metrics.metrics                  248     15    94%   287, 292, 610, 657, 679, 722, 844-845, 847, 917, 952-957, 982-983
sklearn.metrics.pairwise                 156     16    90%   90, 92, 167-169, 239, 244, 385, 489, 513, 552, 631, 640, 650, 655, 658
--------------------------------------------------------------------
TOTAL                                    595     35    94%   
----------------------------------------------------------------------

I just rebuilt and ran the test suite, everything looks fine, I'd say this is good to go (just a featureunion failure that is fixed in master). Sorry but if the history is not pretty enough, I'll leave the squashing to someone else as I will have a busy weekend.

[amueller]

There seems to be a linker problem. I get
ImportError: /home/andy/checkout/scikit-learn/sklearn/metrics/_euclidean_fast.so: undefined symbol: cblas_dsyrk.
Sorry, don't have time to investigate atm.

[amueller]

I couldn't reproduce the error and merged. Thanks for all the work! :)

[amueller]

Ok, so jenkins segfaulted. Thats not good :(

[GaelVaroquaux]

Ok, so jenkins segfaulted. Thats not good :(

Maybe it just needs a cleanup. I think that I got that once on that
branch.

I am trying to see if I can spot a problem on one of my computers.

[GaelVaroquaux]

On Tue, Oct 23, 2012 at 04:03:02PM +0200, Gael Varoquaux wrote:

I am trying to see if I can spot a problem on one of my computers.

Can't spot any problem of course.

How do I trigger a new build (without making a commit) to see if this
goes away :$ ?

[GaelVaroquaux]

On Tue, Oct 23, 2012 at 04:21:17PM +0200, Gael Varoquaux wrote:

On Tue, Oct 23, 2012 at 04:03:02PM +0200, Gael Varoquaux wrote:

I am trying to see if I can spot a problem on one of my computers.

Can't spot any problem of course.

Good news. I can reproduce on one of my computers.

[amueller]

for some definition of "good".
Valgrind?

[GaelVaroquaux]

On Tue, Oct 23, 2012 at 04:49:05PM +0200, Gael Varoquaux wrote:

Good news. I can reproduce on one of my computers.

I think that I found the problem:
#1268

I can't believe that we didn't fall in this trap before...

[GaelVaroquaux]

On Tue, Oct 23, 2012 at 05:19:17PM +0200, Gael Varoquaux wrote:

I think that I found the problem:
#1268

I was wrong.

To reproduce, 'python -c import sklearn.metrics.pairwise' is sufficient.

[GaelVaroquaux]

A backtrace, just in case people have great ideas looking at it:

#0  0x0000000000568af0 in visit_decref ()
#1  0x000000000056953e in list_traverse.16457 ()
#2  0x0000000000543015 in collect.48496 ()
#3  0x00000000005439b9 in _PyObject_GC_Malloc ()
#4  0x00000000004df8cd in list_iter ()
#5  0x00000000004f5602 in PyObject_GetIter ()
#6  0x0000000000440e60 in builtin_zip.32579 ()
#7  0x0000000000497ea4 in PyEval_EvalFrameEx ()
#8  0x000000000049f1c0 in PyEval_EvalCodeEx ()
#9  0x00000000004983b8 in PyEval_EvalFrameEx ()
#10 0x000000000049f1c0 in PyEval_EvalCodeEx ()
#11 0x00000000004983b8 in PyEval_EvalFrameEx ()
#12 0x000000000049f1c0 in PyEval_EvalCodeEx ()
#13 0x00000000004a7f18 in PyImport_ExecCodeModuleEx ()
#14 0x000000000053cde1 in load_source_module.39052 ()
#15 0x000000000056211a in imp_load_module.39102 ()
#16 0x0000000000497ea4 in PyEval_EvalFrameEx ()
#17 0x0000000000498602 in PyEval_EvalFrameEx ()
#18 0x0000000000498602 in PyEval_EvalFrameEx ()
#19 0x000000000049f1c0 in PyEval_EvalCodeEx ()
#20 0x00000000004983b8 in PyEval_EvalFrameEx ()
#21 0x00000000004a7548 in gen_send_ex.isra.0 ()
#22 0x000000000043f6be in listextend.16657 ()
#23 0x0000000000498546 in PyEval_EvalFrameEx ()
#24 0x000000000049f1c0 in PyEval_EvalCodeEx ()

It clearly happens during a collect.