Search code examples
sparse-matrixarmadilloeigenvalue

Armadillo: eigs_gen for smallest eigenvalue

I'm using armadillo's eigs_gen to find the smallest algebraic eigenvalue of a sparse matrix.

If I request the function for just the smallest eigenvalue the result is incorrect but if I request it for the 2 smallest eigenvalues the result is correct. The code is:

#include <iostream>
#include <armadillo>

using namespace std;
using namespace arma;


int
main(int argc, char** argv)
  {
  cout << "Armadillo version: " << arma_version::as_string() << endl;

  sp_mat A(5,5);

  A(1,2) = -1;
  A(2,1) = -1;
  A(3,4) = -1;
  A(4,3) = -1;

  cx_vec eigval;
  cx_mat eigvec;

  eigs_gen(eigval, eigvec, A, 1, "sr");  // find smallest eigenvalue ---> INCORRECT RESULTS
  eigval.print("Smallest real eigval:");

  eigs_gen(eigval, eigvec, A, 2, "sr");  // find 2 smallest eigenvalues ---> ALMOST CORRECT RESULTS 
  eigval.print("Two smallest real eigvals:");

  return 0;
  }

My compile command is: 

 g++ file.cpp -o file.exe -O2 -I/path-to-armadillo/armadillo-4.600.3/include -DARMA_DONT_USE_WRAPPER -lblas -llapack -larpack

The output is:

Armadillo version: 4.600.3 (Off The Reservation)
Smallest real eigval:
    (+1.000e+00,+0.000e+00)
Two smallest real eigvals:
    (-1.000e+00,+0.000e+00)
    (-1.164e-17,+0.000e+00)

Any idea on why this is happening and how to overcome this is appreciated.

Note: second result is only almost correct because we expect -1, -1 as the two lowest eigenvalues but perhaps repeated eigenvalues are ignored.

Update: including a test matrix construction which, after ryan's changes to include the "sa" option to the library, doesn't seem to converge:

    #define ARMA_64BIT_WORD
    #include <armadillo>
    #include <iostream>
    #include <vector>
    #include <stdio.h>

    using namespace arma;
    using namespace std;

    int main(){

    size_t l(3), ls(l*l*l);
    sp_mat A = sprandn<sp_mat>(ls, ls, 0.01);
    sp_mat B = A.t()*A;

    vec eigval;
    mat eigvec;
    eigs_sym(eigval, eigvec, B, 1, "sa");

    return 0;

    }

The matrix sizes of interest are much larger e.g. ls = 8000 - 27000, and is not quite the matrix constructed here but I presume the problem should be the same.

Solution

  • I believe that the issue here is that you are running eigs_gen() (which calls DNAUPD) on a symmetric matrix. ARPACK notes that DNAUPD is not meant for symmetric matrices, but does not specify what will happen if you use symmetric matrices anyway:

    NOTE: If the linear operator "OP" is real and symmetric with respect to the real positive semi-definite symmetric matrix B, i.e. B*OP = (OP')*B, then subroutine ssaupd should be used instead.

    (from http://www.mathkeisan.com/usersguide/man/dnaupd.html )

    I modified the internal Armadillo code to pass "sa" (smallest algebraic) to the ARPACK calls in eigs_sym() (sp_auxlib_meat.hpp), and I was able to obtain the correct eigenvalues. I've submitted a patch upstream to make "sa" and "la" support available for eigs_sym(), which I think should solve your problem once a new version is released (or at some point in the future).