获得叉积的最快方法

Question

看起来显式计算向量数组的叉积比使用 快得多np.cross。我尝试过向量优先和向量最后，它似乎没有什么区别，尽管这是在类似问题的答案中提出的。是我用错了，还是速度慢了？

在笔记本电脑上，显式计算每个叉积大约需要 60 纳秒。这大概是最快的速度吗？在这种情况下，似乎还没有任何理由去使用 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

Answer 1

首先，如果您希望加快代码速度，您可能应该尝试完全摆脱交叉产品。在许多情况下这是可能的，例如，当与点积结合使用时<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)

这两个相当于，其中第二个使用两个 einsum，每个 einsum 有两个参数，这是类似问题np.cross中建议的技术。

但结果令人失望：这两个变体都慢于np.cross（tiny 除外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


b = perfplot.bench(
    setup=lambda n: np.random.rand(2, n, 3),
    n_range=[2 ** k for k in range(23)],
    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)",
)

b.save("out.png")