numpy中的快速標量三元產品

我有大量的向量三元組，我想爲它們計算scalar triple product。我可以做numpy中的快速標量三元產品

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

但numpy.cross是相當緩慢的。問題的對稱性（Levi-Civita表達式爲eps_{ijk} a_i b_j c_k）表明可能有更好（更快）的方法來計算它，但我似乎無法弄清楚。

任何提示？

來源

2017-02-10 Nico Schlömer

過去的問題已經表明，它做什麼，目前的'cross'是有效的，http://stackoverflow.com/q/ 39662540/901925 – hpaulj

單獨的交叉產品很好，這是一個點產品，使其具有特定的對稱性。這讓我覺得奇怪，應該沒有進一步的優化。 –

如果A，B和C是特定的3向量，那麼當epsilon是Levi-Civita符號時，您正在尋找的愛因斯坦符號是'\ epsilon_ {ijk} A^i B^j C^k'。它不認爲它已經在'numpy.einsum'中實現了。作爲@B的。 M.說，這與'np.dstack（A，B，C）'的設計完全相同。 – tsvikas

我已經做了在答案中提到的方法進行了比較。結果：

@ Divakar的節拍由位einsum交叉之一。

爲了完整起見，讓我注意到還有另一種方法完全依賴點積和sqrt，請參閱here。該方法比雙重交叉和切片總和稍慢。

情節與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, 
     )

來源

2017-02-22 08:32:35

這只是決定因素。

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

但它是較慢....

的其它等效的解決方案是omega=dot(a,cross(b,c)).sum(1)。

但我認爲你必須爲每個det計算約9（對於十字）+3（對於點）+2（對於總和）= 14操作，所以它似乎接近最優。充其量你會贏得雙倍的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

來源

2017-02-10 11:38:31

從'numpy.linalg.det'：'LinAlgError：數組的最後2個維度必須是正方形。請注意，我需要一次計算許多標量三重產品。 –

修復堆棧符號，這是一個非常好的答案 –

我糾正了det版本，但沒有優化。我添加其他解決方案，希望它可以幫助。 –

這裏有一個方法利用的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方法，我們直接進入那個(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

運行測試 -

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

來源

2017-02-10 11:42:05 Divakar

numpy中的快速標量三元產品

回答

相關問題