nipy / PySurfer

Cortical neuroimaging visualization in Python
https://pysurfer.github.io/
BSD 3-Clause "New" or "Revised" License
239 stars 98 forks source link

Fast cross product is not that fast #315

Open akapet00 opened 1 year ago

akapet00 commented 1 year ago

Relevant code: https://github.com/nipy/PySurfer/blob/a5a019ec5c6d25cdef92aca84a5e1ce1f9ba6fef/surfer/utils.py#L179-L213

Hi,

I was browsing through the code and I accidentally stumbled upon the function _fast_cross_3d referenced in the URL above. I don't think that this implementation is faster than numpy's implementation of the cross product - for example, please take a look at the following code:

import numpy as np

x = np.random.rand(1_000_000, 3)
y = np.random.rand(1, 3)

def cross_3d_pysurfer(x, y):
    return np.c_[x[:, 1] * y[:, 2] - x[:, 2] * y[:, 1],
                 x[:, 2] * y[:, 0] - x[:, 0] * y[:, 2],
                 x[:, 0] * y[:, 1] - x[:, 1] * y[:, 0]]

def cross_3d_numpy(x, y):
    return np.cross(x, y)

When I time it, I get the following results:

In [1]: %timeit cross_3d_pysurfer(x, y)
Out[1]: 34.7 ms ± 3.9 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [2]: %timeit cross_3d_numpy(x, y)
Out[2]: 29.4 ms ± 4.63 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

A potential, but very small speed up could be achieved by using np.stack instead of np.c_ because it has some cool compiled numpy's stuff which np.c_ lacks. For example:

def cross_3d_pysurfer_modified(x, y):
    return np.stack((x[:, 1] * y[:, 2] - x[:, 2] * y[:, 1],
                     x[:, 2] * y[:, 0] - x[:, 0] * y[:, 2],
                     x[:, 0] * y[:, 1] - x[:, 1] * y[:, 0])).T
In [3]: %timeit cross_3d_pysurfer_modified(x, y)
Out[3]: 28.4 ms ± 3.95 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

I didn't want to open a PR because I really don't know if this is even relevant anymore or if the difference in a few ms is that important, but still wanted to let you know :)

larsoner commented 1 year ago

FYI I think the motivation was that numpy's cross did not originally broadcast. Now that it does we could replace it. But we'd want to make sure that it doesn't require a super new NumPy (e.g., 1.23) and works with some older ones (back to 1.19 or 1.20). But from looking at the NumPy doc it seems like it was supported all the way back in 1.13, so yes we could change this!