user1942439 - 2 months ago 13

Python Question

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`

`n* vectors (xyz)`

`trivectors2`

`m*vectors(xyz)`

`allhvals1`

`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.

Answer

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 with `a=numpy.random.rand(n,3)`

, `b=numpy.random.rand(n,3)`

:

The nested `einsum`

is about twice as fast as `cross`

for the largest `n`

tested.