为什么用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?谢谢!