为什么GNU科学库的矩阵乘法比numpy.matmul慢？

Question

为什么用Numpy进行矩阵乘法要比GSL中的gsl_blas_sgemm快很多，例如：

import numpy as np
import time 


N = 1000
M = np.zeros(shape=(N, N), dtype=np.float)

for i in range(N):
    for j in range(N):
        M[i, j] = 0.23 + 100*i + j

tic = time.time()
np.matmul(M, M)
toc = time.time()
print(toc - tic)

在C++中，这个过程需要约0.017 - 0.019秒的时间：

#include <chrono>
#include <iostream>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

using namespace std::chrono;

int main(void) {

    int N = 1000;

    gsl_matrix_float* M = gsl_matrix_float_alloc(N, N);
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            gsl_matrix_float_set(M, i, j, 0.23 + 100 * i + j);
        }
    }

    gsl_matrix_float* C = gsl_matrix_float_alloc(N, N); // save the result into C

    auto start = high_resolution_clock::now();

    gsl_blas_sgemm(CblasNoTrans, CblasNoTrans, 1.0, M, M, 0.0, C);

    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<milliseconds>(stop - start);
    std::cout << duration.count() << std::endl;

    return 0;
}

我得到了大约2.7秒的乘法运行时间。我还使用最大速度选项/02进行编译。我正在使用Visual Studio工作。我一定做错了什么。我并没有期望从C ++代码中获得更好的性能，因为我知道Numpy是优化的C代码，但我也没有预料到它比Python慢约150倍。为什么会这样？如何相对于Numpy改善乘法运行时间？

问题背景： 我需要使用蒙特卡罗方法评估1000到2000个维度的积分。为此，我几乎将整个被积函数写成Numpy数组操作，这很快，但我需要更快的速度才能计算相同的被积函数100,000至500,000次，所以任何小的改进都会有所帮助。是否值得用C / C ++编写相同的代码，还是应该坚持使用Numpy？谢谢！

- ibroketheinternet

Answer 1

\{\{简短概括：\}\} C++代码和Numpy使用的矩阵乘法库不一样。

\{\{GSL库的矩阵乘法未经过优化\}\}，在我的计算机上以串行方式运行，没有使用SSE/AVX指令，没有有效地展开循环以执行寄存器分块，我还怀疑由于缺乏分块而无法有效地使用CPU缓存。这些优化对于实现高性能至关重要，并且在快速线性代数库中广泛使用。 Numpy使用安装在计算机上的BLAS库进行矩阵乘法 。在许多Linux平台上，它使用OpenBLAS或Intel MKL。两者都非常快（它们使用上述所有方法）并且应该并行运行。

您可以在此处找到Numpy使用的BLAS实现。在我的Linux机器上，默认情况下，Numpy使用CBLAS，其内部使用OpenBLAS（奇怪的是Numpy不能直接检测到OpenBLAS）。

有许多快速的并行BLAS实现（GotoBLAS、ATLAS、BLIS等）。开源BLIS库非常棒，因为它的矩阵乘法在许多不同体系结构上运行得非常快。

因此，改进C++代码的最简单方法是使用cblas_sgemm CBLAS函数，并链接快速的BLAS库，例如OpenBLAS或BLIS等。

---

更多信息：

了解GSL性能的一种简单方法是使用分析器（如Linux上的perf或Windows上的VTune）。在您的情况下，Linux perf报告超过99％的时间花费在libgslcblas.so（即GSL库）中。更具体地说，大部分执行时间都花费在以下汇编循环中：

250:   movss   (%rdx),%xmm1
       add     $0x4,%rax
       add     $0x4,%rdx
       mulss   %xmm2,%xmm1           # scalar instructions
       addss   -0x4(%rax),%xmm1
       movss   %xmm1,-0x4(%rax)
       cmp     %rax,%r9
     ↑ jne     250

关于 Numpy，它的 99% 时间都花在 libopenblasp-r0.3.13.so（即 OpenBLAS 库）上。更具体地说，是在函数 dgemm_kernel_HASWELL 的以下汇编代码中：

110:   lea          0x80(%rsp),%rsi 
       add          $0x60,%rsi 
       mov          %r12,%rax 
       sar          $0x3,%rax 
       cmp          $0x2,%rax 
     ↓ jl           d26 
       prefetcht0   0x200(%rdi)          # Data prefetching
       vmovups      -0x60(%rsi),%ymm1 
       prefetcht0   0xa0(%rsi)
       vbroadcastsd -0x80(%rdi),%ymm0    # Fast SIMD instruction (AVX)
       prefetcht0   0xe0(%rsi)
       vmovups      -0x40(%rsi),%ymm2 
       prefetcht0   0x120(%rsi)
       vmovups      -0x20(%rsi),%ymm3 
       vmulpd       %ymm0,%ymm1,%ymm4
       prefetcht0   0x160(%rsi)
       vmulpd       %ymm0,%ymm2,%ymm8 
       vmulpd       %ymm0,%ymm3,%ymm12 
       prefetcht0   0x1a0(%rsi)
       vbroadcastsd -0x78(%rdi),%ymm0 
       vmulpd       %ymm0,%ymm1,%ymm5 
       vmulpd       %ymm0,%ymm2,%ymm9 
       [...]

我们可以清楚地看到，GSL代码没有被优化（因为标量代码和简单循环），而OpenBLAS代码是被优化过的，因为它至少使用了宽SIMD指令、数据预取和循环展开。请注意，执行的OpenBLAS代码并不是最优的，因为它可以使用我处理器上可用的FMA指令。

- Jérôme Richard