Armadillo eigs_sym: decomposition failed

Question

I am writing C++ code within Visual Studio and using Armadillo 7.900.1

I am having no luck getting eigs_sym function working under Armadillo (I am using the lapack and blas versions that came with Armadillo). The typical error message is as below:

error: lapack::stedc(): failed to compute all eigenvalues
warning: eigs_sym(): decomposition failed

The code that is producing this output is:

sp_mat T(locations, values);
arma::eigs_sym(eigval, eigvec, T, num_eigs_wanted, "sm", tol);

where

T is a 480,000x480,000 sparse matrix. 

I can get the code working if T is small (ie. 2000x2000) and tol is high (tol = 5). But once T is about 20000x20000, then no matter what values of tol or num_eigs_wanted are used then the above error occurs. (I obviously change "values" and locations" to alter the size of T).

The matrix T is symmetric, real, positive definite.

The exact same matrix has no problems when I call the eigs function within Matlab.

Any ideas? To me it seems that eigs_sym in Armadillo is broken...any alternatives that people have used?

Cheers

Answer 1

Not an answer but way too long for a comment. I can reproduce with the code below. The matrix is just a diagonal matrix with ascending values on the diagonal starting at 1. The eigenvalues are obviously 1,2,3,... and the matrix is positive definite and symmetric.

#include <iostream>
#include <armadillo>

int main()
{
  constexpr int N = 20000;
  arma::umat locations(2,N);
  arma::vec values(N);

  for (int i = 0; i < N; ++i)
  {
    locations(0,i) = i;
    locations(1,i) = i;
    values(i) = i+1;
  }

  arma::sp_mat T(locations, values);

  int num_eigs_wanted = 1;
  arma::vec eigval(num_eigs_wanted);
  arma::mat eigvec(N,num_eigs_wanted);
  double tol = 1e-6;
  arma::eigs_sym(eigval, eigvec, T, num_eigs_wanted, "sm", tol);

  std::cout << eigval << '\n';
}

Compiler flags are

clang++-5.0 -std=c++11 test.cpp -larmadillo

I'm using Armadillo with MKL backend and I receive:

warning: eigs_sym(): decomposition failed
[matrix size: 0x1]

---

It seems to be a problem of Armadillo, because when I use the unsupported Arpack support of Eigen it works just fine.

#include <iostream>
#include <Eigen/Sparse>
#include <unsupported/Eigen/ArpackSupport>

int main()
{
  constexpr int N = 20000;

  std::vector< Eigen::Triplet<double> > triplets;
  triplets.reserve(N);
  for (int i = 0; i < N; ++i)
    triplets.push_back( {i,i,i+1} );

  Eigen::SparseMatrix<double> T(N,N);
  T.setFromTriplets(triplets.begin(), triplets.end());

  int num_eigs_wanted = 1;
  double tol = 1e-6;
  Eigen::ArpackGeneralizedSelfAdjointEigenSolver< Eigen::SparseMatrix<double> > eigs_sym;
  eigs_sym.compute(T, num_eigs_wanted, "SM", Eigen::ComputeEigenvectors, tol);
  std::cout << eigs_sym.eigenvalues() << '\n';
}

You might also want to take a look at https://github.com/yixuan/spectra/