如何解决Numpy vs Eigen vs Xtensor 线性代数基准奇数
我最近试图比较不同的 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 毫秒
- Xtensor/Xtensor-blas 和 OpenBlas:378.156 毫秒
- Numpy 1.17.4:172.582 毫秒
这是一个合理的预期结果吗?为什么 C++ 库比 Numpy 慢?所有 3 个包都使用某种 Lapack/BLAS 后端,但 3 个包之间存在显着差异。特别是,Xtensor 将使用 OpenBlas 线程将我的 cpu 固定到 100% 使用率,但仍然设法获得更差的性能。
我想知道 C++ 库是否真的在执行矩阵的逆/伪逆,以及这是否是导致这些结果的原因。在 C++ 测试代码的注释部分中,我注意到当我对 Eigen 和 Xtensor 的结果进行完整性检查时,矩阵与其逆矩阵之间的结果矩阵乘积甚至不接近单位矩阵。我尝试使用较小的矩阵(10x10),认为这可能是一个精度错误,但问题仍然存在。在另一个测试中,我测试秩,这些矩阵是满秩的。为了确保我没有发疯,我在两种情况下都尝试使用 inv() 而不是 pinv() ,结果是一样的。我是否在这个线性代数基准测试中使用了错误的函数,或者这个 Numpy 是否在 2 个功能失调的低级库上动摇了刀?
编辑: 感谢大家对这个问题的关注。我想我已经弄清楚了这个问题。我怀疑 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 和 OpenBlas:337.299 毫秒。
- Numpy 1.17.4:(之前)172.582 毫秒
实际上,有点奇怪的是,添加 -O3 和 -ffast-math 实际上稍微减慢了代码的速度。当我尝试时,-march=native 对我来说性能提升最大。此外,对于这些问题,OpenBLAS 比 CBLAS 快 2-3 倍。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。