如何加快矢量叉积计算

Question

嗨,我在这里比较新,并尝试用numpy做一些计算.我从一个特定的计算中经历了很长的时间,并且无法以任何更快的方式实现同​​样的事情.

基本上它是光线三角形交叉算法的一部分,我需要从两个不同大小的矩阵计算所有矢量交叉产品.

我使用的代码是:

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

其中dirvectors是的阵列n* vectors (xyz)和trivectors2是的阵列m*vectors(xyz).allhvals1是一个大小的交叉产品的数组n*M*vector (xyz).这有效,但速度很慢.它本质上是来自每个阵列的每个向量的n*m矩阵.希望你明白.根据参数的不同,每个的大小从大约1-4000变化(我基本上根据大小对dirvectors进行分块).

任何建议表示赞赏 不幸的是,我的矩阵数学有点不稳定.

Answer 1

如果你看的源代码的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

我们可以尝试计算它使用列维-奇维塔符号和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

答对了!接近一个数量级的改进......

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

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