Numpy vs Eigen vs Xtensor线性代数基准奇数

Question

我最近试图比较不同的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 << "nnn";
}
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 << "nnn";
}
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次运行会带来一些差异，但我已经运行了更多，结果大致如下：

• 特征3.3.9:196.804毫秒
• Xtensor/Xtensor blas w/OpenBlas:378.156毫秒
• Numpy 1.17.4:172.582毫秒

这是一个合理的结果吗？为什么C++库比Numpy慢？所有3个包都使用某种Lapack/BLAS后端，但这3个包之间存在显著差异。特别是，Xtensor会将我的CPU与OpenBlas的线程绑定到100%的使用率，但性能仍然较差。

我想知道C++库是否真的在执行矩阵的逆/伪逆，以及这是否是导致这些结果的原因。在C++测试代码的评论部分，我注意到，当我健全地检查Eigen和X张量的结果时，矩阵及其逆矩阵之间的矩阵乘积甚至不接近单位矩阵。我尝试使用较小的矩阵(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;

然而，时间没有太大变化：

• 特征3.3.9:201.386毫秒
• Xtensor/Xtensor blas w/OpenBlas：337.299毫秒
• Numpy 1.17.4:(从前(172.582毫秒

实际上，有点奇怪的是，添加-O3和-fast数学实际上使代码慢了一点-march=native在我尝试它时对我来说性能提升最大。此外，OpenBLAS在这些问题上比CBLAS快2-3倍。

Answer 1

首先，您不是在计算相同的东西。

要计算l矩阵的逆，请对Eigen使用l.inverse((，对extensor 使用xt:：linalg:：inv((

当您将Blas链接到Eigen或extensor时，这些操作会自动分派到您选择的Blas。

我尝试替换反函数，用MatrixXd和xt:：xtensor替换auto以避免延迟求值，将openblas链接到Eigen、xtensor和numpy，并仅使用-O3标志编译，以下是我的Macbook pro M1:上的结果

特征-3.3.9(带openblas(-~38ms

特征-3.3.9(无openblas(-~85 ms

扩展器主控器(带openblas(-~41毫秒

Numpy-1.21.2(带openblas(-约35毫秒