Fast cross product is not that fast

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

[larsoner]

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!