ARPACK-Armadillo
SymEigsShiftSolver< Scalar, SelectionRule, OpType > Class Template Reference

#include <SymEigsSolver.h>

Inheritance diagram for SymEigsShiftSolver< Scalar, SelectionRule, OpType >:
SymEigsSolver< Scalar, SelectionRule, OpType >

Public Member Functions

 SymEigsShiftSolver (OpType *op_, int nev_, int ncv_, Scalar sigma_)
 
- Public Member Functions inherited from SymEigsSolver< Scalar, SelectionRule, OpType >
 SymEigsSolver (OpType *op_, int nev_, int ncv_)
 
void init (Scalar *init_resid)
 
void init ()
 
int compute (int maxit=1000, Scalar tol=1e-10)
 
int num_iterations ()
 
int num_operations ()
 
Vector eigenvalues ()
 
Matrix eigenvectors (int nvec)
 
Matrix eigenvectors ()
 

Detailed Description

template<typename Scalar = double, int SelectionRule = LARGEST_MAGN, typename OpType = DenseSymShiftSolve<double>>
class SymEigsShiftSolver< Scalar, SelectionRule, OpType >

This class implements the eigen solver for real symmetric matrices using the shift-and-invert mode. The background information of the symmetric eigen solver is documented in the SymEigsSolver class. Here we focus on explaining the shift-and-invert mode.

The shift-and-invert mode is based on the following fact: If \(\lambda\) and \(x\) are a pair of eigenvalue and eigenvector of matrix \(A\), such that \(Ax=\lambda x\), then for any \(\sigma\), we have

\[(A-\sigma I)^{-1}x=\nu x\]

where

\[\nu=\frac{1}{\lambda-\sigma}\]

which indicates that \((\nu, x)\) is an eigenpair of the matrix \((A-\sigma I)^{-1}\).

Therefore, if we pass the matrix operation \((A-\sigma I)^{-1}y\) (rather than \(Ay\)) to the eigen solver, then we would get the desired values of \(\nu\), and \(\lambda\) can also be easily obtained by noting that \(\lambda=\sigma+\nu^{-1}\).

The reason why we need this type of manipulation is that the algorithm of ARPACK-Armadillo (and also ARPACK) is good at finding eigenvalues with large magnitude, but may fail in looking for eigenvalues that are close to zero. However, if we really need them, we can set \(\sigma=0\), find the largest eigenvalues of \(A^{-1}\), and then transform back to \(\lambda\), since in this case largest values of \(\nu\) implies smallest values of \(\lambda\).

To summarize, in the shift-and-invert mode, the selection rule will apply to \(\nu=1/(\lambda-\sigma)\) rather than \(\lambda\). So a selection rule of LARGEST_MAGN combined with shift \(\sigma\) will find eigenvalues of \(A\) that are closest to \(\sigma\). But note that the eigenvalues() method will always return the eigenvalues in the original problem (i.e., returning \(\lambda\) rather than \(\nu\)), and eigenvectors are the same for both the original problem and the shifted-and-inverted problem.

Template Parameters
ScalarThe element type of the matrix. Currently supported types are float and double.
SelectionRuleAn enumeration value indicating the selection rule of the shifted-and-inverted eigenvalues. The full list of enumeration values can be found in SelectionRule.h .
OpTypeThe name of the matrix operation class. Users could either use the DenseSymShiftSolve wrapper class, or define their own that impelemnts all the public member functions as in DenseSymShiftSolve.

Below is an example that illustrates the use of the shift-and-invert mode:

#include <armadillo>
#include <SymEigsSolver.h> // Also includes <MatOp/DenseSymShiftSolve.h>
int main()
{
// A size-10 diagonal matrix with elements 1, 2, ..., 10
arma::mat M(10, 10, arma::fill::zeros);
for(int i = 0; i < M.n_rows; i++)
M(i, i) = i + 1;
// Construct matrix operation object using the wrapper class
// Construct eigen solver object with shift 0
// This will find eigenvalues that are closest to 0
DenseSymShiftSolve<double> > eigs(&op, 3, 6, 0.0);
eigs.init();
eigs.compute();
arma::vec evalues = eigs.eigenvalues();
evalues.print("Eigenvalues found:"); // Will get (3.0, 2.0, 1.0)
return 0;
}

Also an example for user-supplied matrix shift-solve operation class:

#include <armadillo>
#include <SymEigsSolver.h>
// M = diag(1, 2, ..., 10)
class MyDiagonalTenShiftSolve
{
private:
double sigma_;
public:
int rows() { return 10; }
int cols() { return 10; }
void set_shift(double sigma) { sigma_ = sigma; }
// y_out = inv(A - sigma * I) * x_in
// inv(A - sigma * I) = diag(1/(1-sigma), 1/(2-sigma), ...)
void perform_op(double *x_in, double *y_out)
{
for(int i = 0; i < rows(); i++)
{
y_out[i] = x_in[i] / (i + 1 - sigma_);
}
}
};
int main()
{
MyDiagonalTenShiftSolve op;
// Find three eigenvalues that are closest to 3.14
MyDiagonalTenShiftSolve> eigs(&op, 3, 6, 3.14);
eigs.init();
eigs.compute();
arma::vec evalues = eigs.eigenvalues();
evalues.print("Eigenvalues found:"); // Will get (4.0, 3.0, 2.0)
return 0;
}

Definition at line 426 of file SymEigsSolver.h.

Constructor & Destructor Documentation

template<typename Scalar = double, int SelectionRule = LARGEST_MAGN, typename OpType = DenseSymShiftSolve<double>>
SymEigsShiftSolver< Scalar, SelectionRule, OpType >::SymEigsShiftSolver ( OpType *  op_,
int  nev_,
int  ncv_,
Scalar  sigma_ 
)
inline

Constructor to create a eigen solver object using the shift-and-invert mode.

Parameters
op_Pointer to the matrix operation object, which should implement the shift-solve operation of \(A\): calculating \((A-\sigma I)^{-1}y\) for any vector \(y\). Users could either create the object from the DenseSymShiftSolve wrapper class, or define their own that impelemnts all the public member functions as in DenseSymShiftSolve.
nev_Number of eigenvalues requested. This should satisfy \(1\le nev \le n-1\), where \(n\) is the size of matrix.
ncv_Parameter that controls the convergence speed of the algorithm. Typically a larger ncv_ means faster convergence, but it may also result in greater memory use and more matrix operations in each iteration. This parameter must satisfy \(nev < ncv \le n\), and is advised to take \(ncv \ge 2\cdot nev\).
sigma_The value of the shift.

Definition at line 459 of file SymEigsSolver.h.


The documentation for this class was generated from the following file: