获得叉积的最快方法答案

【问题标题】：fastest way to obtain cross product获得叉积的最快方法
【发布时间】：2015-06-25 03:19:05
【问题描述】：

看起来显式计算向量数组的叉积比使用np.cross 快得多。我已经尝试过矢量优先和矢量最后，它似乎没有什么区别，尽管这是在类似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

【问题讨论】：

查看np.cross的代码。它正在做你正在做的事情，加上一些处理大小为 2 的情况的掩护，以及一些轴交换，因此它可以使用像 a[1]*b[2] - a[2]*b[1] 这样的表达式。只要对大维度进行矢量化处理，在小维度（大小 3）上执行一些显式步骤就不会影响您的速度。
（其中一个）我的问题是：为什么 np.cross 几乎慢了 10 倍，与大小或顺序无关？
正如@Jaime 暗示的那样，更新numpy 可能会解决这个问题。我在1.9.2 看到非常相似的时间。
之前的跨积速度问题：stackoverflow.com/questions/20908754/…; stackoverflow.com/questions/28132147/…;
swapaxes 对速度没有任何帮助，因为内存布局仍然相同。如果从一开始就以这种方式生成数组，vfirst 会稍微快一些。

标签： python-2.7 numpy numpy-ufunc

【解决方案1】：

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，吸取了教训！
这曾经是公认的答案 - 实际上是我需要的答案（在改进 np.cross() 之前，我使用的是旧版本的 numpy）。但是我已经切换到@NicoSchlömer 的答案 - 我认为对于正在寻找实现交叉产品的最快方法的人来说，这是最有用的信息找到这个问题。再次感谢您的帮助！

【解决方案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，每个带有两个参数，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


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

【讨论】：

非常丰富的情节！当我用大量单个 3 向量的叉积进行蒙特卡罗时，np.einsum() 似乎至少有一个有用的优势。由于我最初的问题中提到的缓慢部分是因为在 np.cross() 速度提高之前的旧版本的 numpy，你能注意（记录在案）你测试了哪个版本吗？

【解决方案3】：

只需将您的 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

如果您的显式版本有效，我不会为了获得这个新的cross 而升级numpy。

【讨论】：