C++ 中的 Eigen 比 Python 中的 Scipy 慢

Question

我需要快速生成巨型矩阵并对角化。为了易于使用，我开始使用 Python，但除了代数之外，它做的许多事情在 Python 中都非常慢。 在python中，我使用 scipy.linalg.eigh 来计算一些矩阵的特征值和特征向量（这里使用的矩阵是465x465，矩阵有很多零（大约一半），但大部分不是稀疏的 - 然而我一直在检查两者稀疏和密集方法。编辑：我已更改代码以使用完全随机矩阵，因此不需要创建自己的文件，并且仍然显示结果。） - 特征值计算非常快。这是加载文件并计算特征值的代码。

import numpy as np
from scipy import linalg as ln
import time


data=np.random.rand(465,465)
H0=data
data=data=np.random.rand(465,465)
HI=data

x=np.arange(-10,10)
t=time.time()
for g in x:
    t0=time.time()
    x,v=ln.eigh(H0+g*HI)
    print(time.time()-t0)
print(time.time()-t)

每次对角化的时间在0.02s左右。总时间约为0.4秒。

然后我尝试在 C++ 中使用“Eigen”库：

#include <iostream>
#include <fstream>
#include <Eigen/Dense>
#include <time.h>
using Eigen::MatrixXd;
using Eigen::SelfAdjointEigenSolver;
using namespace std;
int main(int argc, char **argv)
{
    int Nstates=465;
    MatrixXd X=MatrixXd::Random(Nstates,Nstates);
    MatrixXd H0=X+X.transpose();
    X=MatrixXd::Random(Nstates,Nstates);
    MatrixXd HI=X+X.transpose();

    int t=clock();
    SelfAdjointEigenSolver<MatrixXd> es;
    for(double g=-10;g<10;g+=1)
    {
        int t0=clock();
        
        MatrixXd matrix=H0+g*HI;
        es.compute(matrix);
        cout<<double(clock()-t0)/double(CLOCKS_PER_SEC)<<"\n";
    }
    cout<<double(clock()-t)/double(CLOCKS_PER_SEC)<<"\n";
    
    return 0;
}

虽然代码有效，但每次对角化大约需要 3.4 秒 整个考验过程花了65.5秒

然后，我尝试使用 Spectra 来完成此操作，据我了解，这应该有助于计算巨型矩阵。在此代码中，我使用密集矩阵，但我尝试使用稀疏矩阵，时间几乎相同。在光谱中，我只计算前 30 个特征值，并且我尝试了很多参数 - 最快的是第二个参数大约 80 或 90，而对于随机矩阵，大约需要 0.7 秒才能完成一个，15.5 秒才能完成对角化所有。

#include <iostream>
#include <fstream>
#include <Eigen/Core>
#include <Spectra/SymEigsSolver.h>
#include <time.h>
using Eigen::MatrixXd;
using namespace std;
using namespace Spectra;
int main(int argc, char **argv)
{
    int Nstates=465;
    MatrixXd X=MatrixXd::Random(Nstates,Nstates);
    MatrixXd H0=X+X.transpose();
    X=MatrixXd::Random(Nstates,Nstates);
    MatrixXd HI=X+X.transpose();
    int t=clock();
    for(double g=-10;g<10;g+=1)
    {
        int t0=clock();
        
        MatrixXd matrix=H0+g*HI;
        DenseSymMatProd<double> op(matrix);
        SymEigsSolver<DenseSymMatProd<double>> eigs(op, 30, 80);
        eigs.init();
        int nconv = eigs.compute(SortRule::LargestAlge);
    
        cout<<double(clock()-t0)/double(CLOCKS_PER_SEC)<<"\n";
    }
    cout<<double(clock()-t)/double(CLOCKS_PER_SEC)<<"\n";
    
    return 0;
}

我知道 scipy 使用 Lapack 来解决，我真的不想在 C++ 中使用 Lapack - 但正如我读到的，Eigen 应该在速度上与 Lapack 竞争，并在易用性方面阻碍它。 Spectra 原本应该在速度上击败 Eigen。所以我相信我做错了什么。我将非常感谢您的帮助，并建议我应该做什么。

编辑：最初在 C++ 中我使用这个编译器标志：-gdwarf-2;-O0;-Wall

我尝试将 -O0 更改为 -O2，这最终使一切变得更快 - Eigen 使用了 0.3s，而 Spectra 使用了 0.04（这仍然比 Scipy 慢，但我猜是有效的）。

Edit2：我用 chrono::steady_clock 而不是时钟测试了代码，但时间几乎相同。

Edit3：使用 -O3 标志，甚至 Eigen 也接近 Scipy 速度，而 Spectra 是 Scipy 的两倍。所以我觉得还是不错的。

Answer 1

非常感谢，我已经尝试了一些东西，事实上，添加“-O3”标志进行优化，将计算时间从每个矩阵约 3.4 秒降低到每个矩阵 0.02 秒 - 我从未见过该标志带来如此巨大的变化，所以从来没有费心去使用它，因为正如我听说的那样，它使调试变得更加困难。最后，使用 -O3 的代码大约是 python 代码的两倍，这正是我所期望的。谢谢！