numpy比Eigen C++更快、更高效

Question

最近，我和一位同事讨论了在性能方面比较python和C++的问题。我们两个大部分时间都在使用这些语言进行线性代数。所以我写了两个脚本，一个在python3中使用numpy，另一个在C++中使用Eigen。

Python3 numpy版本matmul_numpy.py:

import numpy as np
import time
a=np.random.rand(2000,2000)
b=np.random.rand(2000,2000)
start=time.time()
c=a*b
end=time.time()
print(end-start) 

如果我用运行这个脚本

python3 matmul_numpy.py

这将返回：

0.07 seconds

C++特征版本matmul_eigen.cpp:

#include <iostream>
#include <Eigen/Dense>
#include "time.h"
int main(){
clock_t start,end;
size_t n=2000;
Eigen::MatrixXd a=Eigen::MatrixXd::Random(n,n);
Eigen::MatrixXd b=Eigen::MatrixXd::Random(n,n);
start=clock();
Eigen::MatrixXd c=a*b;
end=clock();
std::cout<<(double)(end-start)/CLOCKS_PER_SEC<<std::endl;
return 0;}

我编译它的方式是，

g++ matmul_eigen.cpp -I/usr/include/eigen3 -O3 -march=native -std=c++17 -o matmul_eigen

这将返回(c++11和c++17)：

0.35 seconds

这对我来说很奇怪，1-为什么这里的numpy比C++快？我是否缺少任何其他用于优化的标志？

我想也许是因为python解释器，它在这里执行程序的速度更快。所以我用cython在堆栈中使用这个线程编译代码。

编译后的python脚本仍然更快(0.11秒)。这又给我增加了两个问题：

2-为什么时间变长了？解释器还会进行优化吗？

3-为什么python脚本的二进制文件(37kb)比c++(57kb)小？

如果有任何帮助，我将不胜感激，

感谢

Answer 1

最大的问题是您正在比较两个完全不同的东西：

• 在Numpy中，a*b执行逐元素乘法，因为a和b是2D阵列，不被视为矩阵。CCD_ 4执行矩阵乘法
• 在Eigen中，a*b执行矩阵乘法，而不是逐元素乘法(请参阅文档)。这是因为a和b是矩阵，而不仅仅是2D阵列

两者给出了完全不同的结果。此外，矩阵乘法在O(n**3)时间内进行，而元素乘法在O(n**2)时间内进行。矩阵乘法核通常是高度优化的，并且有计算限制。它们通常被大多数BLAS库并行化。元素乘法是内存绑定的(尤其是在页面错误的情况下)。因此，这并不奇怪，矩阵乘法比逐元素乘法慢，并且由于后面是内存绑定的，所以差距也不大。

在我的i5-9600KF处理器(有6个内核)上，Numpy执行a*b(按顺序)需要9毫秒，执行a@b(并行，使用OpenBLAS)需要65毫秒。

注意，像这样的Numpy元素乘法不是并行的(至少在Numpy的标准默认实现中不是这样)。Numpy的矩阵乘法使用一个BLAS库，默认情况下通常是OpenBLAS(这取决于实际的目标平台)。Eigen也应该使用BLAS库，但它可能与Numpy的库不同。

还要注意，clock不是测量并行代码的好方法，因为它测量的不是墙上的时钟时间，而是CPU时间(更多信息请参阅本文)。CCD_ 13通常是C++中更好的替代方案。

---

3-为什么python脚本的二进制文件(37kb)比c++(57kb)小？

Python通常编译为字节码，而字节码不是本机汇编代码。C++通常被编译成可执行程序，其中包含汇编的二进制代码、用于运行程序的附加信息以及元信息。字节码通常非常紧凑，因为它们级别更高。本机编译器可以执行优化，使程序变得更大，例如循环展开和内联。这样的优化不是由CPython(默认的Python解释器)完成的。事实上，CPython对字节码不执行(高级)优化。请注意，您可以告诉GCC等本地编译器使用-Os和-s等标志生成较小的代码(尽管通常较慢)。

Answer 2

根据我从@Jérôme Richard的回复和评论@user17732522中学到的内容。我在比较中似乎犯了两个错误，

1-我在python脚本中定义乘法时犯了一个错误，它应该是np.matmul(a，b)或np.dot(a，b)或a@b。而不是a*b，这是一个元素乘法。2-我在C++代码中没有正确测量时间。clock_t不适合此计算，std:：chrono：：steady_clock效果更好。

通过应用这些注释，c++特征代码的速度是python的10倍。

matmul_eigen.cpp:的更新代码

#include <iostream>
#include <Eigen/Dense>
#include <chrono>
int main(){
size_t n=2000;
Eigen::MatrixXd a=Eigen::MatrixXd::Random(n,n);
Eigen::MatrixXd b=Eigen::MatrixXd::Random(n,n);
auto t1=std::chrono::steady_clock::now();
Eigen::MatrixXd c=a*b;
auto t2=std::chrono::steady_clock::now();
std::cout<<(double)std::chrono::duration_cast<std::chrono::microseconds>(t2-t1).count()/1000000.0f<<std::endl;
return 0;}

要进行编译，应同时考虑矢量化和多线程标志。

g++ matmul_eigen.cpp -I/usr/include/eigen3 -O3 -std=c++17 -march=native -fopenmp -o eigen_matmul

要使用多个线程来运行代码：

OMP_NUM_THREADS=4 ./eigen_matmul

其中"；4〃；是openmp可以使用的CPU数量，您可以使用查看您有多少

lscpu | grep "CPU(s):"

这将返回0.104秒。

更新的python脚本matmul_numpy.py:

import numpy as np
import time
a=np.random.rand(2000,2000)
b=np.random.rand(2000,2000)
a=np.array(a, dtype=np.float64)
b=np.array(b, dtype=np.float64)
start=time.time()
c=np.dot(a,b)
end=time.time()
print(end-start)

要运行代码，

python3 matmul_numpy.py

这将返回1.0531秒。

关于这样的原因，我认为@Jérôme Richard的回应是一个更好的参考。