Fast scalar triple product in numpy

Question

I have a large number of vector triples, and I would like to compute the scalar triple product for them. I can do

import numpy

n = 871
a = numpy.random.rand(n, 3)
b = numpy.random.rand(n, 3)
c = numpy.random.rand(n, 3)

# <a, b x c>
omega = numpy.einsum('ij, ij->i', a, numpy.cross(b, c))

but numpy.cross is fairly slow. The symmetry of the problem (its Levi-Civita expression is eps_{ijk} a_i b_j c_k) suggests that there might be a better (faster) way to compute it, but I can't seem to figure it out.

Any hints?

Answer 1

It's just the determinant.

omega=det(dstack([a,b,c]))

But it is slower....

An other equivalent solution is omega=dot(a,cross(b,c)).sum(1) .

But I think you have to compute about 9 (for cross) + 3 (for dot) + 2 (for sum) = 14 operations for each det, so it seems to be near optimal. At best you will win a two factor in numpy.

EDIT:

If speed is critical, you must go at low level. numba is a easy way to do that for a 15X factor here :

from numba import njit

@njit
def multidet(a,b,c):
    n=a.shape[0]
    d=np.empty(n)
    for i in range(n):
        u,v,w=a[i],b[i],c[i]
        d[i]=\
        u[0]*(v[1]*w[2]-v[2]*w[1])+\
        u[1]*(v[2]*w[0]-v[0]*w[2])+\
        u[2]*(v[0]*w[1]-v[1]*w[0])  # 14 operations / det
    return d

some tests:

In [155]: %timeit multidet(a,b,c)
100000 loops, best of 3: 7.79 µs per loop

In [156]: %timeit numpy.einsum('ij, ij->i', a, numpy.cross(b, c))
10000 loops, best of 3: 114 µs per loop


In [159]: allclose(multidet(a,b,c),omega)
Out[159]: True