如何加快向量叉积计算

Question

嗨，我在这里比较新，正在尝试用 numpy 进行一些计算。我从一个特定的计算中经历了很长的时间，无法找到更快的方法来实现同样的目标。

基本上它是射线三角形相交算法的一部分，我需要从两个不同大小的矩阵中计算所有矢量 cros 产品。

我使用的代码是：

allhvals1 = numpy.cross( dirvectors[:,None,:], trivectors2[None,:,:] )

其中dirvectors 是n* vectors (xyz) 的数组，trivectors2 是m*vectors(xyz) 的数组。 allhvals1 是大小为 n*M*vector (xyz) 的叉积数组。 这有效，但速度很慢。它本质上是每个数组中每个向量的 n*m 矩阵。希望你能理解。每个的大小从大约 1 到 4000 不等，具体取决于参数（我基本上根据大小对 dirvector 进行分块）。

任何建议表示赞赏。不幸的是，我的矩阵数学有点古怪。

Answer 1

如果您查看np.cross 的the source code，它基本上将xyz 维度移动到所有数组的形状元组的前面，然后将每个组件的计算拼写如下：

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]

在您的情况下，这些产品中的每一个都需要分配巨大的数组，因此整体行为效率不高。

让我们设置一些测试数据：

u = np.random.rand(1000, 3)
v = np.random.rand(2000, 3)

In [13]: %timeit s1 = np.cross(u[:, None, :], v[None, :, :])
1 loops, best of 3: 591 ms per loop

我们可以尝试计算using Levi-Civita symbols和np.einsum如下：

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

In [14]: %timeit s2 = np.einsum('ijk,uj,vk->uvi', eijk, u, v)
1 loops, best of 3: 706 ms per loop

In [15]: np.allclose(s1, s2)
Out[15]: True

因此，虽然它有效，但性能更差。问题是np.einsum 在有两个以上操作数时会遇到麻烦，但优化了两个或更少的路径。所以我们可以尝试分两步重写，看看是否有帮助：

In [16]: %timeit s3 = np.einsum('iuk,vk->uvi', np.einsum('ijk,uj->iuk', eijk, u), v)
10 loops, best of 3: 63.4 ms per loop

In [17]: np.allclose(s1, s3)
Out[17]: True

宾果！接近一个数量级的改进......

带有a=numpy.random.rand(n,3)、b=numpy.random.rand(n,3) 的 NumPy 1.11.0 的一些性能数据：

对于测试过的最大n，嵌套的einsum 的速度大约是cross 的两倍。

Answer 2

在为水下航行器编写动态模拟时，我发现了这种快速交叉乘积的方法：

https://github.com/simena86/Simulink-Underwater-Robotics-Simulator/blob/master/3rdparty/gnc_mfiles/Smtrx.m

效果很好，它是用 Matlab 编写的，但代码非常简单。只需阅读顶部的 cmets。