如何加快向量叉积计算 [英] how to speed up a vector cross product calculation

问题描述

我在这里相对较新，正在尝试使用numpy进行一些计算.我从一种特定的计算中经历了很长的时间，无法找到更快的方法来实现相同的目的.

Hi I'm relatively new here and trying to do some calculations with numpy. I'm experiencing a long elapse time from one particular calculation and can't work out any faster way to achieve the same thing.

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

Basically its part of a ray triangle intersection algorithm and I need to calculate all the vector cros products from two matrices of different sizes.

我使用的代码是:

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

其中，dirvectors是n* vectors (xyz)的数组，而trivectors2是m*vectors(xyz)的数组. allhvals1是大小为n*M*vector (xyz)的叉积的数组. 这行得通，但是非常慢.本质上，它是每个数组中每个向量的n * m矩阵.希望你能理解.每个参数的大小在大约1-4000之间变化(根据参数，我基本上将目录向量分块，具体取决于大小).

where dirvectors is an array of n* vectors (xyz) and trivectors2 is an array of m*vectors(xyz). allhvals1 is an array of the cross products of size n*M*vector (xyz). This works but is very slow. It's essentially the n*m matrix of each vector from each array. Hope that you understand. The sizes of each varies from approx 1-4000 depending on parameters (I basically chunk the dirvectors dependent on size).

任何建议表示赞赏.不幸的是，我的矩阵数学有点夸张.

Any advice appreciated. Unfortunately my matrix math is somewhat flakey.

推荐答案

如果您查看

If you look at the source code of np.cross, it basically moves the xyz dimension to the front of the shape tuple for all arrays, and then has the calculation of each of the components spelled out like this:

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]

在您的情况下，这些产品中的每一个都需要分配巨大的数组，因此总体行为不是很有效.

In your case, each of those products requires allocating huge arrays, so the overall behavior is not very efficient.

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

Lets set up some test data:

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如下:

We can try to compute it using Levi-Civita symbols and np.einsum as follows:

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在有两个以上操作数时会遇到麻烦，但对两个或更少的操作数具有优化的路径.因此，我们可以尝试分两步重写它，以查看是否有帮助:

So while it works, it has worse performance. The thing is that np.einsum has trouble when there are more than two operands, but has optimized pathways for two or less. So we can try to rewrite it in two steps, to see if it helps:

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

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

Bingo! Close to an order of magnitude improvement...

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

Some performance figures for NumPy 1.11.0 with a=numpy.random.rand(n,3), b=numpy.random.rand(n,3):

嵌套的einsum大约是测试的最大n的cross的两倍.

The nested einsum is about twice as fast as cross for the largest n tested.

这篇关于如何加快向量叉积计算的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！