最近我尝试比较不同的Python和C++矩阵库在线性代数性能方面的差异，以便确定在即将要开始的项目中使用哪一个或哪些。虽然有多种类型的线性代数运算，但我选择主要关注矩阵求逆，因为它似乎是给出奇怪结果的那个。我编写了下面的代码进行比较，但我认为我一定做错了什么。
C++ 代码
#include <iostream>
#include "eigen/Eigen/Dense"
#include <xtensor/xarray.hpp>
#include <xtensor/xio.hpp>
#include <xtensor/xview.hpp>
#include <xtensor/xrandom.hpp>
#include <xtensor-blas/xlinalg.hpp> //-lblas -llapack for cblas, -llapack -L OpenBLAS/OpenBLAS_Install/lib -l:libopenblas.a -pthread for openblas
//including accurate timer
#include <chrono>
//including vector array
#include <vector>
void basicMatrixComparisonEigen(std::vector<int> dims, int numrepeats = 1000);
void basicMatrixComparisonXtensor(std::vector<int> dims, int numrepeats = 1000);
int main()
{
std::vector<int> sizings{1, 10, 100, 1000, 10000, 100000};
basicMatrixComparisonEigen(sizings, 2);
basicMatrixComparisonXtensor(sizings,2);
return 0;
}
void basicMatrixComparisonEigen(std::vector<int> dims, int numrepeats)
{
std::chrono::high_resolution_clock::time_point t1;
std::chrono::high_resolution_clock::time_point t2;
using time = std::chrono::high_resolution_clock;
std::cout << "Timing Eigen: " << std::endl;
for (auto &dim : dims)
{
std::cout << "Scale Factor: " << dim << std::endl;
try
{
//Linear Operations
auto l = Eigen::MatrixXd::Random(dim, dim);
//Eigen Matrix inversion
t1 = time::now();
for (int i = 0; i < numrepeats; i++)
{
Eigen::MatrixXd pinv = l.completeOrthogonalDecomposition().pseudoInverse();
//note this does not come out to be identity. The inverse is wrong.
//std::cout<<l*pinv<<std::endl;
}
t2 = time::now();
std::cout << "Eigen Matrix inversion took: " << std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1).count() * 1000 / (double)numrepeats << " milliseconds." << std::endl;
std::cout << "\n\n\n";
}
catch (const std::exception &e)
{
std::cout << "Error: '" << e.what() << "'\n";
}
}
}
void basicMatrixComparisonXtensor(std::vector<int> dims, int numrepeats)
{
std::chrono::high_resolution_clock::time_point t1;
std::chrono::high_resolution_clock::time_point t2;
using time = std::chrono::high_resolution_clock;
std::cout << "Timing Xtensor: " << std::endl;
for (auto &dim : dims)
{
std::cout << "Scale Factor: " << dim << std::endl;
try
{
//Linear Operations
auto l = xt::random::randn<double>({dim, dim});
//Xtensor Matrix inversion
t1 = time::now();
for (int i = 0; i < numrepeats; i++)
{
auto inverse = xt::linalg::pinv(l);
//something is wrong here. The inverse is not actually the inverse when you multiply it out.
//std::cout << xt::linalg::dot(inverse,l) << std::endl;
}
t2 = time::now();
std::cout << "Xtensor Matrix inversion took: " << std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1).count() * 1000 / (double)numrepeats << " milliseconds." << std::endl;
std::cout << "\n\n\n";
}
catch (const std::exception &e)
{
std::cout << "Error: '" << e.what() << "'\n";
}
}
}
这是使用以下内容编译的：
g++ cpp_library.cpp -O2 -llapack -L OpenBLAS/OpenBLAS_Install/lib -l:libopenblas.a -pthread -march=native -o benchmark.exe
对于OpenBLAS，
g++ cpp_library.cpp -O2 -lblas -llapack -march=native -o benchmark.exe
适用于 cBLAS 的编译器为 g++ 版本 9.3.0。
Python 3 也同样需要。
import numpy as np
from datetime import datetime as dt
#import timeit
start=dt.now()
l=np.random.rand(1000,1000)
for i in range(2):
result=np.linalg.inv(l)
end=dt.now()
print("Completed in: "+str((end-start)/2))
#print(np.matmul(l,result))
#print(np.dot(l,result))
#Timeit also gives similar results
我将关注于在我的计算机上可以运行的最大规模，即1000x1000。虽然我知道只有2个运行会引入一些差异，但是我已经运行了更多的测试，并且结果与以下结果大致相同：
- Eigen 3.3.9: 196.804毫秒 - 使用OpenBlas的Xtensor/Xtensor-blas: 378.156毫秒 - Numpy 1.17.4: 172.582毫秒
这个结果合理吗？为什么C++库比Numpy慢？这三个软件包都使用某种Lapack/BLAS后端，但是它们之间存在显着差异。尤其是，Xtensor将我的CPU占用率锁定在100％并使用OpenBlas的线程，但性能仍然很差。
我想知道C++库是否实际上执行矩阵的逆/伪逆，以及这是否导致了这些结果。在C ++测试代码的注释部分中，我指出当我检查Eigen和Xtensor的结果时，矩阵和其逆矩阵之间的乘积根本不接近单位矩阵。我尝试使用较小的矩阵（10x10）进行测试，以为可能是精度误差，但问题仍然存在。在另一个测试中，我测试了矩阵的秩，这些矩阵的秩都是满秩的。为了确认我没有疯狂，我尝试使用inv（）而不是pinv（），两种情况的结果都相同。我是在为这个线性代数基准测试使用错误的函数，还是Numpy正在对两个失灵的低级库挑衅呢？
编辑： 感谢大家对这个问题的关注。 我想我已经找到了问题所在。我怀疑Eigen和Xtensor具有惰性评估，这实际上会在下游产生错误，并输出随机矩阵而不是逆矩阵。我能够通过以下代码替换来修复奇怪的数值反转失败：
auto temp = Eigen::MatrixXd::Random(dim, dim);
Eigen::MatrixXd l(dim,dim);
l=temp;
并且
auto temp = xt::random::randn<double>({dim, dim});
xt::xarray<double> l =temp;
然而，时间并没有改变太多：
- Eigen 3.3.9: 201.386 毫秒
- Xtensor/Xtensor-blas w/ OpenBlas: 337.299 毫秒。
- Numpy 1.17.4: (之前的数据) 172.582 毫秒
实际上有点奇怪的是，加上-O3和-ffast-math选项竟然使代码变慢了一点。 当我尝试时，-march=native带来了最大的性能提升。 另外，在这些问题上，OpenBLAS比CBLAS快2-3倍。
-O2编译选项中未启用自动向量化，而在-O3中启用。为了在稍微降低二进制文件的可移植性的代价下获得小幅加速，您可以使用-march=native。此外，如果您不关心非常精确的结果，并且确定NaN/Inf/等不应出现，则可以使用-ffast-math。这应该会使您的 C++ 代码稍微快一些。除此之外，三个实现所使用的算法不一定相同，某些算法可能会根据 输入矩阵 更快或更慢。它们也可能更或少数值稳定。 - Jérôme Richard