Search code examples
c++eigenlapackarmadilloeigen3

Does the SymEigsShiftSolver of Spectra not return eigenvectors?


I have implemented the SymEigsShiftSolver for computing the eigenvalues of a large sparse matrix, however it does not return me the eigenvectors. Has it not been implemented as of yet?

void Eigens::computeEigenvectors(Matrices m)
{
SparseSymShiftSolve<double> op(m.Lpl);
SymEigsShiftSolver< double, SMALLEST_MAGN, SparseSymShiftSolve<double> >
            eigs(&op, k, 4, 0.0);
eigs.init();
int nconv = eigs.compute();
VectorXd evalues;
evalues.resize(k);
if(eigs.info() == SUCCESSFUL)
    evalues = eigs.eigenvalues();
cout << "Eigenvalues found:\n" << evalues << endl;
cout <<"\nHere is the matrix whose columns are eigenvectors of the Laplacian Matrix \n"
     <<"corresponding to these eigenvalues: \n"
     <<eigs.eigenvectors()<<endl;
}

Solution

  • I'm not sure why you've added the "armadillo" tag, since you're using the eigen library. I've provided the following response if you are indeed interested in an Armadillo based solution.

    Both Armadillo and Spectra use very similar underlying code for sparse eigendecomposition (the code was written by the same author), but Armadillo has a simplified user interface. To compute the eigenvectors of a symmetric sparse matrix using Armadillo, use the eigs_sym() function:

    // generate sparse symmetric matrix
    sp_mat A = sprandu<sp_mat>(5000, 5000, 0.1);
    sp_mat B = A.t()*A;
    
    vec eigval;
    mat eigvec;
    
    eigs_sym(eigval, eigvec, B, 5);       // find 5 eigenvectors with largest magnitude
    eigs_sym(eigval, eigvec, B, 5, "sm"); // find 5 eigenvectors with smallest magnitude
    

    If you have a non-symmetric matrix, the eigs_gen() function can be used in a similar manner:

    sp_mat A = sprandu<sp_mat>(5000, 5000, 0.1);  
    
    cx_vec eigval;
    cx_mat eigvec;
    
    eigs_gen(eigval, eigvec, A, 5);  // find 5 eigenvalues/eigenvectors
    

    A sparse matrix can be constructed from a 1D array in a simple manner:

    double data[nrows * ncols]; // 1D array representation of your matrix
    
    sp_mat X = sp_mat( mat(data, nrows, ncols, false) );