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

【问题标题】：Why is the GNU scientific library matrix multiplication slower than numpy.matmul?为什么 GNU 科学库矩阵乘法比 numpy.matmul 慢？
【发布时间】：2021-08-05 11:53:45
【问题描述】：

为什么用 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, 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 维的积分，我正在使用 Monte-Carlo 方法进行评估。为此，我几乎将整个被积函数编写为 Numpy 数组操作，它的工作速度非常快，但我需要它更快，以便将相同的被积函数评估 100.000 到 500.000 次，因此任何一点改进都会有所帮助。用 C/C++ 编写相同的代码有意义还是我应该坚持使用 Numpy？谢谢！

【问题讨论】：

标签： python c++ performance numpy gsl

【解决方案1】：

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 性能有多差的一种简单方法是使用 profiler（如 Linux 上的 perf 或 Windows 上的 VTune）。在您的情况下，Linux 性能，报告 > 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 instructions。

【讨论】：