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

如何解决为什么 GNU 科学库矩阵乘法比 numpy.matmul 慢？

为什么用 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)

给出 0.017 - 0.019 秒之间的值，而在 C++ 中：

#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,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？谢谢！

解决方法

TL;DR：C++ 代码和 Numpy 不使用相同的矩阵乘法库。

GSL 库的矩阵乘法没有优化。在我的机器上，它按顺序运行，不使用 SIMD 指令 (SSE/AVX)，不能有效地展开循环以执行寄存器平铺。我还怀疑由于缺乏平铺，它也没有有效地使用 CPU 缓存。这些优化对于实现高性能和广泛应用于快速线性代数库至关重要。

Numpy 使用安装在您机器上的 BLAS library。在许多 Linux 平台上，它使用 OpenBLAS 或 Intel MKL。两者都非常快（它们使用上述所有方法）并且应该并行运行。

您可以找到 Numpy here 使用的 BLAS 实现。在我的 Linux 机器上，Numpy 默认使用内部使用 OpenBLAS 的 CBLAS（奇怪的是，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,%ymm5 
       vmulpd       %ymm0,%ymm9 
       [...]

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