Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

how to speed up a vector cross product calculation

Hi I'm relatively new here and trying to do some calculations with numpy. I'm experiencing a long elapse time from one particular calculation and can't work out any faster way to achieve the same thing.

Basically its part of a ray triangle intersection algorithm and I need to calculate all the vector cros products from two matrices of different sizes.

The code I was using was :

allhvals1 = numpy.cross( dirvectors[:,None,:], trivectors2[None,:,:] )

where dirvectors is an array of n* vectors (xyz) and trivectors2 is an array of m*vectors(xyz). allhvals1 is an array of the cross products of size n*M*vector (xyz). This works but is very slow. It's essentially the n*m matrix of each vector from each array. Hope that you understand. The sizes of each varies from approx 1-4000 depending on parameters (I basically chunk the dirvectors dependent on size).

Any advice appreciated. Unfortunately my matrix math is somewhat flakey.

like image 919
user1942439 Avatar asked Jan 03 '14 16:01

user1942439


People also ask

How is AxB calculated?

(See section 3-7 in the text for more review.) Magnitude: |AxB| = A B sinθ. Just like the dot product, θ is the angle between the vectors A and B when they are drawn tail-to-tail. Direction: The vector AxB is perpendicular to the plane formed by A and B.

Does order matter for cross product of vectors?

1. a × b = −(b × a). In words, the order does matter! Calculating b × a gives us a vector in the opposite direction as a × b.


1 Answers

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.

like image 106
Jaime Avatar answered Oct 24 '22 18:10

Jaime