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

给出 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 维的积分，我正在使用蒙特卡洛方法进行计算。为此，我几乎将整个被积函数写为 Numpy 数组操作，这工作得非常快，但我需要它更快，以便评估相同的被积函数 100.000 到 500.000 次，所以任何小的改进都会有所帮助。在 C/C++ 中编写相同的代码是否有意义，还是我应该坚持使用 Numpy？谢谢！

Answer 1

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

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

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

您可以找到 Numpy 使用的 BLAS 实现。在我的 Linux 机器上，Numpy 默认使用内部使用 OpenBLAS 的 CBLAS（奇怪的是，Numpy 没有直接检测到 OpenBLAS）。

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

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

更多信息：

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

250:   movss   (%rdx),%xmm1
       add     [=10=]x4,%rax
       add     [=10=]x4,%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          [=11=]x60,%rsi 
       mov          %r12,%rax 
       sar          [=11=]x3,%rax 
       cmp          [=11=]x2,%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。

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

Why is the GNU scientific library matrix multiplication slower than numpy.matmul?

c++

python

performance

numpy

gsl