Fast scalar 3D product in numpy

Question

I have a large number of vector triples, and I would like to calculate a scalar three-local 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.crossrather slow. The symmetry of the problem (its expression is Levi-Civita eps_{ijk} a_i b_j c_k) suggests that there may be a better (faster) way to calculate it, but I cannot figure out how to understand it.

Any clues?

Answer 1

I made a comparison of the methods mentioned in the answers. Results:

@Divakar bits einsum-cross one per bit.

, , ​​sqrt, . . , einsum-cross, slice-sum.

---

perfplot,

import numpy
import perfplot


def einsum_cross(data):
    a, b, c = data
    return numpy.einsum('ij, ij->i', a, numpy.cross(b, c))


def det(data):
    a, b, c = data
    return numpy.linalg.det(numpy.dstack([a, b, c]))


def slice_sum(data):
    a, b, c = data
    c0 = b[:, 1]*c[:, 2] - b[:, 2]*c[:, 1]
    c1 = b[:, 2]*c[:, 0] - b[:, 0]*c[:, 2]
    c2 = b[:, 0]*c[:, 1] - b[:, 1]*c[:, 0]
    return a[:, 0]*c0 + a[:, 1]*c1 + a[:, 2]*c2


perfplot.show(
        setup=lambda n: (
            numpy.random.rand(n, 3),
            numpy.random.rand(n, 3),
            numpy.random.rand(n, 3)
            ),
        n_range=[2**k for k in range(1, 20)],
        kernels=[einsum_cross, det, slice_sum],
        logx=True,
        logy=True,
        )

Answer 2

.

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

....

omega=dot(a,cross(b,c)).sum(1).

, 9 ( ) + 3 ( ) + 2 ( ) = 14 det, . numpy.

EDIT:

, . numba - 15X:

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

:

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

Answer 3

, slicing -

def slicing_summing(a,b,c):
    c0 = b[:,1]*c[:,2] - b[:,2]*c[:,1]
    c1 = b[:,2]*c[:,0] - b[:,0]*c[:,2]
    c2 = b[:,0]*c[:,1] - b[:,1]*c[:,0]
    return a[:,0]*c0 + a[:,1]*c1 + a[:,2]*c2

We can replace the first three steps that calculate c0, c1, c2its complex version with a single-line layer, for example:

b[:,[1,2,0]]*c[:,[2,0,1]] - b[:,[2,0,1]]*c[:,[1,2,0]]

This will create another array (n,3)that should be used with afor summing, resulting in an array (n,). Using the proposed method, slicing_summingwe directly fall into this array (n,)with the summation of these three slices and thereby avoid this intermediate array (n,3).

Run Example -

In [86]: # Setup inputs    
    ...: n = 871
    ...: a = np.random.rand(n, 3)
    ...: b = np.random.rand(n, 3)
    ...: c = np.random.rand(n, 3)
    ...: 

In [87]: # Original approach
    ...: omega = np.einsum('ij, ij->i', a, np.cross(b, c))

In [88]: np.allclose(omega, slicing_summing(a,b,c))
Out[88]: True

Runtime Test -

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

In [91]: %timeit slicing_summing(a,b,c)
1000 loops, best of 3: 63 µs per loop