Skip to content

Fast cross product is not that fast #315

Open

Description

@akapet00

Relevant code:

PySurfer/surfer/utils.py

Lines 179 to 213 in a5a019e

def _fast_cross_3d(x, y):
"""Compute cross product between list of 3D vectors
Much faster than np.cross() when the number of cross products
becomes large (>500). This is because np.cross() methods become
less memory efficient at this stage.
Parameters
----------
x : array
Input array 1.
y : array
Input array 2.
Returns
-------
z : array
Cross product of x and y.
Notes
-----
x and y must both be 2D row vectors. One must have length 1, or both
lengths must match.
"""
assert x.ndim == 2
assert y.ndim == 2
assert x.shape[1] == 3
assert y.shape[1] == 3
assert (x.shape[0] == 1 or y.shape[0] == 1) or x.shape[0] == y.shape[0]
if max([x.shape[0], y.shape[0]]) >= 500:
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]]
else:
return np.cross(x, y)

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 :)

Activity

  1. larsoner commented on Sep 2, 2022

    @larsoner
    Contributor

    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!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions