如何加速向量叉积的计算?

11

嗨,我是相对新的,尝试使用numpy进行一些计算。我在一个特定的计算中遇到了很长的经过时间,无法找到更快的方法来完成同样的事情。

基本上,这是光线三角形相交算法的一部分,我需要从两个不同大小的矩阵中计算所有向量叉积。

我使用的代码是:

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

其中dirvectors是一个包含n*向量(xyz)的数组,trivectors2是一个包含m*向量(xyz)的数组。allhvals1是一个大小为n*M*向量(xyz)的叉积数组。

这个方法能够工作,但速度很慢。本质上,它是每个数组中的每个向量形成的n*m矩阵。希望你明白。每个数组的大小因参数不同而变化,大约在1-4000之间(我基本上根据大小分块dirvectors)。

欢迎任何建议。不幸的是,我的矩阵数学有点薄弱。


1
虽然不想唱反调,但这并不是一个论坛 :) 我之所以提到这个问题是因为有太多人把这个网站当作一个论坛。不过,你的问题没有任何问题。 - keyser
有可能大部分时间都花在提取向量上,而不是叉积运算上。我建议在进行叉积运算之前,先将它们提取到变量中。然后可以使用此技术来获得更好的洞察力。 - Mike Dunlavey
2个回答

20
如果您查看 源代码 中的 np.cross,它基本上将所有数组的 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

我们可以尝试使用Levi-Civita符号和np.einsum来计算它,具体方法如下:

使用Levi-Civita符号

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

太棒了!接近一个数量级的提升...

以下是NumPy 1.11.0的一些性能数据,使用a=numpy.random.rand(n,3)b=numpy.random.rand(n,3)

enter image description here

对于测试过的最大的n值,嵌套的einsumcross快约两倍。


顺便提一下,感谢@hpaulj昨天的评论,提醒我Levi-Civita符号的相关知识。 - Jaime
直到现在我才知道np.einsum,非常酷。现在我感到失望的是,我的大学课程中唯一提到Levi-Civita符号的只是在讨论表示叉积的各种方法时顺便提到过(或者可能是关于矩阵行列的排列组合之类的东西;我们之后从未回到它们,所以我记得不太清楚)。这就是我成为工程师而不是数学家的代价,我想。 - JAB
我会尝试一下并看看它的表现。我曾经看过einsum,但是对它不太了解。 - user1942439
只是确认这个很好地工作了。在这个例程的其余部分仍然有一些加速工作要做,但这确实帮了很多忙。 - user1942439
如果您发现这个回答解决了您的问题,您可以考虑点击其旁边的复选框接受它... - Jaime
显示剩余3条评论

0

1
最好将Matlab代码翻译成Python代码。 - zangw
虽然提供链接很好,但请为您的链接提供上下文,因为一个答案实际上应该包含一个答案。请注意,仅仅是一个指向外部网站的链接可能是为什么和如何删除某些答案?的原因之一。 - Jan

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