看起来显式计算向量数组的叉积比使用 np.cross
快得多。我已经尝试过 vector-first 和 vector-last ,它似乎没有什么区别,尽管这是在对类似 question 的回答中提出的。是我用错了,还是速度变慢了?
在笔记本电脑上,每个交叉产品的显式计算似乎需要大约 60ns。这是〜大约〜一样快吗?在这种情况下,似乎还没有任何理由去 Cython 或 PyPy 或编写特殊的 ufunc
。
我也看到了使用 einsum 的引用资料,但我真的不明白如何使用它,并且怀疑它不是更快。
a = np.random.random(size=300000).reshape(100000,3) # vector last
b = np.random.random(size=300000).reshape(100000,3)
c, d = a.swapaxes(0, 1), b.swapaxes(0, 1) # vector first
def npcross_vlast(): return np.cross(a, b)
def npcross_vfirst(): return np.cross(c, d, axisa=0, axisb=0)
def npcross_vfirst_axisc(): return np.cross(c, d, axisa=0, axisb=0, axisc=0)
def explicitcross_vlast():
e = np.zeros_like(a)
e[:,0] = a[:,1]*b[:,2] - a[:,2]*b[:,1]
e[:,1] = a[:,2]*b[:,0] - a[:,0]*b[:,2]
e[:,2] = a[:,0]*b[:,1] - a[:,1]*b[:,0]
return e
def explicitcross_vfirst():
e = np.zeros_like(c)
e[0,:] = c[1,:]*d[2,:] - c[2,:]*d[1,:]
e[1,:] = c[2,:]*d[0,:] - c[0,:]*d[2,:]
e[2,:] = c[0,:]*d[1,:] - c[1,:]*d[0,:]
return e
print "explicit"
print timeit.timeit(explicitcross_vlast, number=10)
print timeit.timeit(explicitcross_vfirst, number=10)
print "np.cross"
print timeit.timeit(npcross_vlast, number=10)
print timeit.timeit(npcross_vfirst, number=10)
print timeit.timeit(npcross_vfirst_axisc, number=10)
print all([npcross_vlast()[7,i] == npcross_vfirst()[7,i] ==
npcross_vfirst_axisc()[i,7] == explicitcross_vlast()[7,i] ==
explicitcross_vfirst()[i,7] for i in range(3)]) # check one
explicit
0.0582590103149
0.0560920238495
np.cross
0.399816989899
0.412983894348
0.411231040955
True
最佳答案
首先,如果您希望加快代码速度,您可能应该尝试完全摆脱交叉产品。在许多情况下这是可能的,例如,当与点积 <a x b, c x d> = <a, c><b, d> - <a, d><b, c>
结合使用时。
无论如何,如果您真的需要明确的交叉产品,请查看
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
np.einsum('ijk,aj,ak->ai', eijk, a, b)
np.einsum('iak,ak->ai', np.einsum('ijk,aj->iak', eijk, a), b)
这两个等效于
np.cross
,其中第二个使用两个 einsum,每个 einsum 带有两个参数,这是 a similar question 中建议的一种技术。然而,结果令人失望:这两种变体都比
np.cross
慢(除了微小的 n
):该情节是用
import numpy as np
import perfplot
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
perfplot.show(
setup=lambda n: np.random.rand(2, n, 3),
n_range=[2**k for k in range(14)],
kernels=[
lambda X: np.cross(X[0], X[1]),
lambda X: np.einsum('ijk,aj,ak->ai', eijk, X[0], X[1]),
lambda X: np.einsum('iak,ak->ai', np.einsum('ijk,aj->iak', eijk, X[0]), X[1]),
],
labels=['np.cross', 'einsum', 'double einsum'],
xlabel='len(a)',
logx=True,
logy=True,
)
关于python-2.7 - 获得叉积的最快方法,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/29713321/