how to speed up a vector cross product calculation

by

If you look at the source code of np.cross, it basically moves the xyz dimension to the front of the shape tuple for all arrays, and then has the calculation of each of the components spelled out like this:

x = a[1]*b[2] - a[2]*b[1]
y = a[2]*b[0] - a[0]*b[2]
z = a[0]*b[1] - a[1]*b[0]

In your case, each of those products requires allocating huge arrays, so the overall behavior is not very efficient.

Lets set up some test data:

u = np.random.rand(1000, 3)
v = np.random.rand(2000, 3)

In [13]: %timeit s1 = np.cross(u[:, None, :], v[None, :, :])
1 loops, best of 3: 591 ms per loop

We can try to compute it using Levi-Civita symbols and np.einsum as follows:

eijk = np.zeros((3, 3, 3))
eijk[0, 1, 2] = eijk[1, 2, 0] = eijk[2, 0, 1] = 1
eijk[0, 2, 1] = eijk[2, 1, 0] = eijk[1, 0, 2] = -1

In [14]: %timeit s2 = np.einsum('ijk,uj,vk->uvi', eijk, u, v)
1 loops, best of 3: 706 ms per loop

In [15]: np.allclose(s1, s2)
Out[15]: True

So while it works, it has worse performance. The thing is that np.einsum has trouble when there are more than two operands, but has optimized pathways for two or less. So we can try to rewrite it in two steps, to see if it helps:

In [16]: %timeit s3 = np.einsum('iuk,vk->uvi', np.einsum('ijk,uj->iuk', eijk, u), v)
10 loops, best of 3: 63.4 ms per loop

In [17]: np.allclose(s1, s3)
Out[17]: True

Bingo! Close to an order of magnitude improvement…

Some performance figures for NumPy 1.11.0 with a=numpy.random.rand(n,3), b=numpy.random.rand(n,3):

enter image description here

The nested einsum is about twice as fast as cross for the largest n tested.

More Related Contents:

  1. Most efficient way to map function over numpy array
  2. Performance of Pandas apply vs np.vectorize to create new column from existing columns
  3. numpy: most efficient frequency counts for unique values in an array
  4. Why is numpy’s einsum faster than numpy’s built in functions?
  5. Roll rows of a matrix independently
  6. Numpy sum elements in array based on its value
  7. Efficiently return the index of the first value satisfying condition in array
  8. Bin elements per row – Vectorized 2D Bincount for NumPy
  9. Fastest way to grow a numpy numeric array
  10. Vectorized searchsorted numpy
  11. Frequency counts for unique values in a NumPy array
  12. Most efficient way to forward-fill NaN values in numpy array
  13. Numpy: Fix array with rows of different lengths by filling the empty elements with zeros
  14. Pandas pd.Series.isin performance with set versus array
  15. Why is numpy.array so slow?
  16. Efficient evaluation of a function at every cell of a NumPy array
  17. Efficiently detect sign-changes in python
  18. Efficiently Calculating a Euclidean Distance Matrix Using Numpy
  19. Efficient way to take the minimum/maximum n values and indices from a matrix using NumPy
  20. why isn’t numpy.mean multithreaded?
  21. Numpy vs Cython speed
  22. numpy float: 10x slower than builtin in arithmetic operations?
  23. Fastest save and load options for a numpy array
  24. No speedup when summing uint16 vs uint64 arrays with NumPy?
  25. Efficiently replace elements in array based on dictionary – NumPy / Python
  26. Python Numpy Data Types Performance
  27. Getting days since last occurence in Pandas DataFrame?
  28. What is the most efficient way to check if a value exists in a NumPy array?
  29. Efficient dot products of large memory-mapped arrays
  30. What is the most efficient way of doing square root of sum of square of two numbers?

Leave a Comment