我现在需要求解一类含参哈密顿矩阵的特征值，矩阵形式是对角矩阵mat1加上参数g乘以另一个矩阵mat2：

matrix = mat1 + g*mat2

我的目标是得到随参数g变化的最小的n个特征值（不是绝对值最小，是数值最小），最后画出横轴为g、纵轴为特征值的曲线图。我用了Eigen和Spectra库的几种方法，遇到了不同的问题，想请教怎么快速准确地得到结果，解决特征值曲线交叉时的错位和检测缺失问题。

方法1：Eigen的SelfAdjointEigenSolver（稠密矩阵） #

这个方法的结果是最好的，因为它会计算所有特征值，我总能轻松提取到最小的那些。唯一的小问题是特征值曲线交叉时，因为Eigen会自动排序特征值，导致曲线看起来“交换了位置”，不过这个应该是排序的正常现象，我觉得没法避免。但问题在于，稠密矩阵计算所有特征值的速度太慢了，矩阵规模大一点就跑不动了。

下面是用这个方法得到的特征值曲线图：

方法2：Spectra的SymEigsShiftSolver #

这个方法是找最接近设定偏移值的特征值，但这就带来了问题：当某些特征值变得非常负时，曲线会出现错位。我试过去把偏移值设为负数，勉强能用，但这要求我提前知道最小特征值的大概范围（而我其实不知道），而且对上层的特征值会产生奇怪的异常。

下面是用这个方法得到的曲线图：

方法3：Spectra的GenEigsSolver（设置SortRule::SmallestReal） #

我本来想通过这个排序规则直接提取最小的实特征值，但结果更糟：不仅特征值曲线会错位，还有些特征值直接没检测到，完全没法得到连续的能级曲线。

下面是这个方法的结果图，以及其中一次特征值错位的细节：

现在的核心问题是：有没有办法既快速（不用算所有特征值）又准确地得到最小的n个特征值，同时解决这些曲线错位、特征值漏检的问题？单纯排序特征值没用，因为有时候某个能级直接没被检测到，结果完全对不上。

为了复现问题，我用的测试矩阵很简单：mat1是对角元从0到100的对角矩阵，mat2是对称随机矩阵（元素范围-1到1）。以下是生成这些图的测试代码：

#include<iostream>
#include <Eigen/Sparse>
#include <Spectra/SymEigsShiftSolver.h>
#include <Spectra/MatOp/SparseSymShiftSolve.h>
#include <chrono>
#include <cstdlib> 
#include<fstream>
#include <Eigen/Dense>
#include <Spectra/GenEigsSolver.h>
#include <Spectra/MatOp/SparseGenMatProd.h>
using Eigen::MatrixXd;
using Eigen::SelfAdjointEigenSolver;
using namespace Spectra;
using namespace Eigen;
using namespace std;
int main()
{
    int N=100;
    SparseMatrix<double>mat1(N,N);
    SparseMatrix<double>mat2(N,N);
    for(int i=0;i<N;i++)
    {
        mat1.insert(i,i)=i;
        for(int j=i;j<N;j++)
        {
            double r=rand()/double(RAND_MAX)*2.0-1;
            mat2.insert(i,j)=r;
            if(i!=j)mat2.insert(j,i)=r;
        }
    }

    cout<<mat1<<"\n\n"<<mat2<<"\n";
    
    // Calculating ALL Eigenvalues from dense matrix using SelfAdjointEigenSolver
    ofstream Energies("output1.csv");
    SelfAdjointEigenSolver<MatrixXd> es;
    for(double g=-10;g<10;g+=0.1)
    {cout<<g<<"    \n";
        MatrixXd matrix=mat1+g*mat2;
        es.compute(matrix);
        Energies<<g<<"\t";
        for(int i=0;i<es.eigenvalues().rows();i++)
        {
            Energies<<es.eigenvalues()[i]<<"\t";
        }
        Energies<<"\n";
        Energies.flush();
    }
    Energies.close();
    
    // Calculating first 20 Eigenvalues from sparse matrix using SymEigsShiftSolver
    Energies.open("output2.csv");
    for(double g=-10;g<10;g+=0.1)
    {
        cout<<g<<"    \n";
        
        SparseSymShiftSolve<double> op(mat1+g*mat2);
        SymEigsShiftSolver<SparseSymShiftSolve<double>> eig(op,20,50,0.0);
        eig.init();
        cout<<eig.compute()<<"\n";
        VectorXd evalues;
        if(eig.info() == CompInfo::Successful)
            evalues = eig.eigenvalues();
        int N=evalues.rows();
        Energies<<g<<"\t";
        for(int i=0;i<N;i++)
        {
            Energies<<evalues(i)<<"\t";
        }
        Energies<<"\n";
        Energies.flush();
    }
    Energies.close();
    
    // Calculating first 20 Eigenvalues from dense matrix using GenEigSolver
    Energies.open("output3.csv");
    
    for(double g=-10;g<10;g+=0.1)
    {
        cout<<g<<"    \n";
        SparseMatrix<double> matrix=mat1+g*mat2;
        SparseGenMatProd<double> op(matrix);
        GenEigsSolver<SparseGenMatProd<double>> eig(op, 20, 50);
        eig.init();
        cout<<eig.compute(SortRule::SmallestReal )<<"\n";
        VectorXd evalues;
        if(eig.info() == CompInfo::Successful)
            evalues = eig.eigenvalues().real();
        int N=evalues.rows();
        Energies<<g<<"\t";
        for(int i=0;i<N;i++)
        {
            Energies<<evalues(N-i-1)<<" ";
        }
        Energies<<"\n";
        Energies.flush();
    }
    Energies.close();
}

备注：内容来源于stack exchange，提问作者martes