"TridiagEigen: eigen decomposition failed" when switching from 0.9.0 to 1.x

[unclejimbo]

Hi, I'm computing the smallest eigenvalues for a Laplacian matrix.

In 0.9.0 version the following code works fine:

using Operator = Spectra::SparseSymShiftSolve<double>;
using Solver = Spectra::SymEigsShiftSolver<double, Spectra::LARGEST_MAGN, Operator>;
Eigen::SparseMatrix<double> L; // the Laplacian matrix
Operator op(L);
Solver eigensolver(&op, k, 2*k+1, 0.0); // for some k
eigensolver.init();
eigensolver.compute();

But when I switch to 1.x (master branch), it throws an exception

TridiagEigen: eigen decomposition failed

with the following code:

using Operator = Spectra::SparseSymShiftSolve<double>;
using Solver = Spectra::SymEigsShiftSolver<Operator>;
Eigen::SparseMatrix<double> L; // the Laplacian matrix
Operator op(L);
Solver eigensolver(op, k, 2*k+1, 0.0); // for some k
eigensolver.init();
eigensolver.compute(Spectra::SortRule::LargestMagn);

Am I getting it wrong? Any help is appreciated, thanks.

[yixuan]

Hi @unclejimbo, do you have any example matrix L that causes this error?

[unclejimbo]

@yixuan Hi, the matrix I'm testing is the cotangent matrix for the Laplace-Beltrami operator on triangular meshes, which should be symmetric. Below is the test code, you can comment out the igl part I used to generate the cotagent matrix, and read in the file I uploaded (see bottom). Check out the matrix generation code in details here if you have time.

Intestringly enough, computing 5 eigenvalues as done in the code is ok. But when I try to compute >=20 eigenvalues it raise the exception.

#include <Eigen/Core>
#include <Eigen/Sparse>
#include <Eigen/SparseExtra>
#include <igl/cotmatrix.h>
#include <igl/readOFF.h>
#include <Spectra/MatOp/SparseSymMatProd.h>
#include <Spectra/MatOp/SparseSymShiftSolve.h>
#include <Spectra/MatOp/SymShiftInvert.h>
#include <Spectra/SymEigsShiftSolver.h>
#include <Spectra/SymGEigsShiftSolver.h>

int main(int argc, char* argv[])
{
    Eigen::MatrixXd V;
    Eigen::MatrixXi F;
    igl::readOFF("mesh.off", V, F);

    Eigen::SparseMatrix<double> L;
    igl::cotmatrix(V, F, L);
    L = -L;

    Eigen::saveMarket(L, "cotmat");

    using Operator = Spectra::SparseSymShiftSolve<double>;
    using Solver = Spectra::SymEigsShiftSolver<Operator>;
    Operator op(L);
    Solver eigensolver(op, 5, 15, 0.0f);
    eigensolver.init();
    eigensolver.compute(Spectra::SortRule::LargestMagn);

    return 0;
}

cotmat.txt

[yixuan]

Hi @unclejimbo, sorry for the delay. The issue is actually easy to fix. The point is that the matrix L has an eigenvalue of zero, so it cannot be used as a shift. If you already know that the matrix is nonnegative positive and you want to find the smallest eigenvalues, then you can use a shift slightly smaller than zero. See the example below.

#include <iostream>
#include <fstream>
#include <sstream>
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <Eigen/Eigenvalues>
#include <Spectra/MatOp/SparseSymShiftSolve.h>
#include <Spectra/SymEigsShiftSolver.h>

int main(int argc, char* argv[])
{
    std::ifstream file("cotmat.txt");
    std::string line;

    std::getline(file, line);
    std::cout << "read 1st line: " << line << std::endl;
    std::getline(file, line);
    std::cout << "read 2nd line: " << line << std::endl;
    std::istringstream iss(line);
    int i, j;
    iss >> i >> j;
    Eigen::SparseMatrix<double> L(i, j);

    double val;
    while(std::getline(file, line))
    {
        std::istringstream iss(line);
        iss >> i >> j >> val;
        L.insert(i - 1, j - 1) = val;
    }

    std::cout << L.rows() << " rows, " << L.cols() << " columns, " << L.nonZeros() << " nnz" << std::endl;

    // Dense eigen solver to verify
    Eigen::MatrixXd mat = L;
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> sol(mat);
    std::cout << sol.eigenvalues().head(20) << std::endl << std::endl;

    using Operator = Spectra::SparseSymShiftSolve<double>;
    using Solver = Spectra::SymEigsShiftSolver<Operator>;
    Operator op(L);
    Solver eigensolver(op, 20, 41, -0.1);
    eigensolver.init();
    eigensolver.compute(Spectra::SortRule::LargestMagn);
    if(eigensolver.info() == Spectra::CompInfo::Successful)
        std::cout << eigensolver.eigenvalues() << std::endl;

    return 0;
}

read 1st line: %%MatrixMarket matrix coordinate  real general
read 2nd line: 1250 1250 8738
1250 rows, 1250 columns, 8738 nnz
6.41554e-15
  0.0189065
  0.0189065
  0.0189065
  0.0383656
  0.0383656
  0.0711537
  0.0711537
  0.0711537
  0.0766289
  0.0766289
  0.0766289
  0.0854201
   0.118646
   0.164199
   0.164199
   0.164199
   0.179575
   0.179575
   0.222099

    0.222099
    0.179575
    0.179575
    0.164199
    0.164199
    0.164199
    0.118646
   0.0854201
   0.0766289
   0.0766289
   0.0766289
   0.0711537
   0.0711537
   0.0711537
   0.0383656
   0.0383656
   0.0189065
   0.0189065
   0.0189065
-6.38378e-16