- html - 出于某种原因,IE8 对我的 Sass 文件中继承的 html5 CSS 不友好?
- JMeter 在响应断言中使用 span 标签的问题
- html - 在 :hover and :active? 上具有不同效果的 CSS 动画
- html - 相对于居中的 html 内容固定的 CSS 重复背景?
我最近试图比较不同的 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 次运行会引入一些差异,但我运行了更多次,结果大致如下:
这是一个合理的预期结果吗?为什么 C++ 库比 Numpy 慢?所有 3 个包都使用某种 Lapack/BLAS 后端,但 3 个包之间存在显着差异。特别是,Xtensor 会将我的 CPU 固定为 100% 使用 OpenBlas 的线程,但仍然设法获得更差的性能。
我想知道 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;
但是,时间并没有太大变化:
实际上,有点奇怪的是,添加 -O3 和 -ffast-math 实际上会稍微减慢代码速度。 -march=native 在我尝试时对我来说性能提升最大。此外,对于这些问题,OpenBLAS 比 CBLAS 快 2-3 倍。
最佳答案
首先,您计算的不是相同的东西。
要计算 l 矩阵的逆,对 Eigen 使用 l.inverse(),对 xtensor 使用 xt::linalg::inv()
当您将 Blas 链接到 Eigen 或 xtensor 时,这些操作会自动分派(dispatch)给您选择的 Blas。
我尝试替换反函数,用 MatrixXd 和 xt::xtensor 替换 auto 以避免延迟计算,将 openblas 链接到 Eigen、xtensor 和 numpy 并仅使用 -O3 标志编译,以下是我的 Macbook pro M1 上的结果:
Eigen-3.3.9(带 openblas) - ~ 38 毫秒
Eigen-3.3.9(无 openblas) - ~ 85 毫秒
xtensor-master(带 openblas) - ~41 毫秒
Numpy- 1.21.2(带 openblas) - ~35 毫秒。
关于Numpy vs Eigen vs Xtensor 线性代数基准奇数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/67463971/
我想了解表达式模板如何在 Eigen 中工作。 我知道两个动态双向量的总和是由如下所示的东西执行的: CwiseBinaryOp, VectorXd const, VectorXd const > o
在本征中,我们可以创建一个矩阵 Matrix3f m; m m(3, 8, 6); // Input after construction m.diagonal() << 3, 8, 6; Diag
我想知道是否有比使用 for 循环更简单的方法来解决我的问题。所以情况是这样的: 一般来说,我想从我的传感器收集数据点(消息类型为 Eigen::Vector3d,我无法更改它,因为它是一个巨大的框架
简而言之,问题是如何通过一个 Eigen::Map 反对一个期待一个函数 Eigen::MatrixXd 对象。 长话短说: 我有这个 C++ 函数声明 void npMatrix(const Eig
考虑以下代码。 const int N = 100; const float alpha = 0.9; Eigen::MatrixXf myVec = Eigen::MatrixXf::Random(
我试图获得两个张量的矩阵乘积,其中一个张量应该在相乘之前转置( At*B )。 到目前为止,我在 eigen documentation 中发现了什么是没有任何转置和两个矩阵都转置的矩阵乘积。 我正在
我的问题很简单,希望也有一个很好的答案:当我构造了 Eigen::MatrixXd 矩阵时,我可以使用多个线程同时填充矩阵中的行吗(如果我可以确保没有行被同时写入),或者我必须在每个线程中创建临时行对
#include #include namespace Eigen { // float op double -> double template struct ScalarBinaryOpTr
我在 Eigen 库的 API 上遇到了一些困难,即用于稀疏矩阵 Cholesky 分解的 SimplicialLLT 类。我需要分解三个矩阵,然后用它们来求解许多方程组(仅更改右侧) - 因此我只想
我正在尝试使用 Eigen 的 JacobiSVD 进行实验。特别是我试图从奇异值分解中重建输入矩阵。 http://eigen.tuxfamily.org/dox/classEigen_1_1Jac
我刚刚阅读了 Structures having static members Eigen 页。后者陈述如下: If you define a structure having members of
我正在使用 C++ 中的 Eigen 库来获取浮点方阵的平方根: MatrixXf gPrime(QUAD_EKF_NUM_STATES, QUAD_EKF_NUM_STATES); gPri
我正在尝试使用 vector 中的相应值来缩放矩阵中的所有列。如果此值为 0,我想用另一个按常数缩放的矩阵中的列替换该列。听起来很复杂,但在 Matlab 中它非常简单(但可能没有完全优化): a(:
使用 Eigen 将仿射变换转换为等距变换(即仅由旋转和平移组成)的最简单方法是什么? 图书馆? 两种变换都是 3D 的。仿射矩阵在左上象限有一个通用的 3x3 矩阵(即旋转、缩放和剪切),而等轴测在
我正在尝试生成 PCL 点云。我所有的点都在以下容器类型中: std::vector > 我想创建一个指向 PCL 点云的指针: pcl::PointCloud::Ptr pc 创建此点云的最有效方法
最近,我编译了 Eigen3 并用它在 Windows 上做一些线性代数任务。 但我想知道为什么 Eigen3 不需要链接额外的 lib 或 DLL(我只需要包含它的头文件) Eigen 是否在编译时
我正在使用特征张量库,在我看来,shuffle() 方法需要一个对象 Eigen::array 作为输入。 在我的实现中,我有一个 std::list 需要传递给 shuffle (当然我只在运行时知
我有两个问题。 Eigen::Affine3f 和 Eigen::Matrix4f 中表示的刚性变换是否相同? 我们如何从 Eigen::Affine3f 转换为 Eigen::Matrix4f,反之
我想用单个 3-D Eigen::Tensor 替换代码中的一系列矩阵。考虑到这一点,我尝试比较张量和矩阵的性能。 下面的函数“tensorContractTest”执行 (n,n,n) 3 阶张量与
如何将 ArrayXXd 数组转换为 MatrixXd?到目前为止,我已经完成了 MatrixXd temp_mat = my_array; 并且隐式转换似乎工作正常,但这是应该完成的方式吗?或者我应
我是一名优秀的程序员,十分优秀!