Быстрый скалярный трехмерный продукт в numpy

У меня большое количество векторных троек, и я хотел бы вычислить для них скалярное трехмерное произведение . я могу сделать

import numpy n = 871 a = numpy.random.rand(n, 3) b = numpy.random.rand(n, 3) c = numpy.random.rand(n, 3) # <a, bx c> omega = numpy.einsum('ij, ij->i', a, numpy.cross(b, c)) 

но numpy.cross довольно медленный. Симметрия проблемы (ее выражение Леви-Чивита – eps_{ijk} a_i b_j c_k ) предполагает, что может быть лучший (более быстрый) способ вычислить ее, но я не могу понять ее.

Любые намеки?

Это просто детерминант.

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

Но это медленнее ….

Другим эквивалентным решением является omega=dot(a,cross(b,c)).sum(1) .

Но я думаю, что вам нужно вычислить около 9 (для креста) + 3 (для точки) + 2 (для суммы) = 14 операций для каждого det, поэтому он кажется почти оптимальным. В лучшем случае вы выиграете два фактора в numpy.

РЕДАКТИРОВАТЬ:

Если скорость критическая, вы должны идти на низком уровне. 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 

Вот один из способов использования 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 

Мы можем заменить первые три шага, которые вычисляют c0, c1, c2 и его сложную версию с одним слоем, например,

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

Это создало бы другой (n,3) массив, который должен использоваться с a для суммирования, приводящий к массиву (n,) . С предлагаемым методом slicing_summing мы непосредственно slicing_summing до этой (n,) формы массива с суммированием этих трех срезов и, таким образом, избегая этого промежуточного (n,3) массива.

Пример прогона –

 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 

Я провел сравнение методов, упомянутых в ответах. Результаты:

введите описание изображения здесь

@ Дивакар бьет бит-эйнсум-крест понемногу.

Для полноты позвольте заметить, что существует другой метод, основанный исключительно на точечных продуктах и ​​sqrt, см. Здесь . Этот метод немного медленнее, чем как einsum-cross, так и slice-sum.


Сюжет был создан с помощью перфолта ,

 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, )