为什么GNU科学库矩阵乘法比numpy.matmul慢？
为什么用 Numpy 的矩阵乘法比用gsl_blas_sgemmGSL快得多，例如：
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？谢谢！
回答
TL;DR： C++ 代码和 Numpy 不使用相同的矩阵乘法库。
GSL 库的矩阵乘法未优化。在我的机器上，它按顺序运行，不使用SIMD 指令（SSE / AVX），不能有效地展开循环以执行寄存器平铺。我还怀疑由于缺乏平铺，它也没有有效地使用 CPU 缓存。这些优化对于实现高性能和广泛用于快速线性代数库至关重要。
Numpy 使用安装在您机器上的BLAS 库。在许多 Linux 平台上，它使用 OpenBLAS 或 Intel MKL。两者都非常快（它们使用上述所有方法）并且应该并行运行。
您可以在此处找到 Numpy 使用的 BLAS 实现。在我的 Linux 机器上，Numpy 默认使用内部使用 OpenBLAS 的 CBLAS（奇怪的是，Numpy 没有直接检测到 OpenBLAS）。
有许多快速并行 BLAS 实现（GotoBLAS、ATLAS、BLIS 等）。开源 BLIS 库很棒，因为它的矩阵乘法在许多不同的架构上都非常快。
因此，改进 C++ 代码的最简单方法是使用cblas_sgemmCBLAS 函数并链接快速 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 指令。