Repository files navigation

NOTE: Spectra 1.0.0 was released in 2021-07-01, with a lot ofAPI-breaking changes. Please see themigration guidefor a smooth transition to the new version.

NOTE: If you are interested in the future development of Spectra, please jointhis thread to share your comments and suggestions.

Spectra stands forSparseEigenvalueComputationToolkitas aRedesignedARPACK. It is a C++ library for large scale eigenvalueproblems, built on top ofEigen,an open source linear algebra library.

Spectra is implemented as a header-only C++ library, whose only dependency,Eigen, is also header-only. HenceSpectra can be easily embedded inC++ projects that require calculating eigenvalues of large matrices.

ARPACK is a software package written inFORTRAN for solving large scale eigenvalue problems. The development ofSpectra is much inspired by ARPACK, and as the full name indicates,Spectra is a redesign of the ARPACK library using the C++ language.

In fact,Spectra is based on the algorithm described in theARPACK Users' Guide,the implicitly restarted Arnoldi/Lanczos method. However,it does not use the ARPACK code, and it isNOT a clone of ARPACK for C++.In short,Spectra implements the major algorithms in ARPACK,butSpectra provides a completely different interface, and it does notdepend on ARPACK.

Spectra is designed to calculate a specified number (k) of eigenvaluesof a large square matrix (A). Usuallyk is much smaller than the size of the matrix(n), so that only a few eigenvalues and eigenvectors are computed, whichin general is more efficient than calculating the whole spectral decomposition.Users can choose eigenvalue selection rules to pick the eigenvalues of interest,such as the largestk eigenvalues, or eigenvalues with largest real parts, etc.

To use the eigen solvers in this library, the user does not need to directlyprovide the whole matrix, but instead, the algorithm only requires certain operationsdefined onA. In the basic setting, it is simply the matrix-vectormultiplication. Therefore, if the matrix-vector productA * x can be computedefficiently, which is the case whenA is sparse,Spectrawill be very powerful for large scale eigenvalue problems.

There are two major steps to use theSpectra library:

1. Define a class that implements a certain matrix operation, for example thematrix-vector multiplicationy = A * x or the shift-solve operationy = inv(A - σ * I) * x.Spectra has defined a number ofhelper classes to quickly create such operations from a matrix object.See the documentation ofDenseGenMatProd,DenseSymShiftSolve, etc.
2. Create an object of one of the eigen solver classes, for exampleSymEigsSolverfor symmetric matrices, andGenEigsSolverfor general matrices. Member functionsof this object can then be called to conduct the computation and retrieve theeigenvalues and/or eigenvectors.

Below is a list of the available eigen solvers inSpectra:

• SymEigsSolver:For real symmetric matrices
• GenEigsSolver:For general real- and complex-valued matrices
• SymEigsShiftSolver:For real symmetric matrices using the shift-and-invert mode
• GenEigsRealShiftSolver:For general real matrices using the shift-and-invert mode,with a real-valued shift
• GenEigsComplexShiftSolver:For general real matrices using the shift-and-invert mode,with a complex-valued shift
• SymGEigsSolver:For generalized eigen solver with real symmetric matrices
• SymGEigsShiftSolver:For generalized eigen solver with real symmetric matrices, using the shift-and-invert mode
• DavidsonSymEigsSolver:Jacobi-Davidson eigen solver for real symmetric matrices, with the DPR correction method
• HermEigsSolver:For complex Hermitian matrices

Below is an example that demonstrates the use of the eigen solver for symmetricmatrices.

#include<Eigen/Core>#include<Spectra/SymEigsSolver.h>// <Spectra/MatOp/DenseSymMatProd.h> is implicitly included#include<iostream>usingnamespaceSpectra;intmain(){// 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;return0;}

Sparse matrix is supported via classes such asSparseGenMatProdandSparseSymMatProd.

#include<Eigen/Core>#include<Eigen/SparseCore>#include<Spectra/GenEigsSolver.h>#include<Spectra/MatOp/SparseGenMatProd.h>#include<iostream>usingnamespaceSpectra;intmain(){// A band matrix with 1 on the main diagonal, 2 on the below-main subdiagonal,// and 3 on the above-main subdiagonalconstint n =10;    Eigen::SparseMatrix<double>M(n, n);    M.reserve(Eigen::VectorXi::Constant(n,3));for(int i =0; i < n; i++)    {        M.insert(i, i) =1.0;if(i >0)            M.insert(i -1, i) =3.0;if(i < n -1)            M.insert(i +1, i) =2.0;    }// Construct matrix operation object using the wrapper class SparseGenMatProd    SparseGenMatProd<double>op(M);// Construct eigen solver object, requesting the largest three eigenvalues    GenEigsSolver<SparseGenMatProd<double>>eigs(op,3,6);// Initialize and compute    eigs.init();int nconv = eigs.compute(SortRule::LargestMagn);// Retrieve results    Eigen::VectorXcd evalues;if(eigs.info() == CompInfo::Successful)        evalues = eigs.eigenvalues();    std::cout <<"Eigenvalues found:\n" << evalues << std::endl;return0;}

And here is an example for user-supplied matrix operation class.

#include<Eigen/Core>#include<Spectra/SymEigsSolver.h>#include<iostream>usingnamespaceSpectra;// M = diag(1, 2, ..., 10)classMyDiagonalTen{public:using Scalar =double;// A typedef named "Scalar" is requiredintrows()const {return10; }intcols()const {return10; }// y_out = M * x_invoidperform_op(constdouble *x_in,double *y_out)const    {for(int i =0; i <rows(); i++)        {            y_out[i] = x_in[i] * (i +1);        }    }};intmain(){    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;    }return0;}

When it is needed to find eigenvalues that are closest to a numberσ,for example to find the smallest eigenvalues of a positive definite matrix(in which caseσ = 0), it is advised to use the shift-and-invert modeof eigen solvers.

In the shift-and-invert mode, selection rules are applied to1/(λ - σ)rather thanλ, whereλ are eigenvalues ofA.To use this mode, users need to define the shift-solve matrix operation. Seethe documentation ofSymEigsShiftSolverfor details.

Spectra provides theHermEigsSolver solver for complex-valued Hermitian matrices,and theGenEigsSolver solver for general complex-valued matrices. See the example below.

#include<Eigen/Core>#include<Spectra/HermEigsSolver.h>#include<Spectra/GenEigsSolver.h>#include<iostream>usingnamespaceSpectra;intmain(){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 DenseGenMatProdusing OpHType = DenseHermMatProd<std::complex<double>>;using OpGType = DenseGenMatProd<std::complex<double>>;    OpHTypeopH(H);    OpGTypeopG(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-valuedif (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);if (eigsG.info() == CompInfo::Successful)    {        Eigen::VectorXcd 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;    }return0;}

TheAPI reference page contains the documentationofSpectra generated byDoxygen,including all the background knowledge, example code and class APIs.

More information can be found in the project pagehttps://spectralib.org.

An optional CMake installation is supported, if you have CMake with at least v3.10 installed. You can install the headers using the following commands:

mkdir build&&cd buildcmake .. -DCMAKE_INSTALL_PREFIX='intended installation directory' -DBUILD_TESTS=TRUEmake all&& maketest&& make install

By installingSpectra in this way, you also create a CMake targetSpectra::Spectra that can be used in subsequent build procedures for other programs.

If you already have Eigen installed, you can specify the installation directory by setting theCMAKE_PREFIX_PATH variable orEigen3_ROOT.For example:

cmake .. -DCMAKE_INSTALL_PREFIX='intended installation directory' -DCMAKE_PREFIX_PATH='path where the installation of Eigen3 can be found' -DBUILD_TESTS=TRUE

A couple of useful environment variables can be set to control the download of Eigen.CPM_DOWNLOAD_ALL=ON will force the download of Eigen, even if an installation is already present on the system.CPM_LOCAL_PACKAGES_ONLY=ON will force the opposite behavior. The download directory can be controlled by setting the variableCPM_SOURCE_CACHE.

Spectra is an open source project licensed underMPL2, the same license used byEigen.