在C++中计算大规模稀疏矩阵特征值有没有更快的实现方法？

优化方案 #

1 替换稀疏Cholesky分解后端 #

你当前用的Spectra默认SparseCholesky底层调用Eigen自带的SimplicialLLT，这是通用稀疏LLT实现，针对大规模稀疏矩阵性能很低。质量矩阵M是对称正定矩阵，完全可以用更高性能的稀疏分解实现：

• 如果你安装了Intel MKL，定义EIGEN_USE_MKL_ALL编译宏，替换为MKL PARDISO作为分解后端，分解3万阶稀疏矩阵的速度可以提升5~10倍。
• 由于是一维杆导出的带状矩阵（带宽极小），直接用Eigen的带状矩阵分解接口，不要用通用稀疏矩阵接口，开销能降低一个数量级。

2 调整Arnoldi子空间参数 #

你当前设置的子空间维度是500，对于仅求解10个特征值的场景完全过剩，子空间过大会大幅提升每次迭代的正交化、矩阵乘开销。常规场景下子空间维度取所求特征值数量的2_{5倍即可，也就是20}50，特征值聚集严重的场景最多调到100，仅此一项就能让求解速度提升数倍。

3 开启编译与运算优化 #

• 编译时开启最高等级优化 -O3 -march=native，同时开启OpenMP支持，Eigen和Spectra的矩阵运算、正交化步骤都能利用多线程加速。
• 启用Eigen的向量指令优化，确保没有定义EIGEN_DONT_VECTORIZE这类禁用向量优化的宏。

4 超大规模场景更换求解库 #

如果要处理百万阶及以上的矩阵，Spectra作为单机优先的库性能上限很低，建议换成SLEPc库，这是专门面向大规模稀疏特征值问题的并行求解库，支持MPI分布式并行，针对广义特征值、极小特征值求解做了深度优化，可支撑百万到千万阶矩阵的求解需求。

修改后的代码示例：

// 编译前定义EIGEN_USE_MKL_ALL宏启用MKL加速
#include <Eigen/Sparse>
#include <Spectra/SymGEigsSolver.h>
#include <Spectra/MatOp/SparseSymMatProd.h>
#include <Spectra/MatOp/SparseCholesky.h>

using namespace Eigen;
using namespace Spectra;

// ... 生成sparsek、sparsem的逻辑 ...

SparseSymMatProd<double> eig1(sparsek);
SparseCholesky<double, Lower, MKL_PARDISO> eig2(sparsem); // 用MKL PARDISO作为分解后端

// 子空间维度从500改为30，仅求10个特征值足够
SymGEigsSolver <SparseSymMatProd<double>, SparseCholesky<double, Lower, MKL_PARDISO>, GEigsMode::Cholesky> eigs(eig1, eig2, 10, 30);

eigs.init();
int nconv = eigs.compute(SortRule::SmallestMagn);

VectorXd evalues;
MatrixXd evecs;
if (eigs.info() == CompInfo::Successful)
{
    evalues = eigs.eigenvalues();
    evecs = eigs.eigenvectors();
}
std::cout << "Generalized Eigenvalues found: \n" << evalues << std::endl;
std::cout << "Generalized Eigenvectors found: \n" << evecs.bottomRows(10) << std::endl;

内容的提问来源于stack exchange，提问作者Mohammad Ali Yaqteen