获得叉积的最快方法

5
看起来,对于一个向量数组的叉积计算而言,显式计算要比使用np.cross快得多。我已经尝试了矢量-先和矢量-后,但似乎没有区别,尽管在类似问题的答案中提出了这种方法。我是在使用错误的方法吗,还是它只是更慢?
在笔记本电脑上,显式计算每个叉积需要大约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

看一下 np.cross 的代码。它正在做你正在做的事情,还加入了一些处理大小为2的情况的覆盖,并进行了一些轴交换,以便可以使用像 a[1]*b[2] - a[2]*b[1] 这样的表达式。只要大维度是向量化的,对小维度(大小为3)进行一些明确的步骤不会影响速度。 - hpaulj
我的一个问题是:为什么np.cross的速度几乎慢了10倍,与大小或顺序无关? - uhoh
1
正如@Jaime所暗示的那样,更新numpy可能会解决这个问题。我在1.9.2上看到非常相似的时间。 - cel
1
swapaxes 对速度没有任何影响,因为内存布局仍然相同。如果数组从一开始就是这样生成的话,vfirst 会稍微快一些。 - hpaulj
如下答案所述,我只有1.7.1版本。@jpaulj感谢您指出swapaxes仅返回视图。 - uhoh
显示剩余2条评论
3个回答

4
np.cross函数在numpy的1.9.x版本中性能有了显著提升。
%timeit explicitcross_vlast()
%timeit explicitcross_vfirst()
%timeit npcross_vlast()
%timeit npcross_vfirst()
%timeit npcross_vfirst_axisc() 

这是我在1.8.0版本中得到的时间记录。
100 loops, best of 3: 4.47 ms per loop
100 loops, best of 3: 4.41 ms per loop
10 loops, best of 3: 29.1 ms per loop
10 loops, best of 3: 29.3 ms per loop
10 loops, best of 3: 30.6 ms per loop

以下是版本号为1.9.0的时间表:
100 loops, best of 3: 4.62 ms per loop
100 loops, best of 3: 4.19 ms per loop
100 loops, best of 3: 4.05 ms per loop
100 loops, best of 3: 4.09 ms per loop
100 loops, best of 3: 4.24 ms per loop

我怀疑这个加速是由合并请求#4338引入的。


谢谢@cel。时间飞逝 - 我回到了1.7.1,吸取了教训! - uhoh
2
这曾经是被接受的答案 - 实际上这就是我需要的答案(在np.cross()被改进之前,我使用的是旧版本的numpy)。但是我已经转向@NicoSchlömer的答案 - 我认为这对于那些正在寻找实现叉积的最快方法的人来说是最有用的信息。再次感谢您的帮助! - uhoh

2
首先,如果您想加快代码速度,您应该尽可能消除叉积。在许多情况下,这是可能的,例如,在与点积一起使用时<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技巧,这是在类似问题中建议使用的
然而,结果令人失望:这两个变体都比np.cross慢(只有当n很小的时候除外): enter image description here 该图由以下代码生成:
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")

非常有用的绘图!当我使用许多个三维向量的叉积进行蒙特卡罗模拟时,np.einsum() 至少似乎具有一定的优势。由于我最初提出的问题的缓慢部分是因为旧版本的 numpy 在 np.cross() 速度上得到了改进,您能否记录一下您测试过的版本(供记录)? - uhoh

1

将您的vlast简单更改为

def stacked_vlast(a,b):
        x = a[:,1]*b[:,2] - a[:,2]*b[:,1]
        y = a[:,2]*b[:,0] - a[:,0]*b[:,2]
        z = a[:,0]*b[:,1] - a[:,1]*b[:,0]
        return np.array([x,y,z]).T

即用堆叠替换列分配,就像(旧的)cross函数所做的那样,会使速度减慢5倍。

当我使用开发中的cross函数的本地副本时,与您的explicit_vlast相比,我获得了轻微的速度提升。该cross使用out参数试图减少临时数组,但我的简单测试表明它在速度上并没有太大的差异。

https://github.com/numpy/numpy/blob/master/numpy/core/numeric.py

如果你的显式版本可以工作,我不会升级numpy来获取这个新的cross

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接