11

我想了解python的一种奇怪行为。让我们考虑一个M形状为 的矩阵6000 x 2000。这个矩阵用有符号整数填充。我想计算np.transpose(M)*M. 两种选择：

  • 当我“自然地”进行操作时（即不指定任何类型），numpy 会选择类型np.int32，操作大约需要 150 秒。
  • 当我强制类型为np.float64（使用dtype=...）时，相同的操作大约需要 2 秒。

我们如何解释这种行为？我天真地认为 int 乘法比 float 乘法便宜。

非常感谢你的帮助。

4

2 回答 2

7

不，整数乘法并不便宜。但稍后会详细介绍。最有可能（我有 99% 的把握）在毯子下numpy调用BLAS例程，其效率可以达到峰值 CPU 性能的 90%。矩阵乘法没有特殊规定int，很可能是在 Python 中完成的，而不是机器编译的版本——我实际上是错的，见下文。

关于int速度float：在大多数架构（英特尔）上，它们大致相同，每条指令大约 3-5 个周期左右，都有串行（X87）和矢量（XMM）版本。在桑迪桥上，PMUL***（整数向量乘法）是 5 个周期，MULP*（浮动乘法）也是如此。使用 Sandy Bridge，您还拥有 256 位 SIMD 向量操作 (YMM) -float每条指令可以获得 8 个操作 - 我不确定是否有int对应的。

这是一个很好的参考： http ://www.agner.org/optimize/instruction_tables.pdf

也就是说，指令延迟并不能解释 75 倍的速度差异。它可能是优化的 BLAS（可能是线程化的）和 int32 在 Python 而不是 C/Fortran 中处理的组合。

我分析了以下片段：

>>> F = (np.random.random((6000,2000))+4)
>>> I = F.astype(np.int32)
>>> np.dot(F, F.transpose()); np.dot(I, I.transpose())

这就是 oprofile 所说的：

CPU_CLK_UNHALT...|
  samples|      %|
------------------
  2076880 51.5705 multiarray.so
  1928787 47.8933 libblas.so.3.0

然而，libblas 是未优化的串行 Netlib Blas。使用良好的 BLAS 实现，47% 会低得多，尤其是在线程化的情况下。

编辑：似乎 numpy确实提供了整数矩阵乘法的编译版本。

于 2013-09-11T14:13:58.283 回答
5

我通过实验发现的一些补充信息。

这是可以规避的。计时在带有 intel mkl for BLAS 的 intel cpu 上。我还使用 fortran 有序数组来保持一切等价，scipy.linalg.blas即 fortran BLAS。

让我们看下面的例子：

from scipy.linalg.blas import sgemm
from scipy.linalg.blas import dgemm

arr_int64 = np.random.randint(-500,500,(6000,2000))

arr_int32 = array_int64.astype(np.int32)

arr_float64 = array_int64.astype(np.float64)+np.random.rand(6000,2000)

arr_float32 = array_int64.astype(np.float32)

首先让我们接听 DGEMM 电话。

%timeit np.dot(arr_float64.T,arr_float64) #400% CPU threaded BLAS
1 loops, best of 3: 969 ms per loop

%timeit np.dot(arr_float32.T,arr_float32) #400% CPU threaded BLAS
1 loops, best of 3: 513 ms per loop

%timeit np.dot(arr_int64.T,arr_int64)     #100% CPU?
1 loops, best of 3: 24.7 s per loop

%timeit np.dot(arr_int32.T,arr_int32)     #100% CPU?
1 loops, best of 3: 21.3 s per loop

直接调用 DGEMM/SGEMM：

%timeit dgemm(alpha=1, a=arr_float64, b=arr_float64, trans_a=True)
1 loops, best of 3: 1.13 s per loop

%timeit dgemm(alpha=1, a=arr_int64, b=arr_int64, trans_a=True)
1 loops, best of 3: 869 ms per loop

%timeit sgemm(alpha=1, a=arr_float32, b=arr_float32, trans_a=True)
1 loops, best of 3: 657 ms per loop

%timeit sgemm(alpha=1, a=arr_int32, b=arr_int32, trans_a=True)
1 loops, best of 3: 432 ms per loop

np.allclose( np.dot(arr_int32.T,arr_int32), sgemm(alpha=1, a=arr_int32, b=arr_int32, trans_a=True))
#True

通话中似乎发生了一些奇怪的事情np.dot。类似于朴素的算法速度：

%timeit np.einsum('ij,jk',arr_int32.T,arr_int32)
1 loops, best of 3: 14.1 s per loop

%timeit np.einsum('ij,jk',arr_int64.T,arr_int64)
1 loops, best of 3: 26 s per loop
于 2013-09-11T16:02:37.180 回答