2向量矩阵点积的性能答案

【问题标题】：Performance of 2 vector-matrix dot product2向量矩阵点积的性能
【发布时间】：2015-09-07 19:52:10
【问题描述】：

这个问题更侧重于计算的性能。

我有 2 个向量矩阵。这意味着它们对于 X、Y、Z 具有 3 个深度维度。矩阵的每个元素都必须与另一个矩阵的相同位置的元素进行点积。

一个简单而低效的代码将是这个：

import numpy as np
a = np.random.uniform(low=-1.0, high=1.0, size=(1000,1000,3))
b = np.random.uniform(low=-1.0, high=1.0, size=(1000,1000,3))
c = np.zeros((1000,1000))
numRow,numCol,numDepth = np.shape(a)

for idRow in range(numRow):
    for idCol in range(numCol):
        # Angle in radians
        c[idRow,idCol] = math.acos(a[idRow,idCol,0]*b[idRow,idCol,0] + a[idRow,idCol,1]*b[idRow,idCol,1] + a[idRow,idCol,2]*b[idRow,idCol,2])

但是，numpy 函数可以加快计算速度，使代码更快：

# Angle in radians
d = np.arccos(np.multiply(a[:,:,0],b[:,:,0]) + np.multiply(a[:,:,1],b[:,:,1]) + np.multiply(a[:,:,2],b[:,:,2]))

但是，我想知道是否还有其他 sintaxis 可以通过其他功能、索引等来改进上面的这一点...

第一个代码需要 4.658s 而第二个代码需要 0.354s

【问题讨论】：

标签： python performance numpy dot-product

【解决方案1】：

您可以使用 np.einsum 执行此操作，它会在任何轴上相乘然后求和：

np.arccos(np.einsum('ijk,ijk->ij', a, b))

执行您在问题中发布的更直接的方法是使用np.sum，沿最后一个轴求和 (-1)：

np.arccos(np.sum(a*b, -1))

他们都给出了相同的答案，但einsum 最快，sum 次之：

In [36]: timeit np.arccos(np.einsum('ijk,ijk->ij', a, b))
10000 loops, best of 3: 20.4 µs per loop

In [37]: timeit e = np.arccos(np.sum(a*b, -1))
10000 loops, best of 3: 29.8 µs per loop

In [38]: %%timeit
   ....: d = np.arccos(np.multiply(a[:,:,0],b[:,:,0]) +
   ....:               np.multiply(a[:,:,1],b[:,:,1]) +
   ....:               np.multiply(a[:,:,2],b[:,:,2]))
   ....:
10000 loops, best of 3: 34.6 µs per loop

【讨论】：

@ askewchan 好答案。但是，我认为 3 个选项之间没有太大区别。确实，在我的情况下，einsum 也提供了最好的解决方案，但只有 5% 的改进。也许是因为以前我以这种方式计算向量？ a[:,:,0] = np.sin(azimA)*np.cos(elevA) a[:,:,1] = -np.cos(azimA)*np.cos(elevA) a[:,:,2] = np.sin(elevA) 来自 size=(1000,1000) 数组？
嗯，我实际上减少了尺寸以适应我的时间......可能差异不会那么大。从循环到乘法解决方案，您确实节省了大部分节省，但标准的做法是我展示的。
您可以通过仅计算 cos 和 sin 每个数组一次来节省您粘贴在评论中的时间：cosazimA = np.cos(azimA) 等，因为您要计算两次相同的东西。但这对以后的计算应该没有任何影响。
@askewchan 好观点！也可以与 einsum 进行交叉乘积，因为它是 c[0] = a[1]*b[2] - a[2]*b[1] c[1] = a[2]*b[0] - a[0]*b[2] c[2] = a[0]*b[1] - a[1]*b[0]

【解决方案2】：

Pythran 编译器可以通过以下方式进一步优化您的原始表达式：

删除临时数组
使用 SIMD 指令
使用多线程

如本例所示：

$ cat cross.py
#pythran export cross(float[][][], float[][][])
import numpy as np
def cross(a,b):
    return np.arccos(np.multiply(a[:, :, 0], b[:, :, 0]) + np.multiply(a[:, :, 1],b[:, :, 1]) + np.multiply(a[:, :, 2], b[:, :, 2]))
$ python -m timeit -s 'import numpy as np; a = np.random.uniform(low=-1.0, high=1.0, size=(1000, 1000, 3)); b = np.random.uniform(low=-1.0, high=1.0, size=(1000, 1000, 3)); c = np.zeros((1000, 1000)); from cross import cross' 'cross(a,b)'
10 loops, best of 3: 35.4 msec per loop
$ pythran cross.py -DUSE_BOOST_SIMD -fopenmp -march=native
$ python -m timeit -s 'import numpy as np; a = np.random.uniform(low=-1.0, high=1.0, size=(1000, 1000, 3)); b = np.random.uniform(low=-1.0, high=1.0, size=(1000, 1000, 3)); c = np.zeros((1000, 1000)); from cross import cross' 'cross(a,b)'
100 loops, best of 3: 11.8 msec per loop

【讨论】：