SymEigsSolver misses some eigenvalues

[pmacg]

I am adapting the first example from the Quick Start page to compute eigenvalues of a graph Laplacian matrix. I'm using the code given at the bottom of the issue, which takes arguments n and k where n is the number of vertices in your graph (dimensions of the Laplacian matrix) and k is the number of eigenvalues to compute.

The code:

• constructs the normalised Laplacian of the cycle graph on n vertices
• uses the example code from the Quick Start page to compute k eigenvalues

When I call the program with different values of k, I get inconsistent results.

>> ./demo 20 5
Eigenvalues found:
      2
1.95106
1.80902
1.58779
1.30902

>> ./demo 20 6
Eigenvalues found:
      2
1.95106
1.95106
1.80902
1.80902
1.58779

Notice that the multiplicities of the eigenvalues seem to be ignored when asking for 5 eigenvalues, but the multiple eigenvalues are returned correctly when we ask for 6.

Smallest Magnitude

If I modify the code to look for the smallest magnitude eigenvalues, then the results are also inconsistent.

>> ./demo 20 5
Eigenvalues found:
        1
 0.690983
 0.412215
 0.190983
0.0489435

>> ./demo 20 6
Eigenvalues found:
    0.690983
    0.412215
    0.190983
   0.0489435
   0.0489435
-9.62807e-23

When asking for 5 eigenvalues, the smallest eigenvalue (which is 0) is not found, and multiplicities are ignored. When asking for 6 eigenvalues, the smallest eigenvalue is found although the multiplicities are still not correct - every eigenvalue other than 0 should have multiplicity 2.

Code

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

// Function declarations
Eigen::MatrixXd construct_cycle_laplacian(int n);

int main(int argc, char *argv[])
{
  // Assume that we have been given two arguments: n and k
  //   n is the dimension of the matrix to construct
  //   k is the number of eigenvalues to compute
  int n = atoi(argv[1]);
  int k = atoi(argv[2]);

  // We are going to calculate the eigenvalues of M
  Eigen::MatrixXd M = construct_cycle_laplacian(n);

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

  // Construct eigen solver object, requesting the largest k eigenvalues
  Spectra::SymEigsSolver<Spectra::DenseSymMatProd<double>> eigs(op, k, fmin(2*k, n));

  // Initialize and compute
  eigs.init();
  eigs.compute(Spectra::SortRule::LargestMagn);

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

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

  return 0;
}

Eigen::MatrixXd construct_cycle_laplacian(int n) {
  // Initialise the laplacian matrix
  Eigen::MatrixXd L(n, n);

  // Add the matrix entries, iterating over the rows
  for (int i = 0; i < n; i++) {
    L(i, i) = 1;
    L(i, (i+n-1) % n) = -0.5;
    L(i, (i+1) % n) = -0.5;
  }

  return L;
}

[yixuan]

Hi @pmacg, thanks for providing these examples. Computing eigenvalues with multiplicities is a known difficult task. I tested the example above with ARPACK, and also occasionally lost the repeated eigenvalues. This is a limitation of the algorithm itself.

Typically you can get more reliable results by increasing ncv and using a smaller tol. For example, with n=20, k=5, you can get the desired results using eigs(op, k, 15) and eigs.compute(Spectra::SortRule::LargestMagn, 1000, 1e-16).

For the smallest eigenvalues, the suggested way is using the shift-and-invert mode. See the SymEigsShiftSolver for details.

[pmacg]

Thanks, I appreciate your response. I've been able to get everything working nicely with the tips you provided, and taking a look at your added tests in e333245.

I'll close this issue.

[yixuan]

Thanks for the feedback!