求助：使用Spectra GenEigsSolver无法获取最小Eigenvalue问题排查

• 程序可正常返回最大特征值，但无法获取最小特征值
• 测试2^3×23矩阵时无此问题，官方文档示例运行正常
• 代码运行耗时不到3秒，编译约10秒

完整代码 #

#include <iostream>
#include <random>
#include <Eigen/Sparse>
#include <Spectra/GenEigsSolver.h>
#include <Spectra/MatOp/SparseGenMatProd.h>
using namespace std;
using uint = unsigned int;

#define ALLOCATE_ARRAY(type, count) ((type*) malloc(sizeof(type) * (count)))
#define REALLOC_ARRAY(array, type, count) ((type*) realloc(array, sizeof(type) * count))
#define MATRIX(array, i,j) array[j+i*total_sites]


uint number_sites = 4;
uint total_sites = 0;
uint total_Hilbert_elements = 0;
uint number_spin_states = 2;
int* computational_base = 0;

default_random_engine ran_generator;

void inizialize_Hamiltonian_matrix (float trasverse_magnetic_field, float longitudinal_magnetic_field, uint boundary_conditions_flag, Eigen::MatrixXf& Hamiltonian_config_matrix);

int main(){

    float beta = 0.3;
    float trasverse_magnetic_field = 1;
    float longitudinal_magnetic_field = 1;

    uint dimensions = 1;
    uint initial_state = 1;
    uint boundary_conditions_flag = 1;
    uint number_eigenvalues = 1;
    total_sites = pow(number_sites, dimensions);

    total_Hilbert_elements = pow(number_spin_states, total_sites);
    computational_base =  ALLOCATE_ARRAY(int, total_Hilbert_elements*total_sites);

    for (int i=0; i<total_Hilbert_elements; ++i){
        uint i_holder = i;
        for (int j=0; j<total_sites; ++j){
            MATRIX(computational_base, i,j) = i_holder%2;
            i_holder = i_holder/2;
        }
    }

    // Build Matrix and print on shell
    Eigen::MatrixXf Hamiltonian_config_matrix(total_Hilbert_elements, total_Hilbert_elements);
    inizialize_Hamiltonian_matrix (trasverse_magnetic_field, longitudinal_magnetic_field, boundary_conditions_flag, Hamiltonian_config_matrix);

    for (int i=0; i<total_Hilbert_elements; ++i){
        for (int j=0; j<total_Hilbert_elements; ++j){
            if (Hamiltonian_config_matrix(i,j)>=0){
                cout << " " << Hamiltonian_config_matrix(i,j);
            }
            else{
                cout << Hamiltonian_config_matrix(i,j);
            }
        }
        cout << endl;
    }

    // Converting in Sparse Matrix and Spectra stuff
    Eigen::SparseMatrix<float> Sparse_Hamiltonian_config;
    Sparse_Hamiltonian_config = Hamiltonian_config_matrix.sparseView();
    Spectra::SparseGenMatProd<float> Operation_object_wrap(Sparse_Hamiltonian_config);

    // Initialize and Compute
    Spectra::GenEigsSolver<Spectra::SparseGenMatProd<float>> Ready_for_Lanczos(Operation_object_wrap, number_eigenvalues, 2*number_eigenvalues+1);
    Ready_for_Lanczos.init();

    uint which_eigenvalues = 0;
    cout << "Type 1 if Largest Eigenvalue, Type 0 if Smallest" << endl;
    cin >> which_eigenvalues;

    if (which_eigenvalues==0){
        int nconv = Ready_for_Lanczos.compute(Spectra::SortRule::SmallestMagn);
    }
    else{
        int nconv = Ready_for_Lanczos.compute(Spectra::SortRule::LargestMagn);
    }

    // Retrieve results
    Eigen::VectorXcf evalues;
    if(Ready_for_Lanczos.info() == Spectra::CompInfo::Successful)
        evalues = Ready_for_Lanczos.eigenvalues();
    cout << "Eigenvalues found:\n" << evalues << endl;

    return 0;
}

void inizialize_Hamiltonian_matrix (float trasverse_magnetic_field, float longitudinal_magnetic_field, uint boundary_conditions_flag, Eigen::MatrixXf& Hamiltonian_config_matrix){
    for (uint j=0; j<total_sites-1; ++j){
        for (uint i=0; i<total_Hilbert_elements; ++i){
            if (MATRIX(computational_base, i,j)==MATRIX(computational_base, i,j+1)){
                Hamiltonian_config_matrix(i,i) = Hamiltonian_config_matrix(i,i) -1.0;
            }
            else{
                Hamiltonian_config_matrix(i,i) = Hamiltonian_config_matrix(i,i) +1.0;
            }
        }
    }
    if (boundary_conditions_flag==1){
        uint j = total_sites-1;
        for (uint i=0; i<total_Hilbert_elements; ++i){
            if (MATRIX(computational_base, i,j)==MATRIX(computational_base, i,0)){
                    Hamiltonian_config_matrix(i,i) = Hamiltonian_config_matrix(i,i) -1.0;
            }
            else{
                Hamiltonian_config_matrix(i,i) = Hamiltonian_config_matrix(i,i) +1.0;
            }
        }
    }
    for (uint j=0; j<total_sites-1; ++j){
        for (uint i=0; i<total_Hilbert_elements; ++i){
            if (MATRIX(computational_base, i,j)==1){
                Hamiltonian_config_matrix(i,i) = Hamiltonian_config_matrix(i,i) -longitudinal_magnetic_field;
            }
            else{
                Hamiltonian_config_matrix(i,i) = Hamiltonian_config_matrix(i,i) +longitudinal_magnetic_field;
            }
        }
    }
    uint new_i_binary = 0;
    for (uint j=0; j<total_sites-1; ++j){
        for (uint i=0; i<total_Hilbert_elements; ++i){
            if (MATRIX(computational_base, i,j)==1){
                new_i_binary = i - pow(2, j);
            }
            else{
                new_i_binary = i + pow(2, j);
            }
            Hamiltonian_config_matrix(new_i_binary,i) = Hamiltonian_config_matrix(new_i_binary,i) -trasverse_magnetic_field;
        }
    }
}

问题修复建议 #

1. 初始化矩阵时清零
   Eigen的MatrixXf默认不会将元素初始化为0，导致矩阵存在随机垃圾值，干扰特征值计算。在创建矩阵后添加清零操作：

   Eigen::MatrixXf Hamiltonian_config_matrix(total_Hilbert_elements, total_Hilbert_elements);
   Hamiltonian_config_matrix.setZero(); // 新增此行
   inizialize_Hamiltonian_matrix(trasverse_magnetic_field, longitudinal_magnetic_field, boundary_conditions_flag, Hamiltonian_config_matrix);
   

2. 增大Lanczos子空间维度
   寻找最小模特征值时，当前设置的子空间维度2*number_eigenvalues+1（即3）过小，不足以让算法收敛。建议增大到更大的值，比如10：

   Spectra::GenEigsSolver<Spectra::SparseGenMatProd<float>> Ready_for_Lanczos(Operation_object_wrap, number_eigenvalues, 10);
   

3. 替换pow为位运算避免精度问题
   pow(2,j)返回浮点数，赋值给uint可能因精度丢失导致索引错误，改用位运算翻转比特位：

   // 替换原new_i_binary计算部分
   new_i_binary = i ^ (1 << j);
   

内容的提问来源于stack exchange，提问作者Matteo