Spectra: Sparse Eigenvalue Computation Toolkit

Spectra calculates a specified number ($k$) of eigenvalues of a large square matrix ($A$). Instead of requiring the full matrix, the algorithm only requires certain operations (like matrix-vector multiplication $Ax$).

To use the library, follow these two steps:

1. Define a matrix operation class: This class must implement the required operation (e.g., $y=Ax$ or the shift-solve $y=(A-\sigma I)^{-1}x$). You can use Spectra's built-in helper classes (like DenseSymMatProd or SparseGenMatProd) or implement your own custom class.
2. Create and run a solver object: Instantiate one of the available solver classes (e.g., SymEigsSolver for symmetric matrices) by passing the operation object, the number of eigenvalues requested, and the number of Lanczos vectors. Call .init() and .compute() to perform the calculation, then use .eigenvalues() or .eigenvectors() to retrieve results.

If you have a custom matrix or a specific way to compute $y=Ax$, you can implement your own operation class. Your class must provide:

• A using Scalar = ...; typedef.
• int rows() const and int cols() const methods.
• A void perform_op(const Scalar *x_in, Scalar *y_out) const method that performs the multiplication.

#include <Eigen/Core>
#include <Spectra/SymEigsSolver.h>
#include <iostream>

using namespace Spectra;

// M = diag(1, 2, ..., 10)
class MyDiagonalTen
{
public:
    using Scalar = double;  // A typedef named "Scalar" is required
    int rows() const { return 10; }
    int cols() const { return 10; }
    // y_out = M * x_in
    void perform_op(const double *x_in, double *y_out) const
    {
        for(int i = 0; i < rows(); i++)
        {
            y_out[i] = x_in[i] * (i + 1);
        }
    }
};

int main()
{
    MyDiagonalTen op;
    SymEigsSolver<MyDiagonalTen> eigs(op, 3, 6);
    eigs.init();
    eigs.compute(SortRule::LargestAlge);
    if(eigs.info() == CompInfo::Successful)
    {
        Eigen::VectorXd evalues = eigs.eigenvalues();
        std::cout << "Eigenvalues found:\n" << evalues << std::endl;
    }

    return 0;
}

This example uses DenseSymMatProd to wrap a dense Eigen matrix and SymEigsSolver to find the largest three eigenvalues.

#include <Eigen/Core>
#include <Spectra/SymEigsSolver.h>
// <Spectra/MatOp/DenseSymMatProd.h> is implicitly included
#include <iostream>

using namespace Spectra;

int main()
{
    // We are going to calculate the eigenvalues of M
    Eigen::MatrixXd A = Eigen::MatrixXd::Random(10, 10);
    Eigen::MatrixXd M = A + A.transpose();

    // Construct matrix operation object using the wrapper class DenseSymMatProd
    DenseSymMatProd<double> op(M);

    // Construct eigen solver object, requesting the largest three eigenvalues
    SymEigsSolver<DenseSymMatProd<double>> eigs(op, 3, 6);

    // Initialize and compute
    eigs.init();
    int nconv = eigs.compute(SortRule::LargestAlge);

    // Retrieve results
    Eigen::VectorXd evalues;
    if(eigs.info() == CompInfo::Successful)
        evalues = eigs.eigenvalues();

    std::cout << "Eigenvalues found:\n" << evalues << std::endl;

    return 0;
}

Spectra supports complex-valued matrices using HermEigsSolver for Hermitian matrices and GenEigsSolver for general complex matrices. Note that for Hermitian matrices, eigenvalues are real-valued while eigenvectors are complex-valued.

#include <Eigen/Core>
#include <Spectra/HermEigsSolver.h>
#include <Spectra/GenEigsSolver.h>
#include <iostream>

using namespace Spectra;

int main()
{
    std::srand(0);

    // We are going to calculate the eigenvalues of H and G
    Eigen::MatrixXcd G = Eigen::MatrixXcd::Random(10, 10);
    // H is Hermitian
    Eigen::MatrixXcd H = G + G.adjoint();

    // Construct matrix operation objects using the wrapper
    // classes DenseHermMatProd and DenseGenMatProd
    using OpHType = DenseHermMatProd<std::complex<double>>;
    using OpGType = DenseGenMatProd<std::complex<double>>;
    OpHType opH(H);
    OpGType opG(G);

    // Construct solver object for H, requesting the largest three eigenvalues
    HermEigsSolver<OpHType> eigsH(opH, 3, 6);

    // Initialize and compute
    eigsH.init();
    int nconvH = eigsH.compute(SortRule::LargestAlge);

    // Retrieve results
    // Eigenvalues are real-valued, and eigenvectors are complex-valued
    if (eigsH.info() == CompInfo::Successful)
    {
        Eigen::VectorXd evaluesH = eigsH.eigenvalues();
        std::cout << "Eigenvalues of H found:\n" << evaluesH << std::endl;
        Eigen::MatrixXcd evecsH = eigsH.eigenvectors();
        std::cout << "Eigenvectors of H:\n" << evecsH << std::endl;
    }

    // Similar procedure for matrix G
    GenEigsSolver<OpGType> eigsG(opG, 3, 6);
    eigsG.init();
    int nconvG = eigsG.compute(SortRule::LargestMagn);

    // Retrieve results
    Eigen::VectorXcd evaluesG;
    if (eigsG.info() == CompInfo::Successful)
    {
        evaluesG = eigsG.eigenvalues();
        std::cout << "Eigenvalues of G found:\n" << evaluesG << std::endl;
        Eigen::MatrixXcd evecsG = eigsG.eigenvectors();
        std::cout << "Eigenvectors of G:\n" << evecsG << std::endl;
    }

    return 0;
}