Speedup for spectra for matrix sizes of around 1M!

[ghost]

Hi Yixuan,

I would like to appreciate your help for taking care of the issue I had with the EigenSymShiftSolver. I have tried to leverage these techniques to try and calculate the eigenvalues of a sparse 40,000 x 40,000 matrix. The spectra library fails to generate the eigenvalues for some reason. The same source code manages to return the eigenvalues for a 30,000 x 30,000 matrix. Is it a limitation of the library?

       int main()
       {
           SparseMatrix<double> Lpl(40000,40000);
           Lpl.reserve(Eigen::VectorXi::Constant(40000, 3));
           for(int i = 0; i < 10; i++)
               {
                   Lpl.insert(i, i) = 1.0;
                 if(i > 0)
                      Lpl.insert(i - 1, i) = 3.0;
                 if(i < 10 - 1)
                        Lpl.insert(i + 1, i) = 2.0;
                  }

             SparseGenMatProd<double> op(Lpl);
           GenEigsSolver< double, LARGEST_MAGN, SparseGenMatProd<double> > eigs(&op, 3, 6);
             eigs.init();
            int nconv = eigs.compute();
            VectorXcd evalues;
            evalues.resize(3);
             if(eigs.info() == SUCCESSFUL)
                   evalues = eigs.eigenvalues();
             cout << "Eigenvalues found:\n" << evalues << endl;

             cout <<"\nHere is the matrix whose columns are eigenvectors of the Laplacian Matrix \n"
            <<"corresponding to these eigenvalues: \n"
               << eigs.eigenvectors() <<endl;

              return 0;
          }

[ghost]

I believe it was a compiler issue. It works perfectly with matrix as big as 80k.

[yixuan]

Great. And also note that the matrix in your example is not symmetric.

[ghost]

Hi Yixuan, may I know which eigen solving algorithm is being leveraged by your library? Also is there anyway I can speed up the whole computation process for sparse matrices as big as 80k or even 1M x 1M? I understand reducing the number of maximum iterations might be a solution, but do you have any suggestion other than that? The issue is that I am computing the Laplacian matrix of a graph and aiming to compute its eigenvectors and for few cases the size is as big as 50k to 80k, in which cases it takes around 20-30 seconds to bring out the first 3 eigenvectors. I was wondering it would be great if this computational time can be sped up.

[yixuan]

It is based on the Implicitly Restarted Arnoldi Method, first implemented by the ARPACK package. I think it is not straightforward to gain general improvement on the algorithm, but if your matrix has special structures, there may be some way to do better. Do you have any information on this? For example, the matrix is banded, or is block-diagonal, etc.

[ghost]

Hi Yixuan,
I was researching a lot on eigensolvers. I also tried to use PETSc/SLEPc but I am failing to obtain major speedup there since I do not wish to run my code on multiple processors.

From what I have observed from the code in spectra/include/MatOp/SparseSymShiftSolve.h, you are using the SimplicialLDLT to calculate (A-\sigma I)^{-1}x. I also noticed from the Eigen library that the SimplificialLDLT is recommended for very sparse and not too large problems (e.g., 2D Poisson eq.). Perhaps SparseLU or even ConjugateGradient might work for the case that I am talking about.(Solve the first 3 or 4 eigenvectors of a sparse symmetric matrix of size around 500k x 500k or even 1 million, that is not block-diagonal or banded). It could lead to significant speedup for large scale problems. I am struggling with trying to incorporate this idea in your library. Would love to hear your suggestion and obtain your guidance with this issue.

Thanks a lot!

[yixuan]

It is possible to use different solvers to implement the inv(A - sigma * I) * x operation. Basically you just modify the SparseSymShiftSolve class to use other solvers such as LU or CG.
I think you can first do some basic benchmarking comparing the speed of sparse solvers. Do you have any success on this?

[ghost]

Yes, I have a list of ISPD benchmark matrices which I have used. The time for computation is around 5-6 times more than that of PETSc/SLEPc for those matrices. I believe it is the linear solver and the preconditioner used in the present algorithm which is slowing things down. A same version of spectra but with the CG and ILU preconditioner implemented would result in considerable speedup, is my firm belief.

[yixuan]

One more question before I try to write a solver for you: do you only need the smallest eigenvalues such that sigma can be assumed to be zero, or do you need a general selection rule that is based on the value of sigma?

[ghost]

Thanks for the wonderful response.
I am looking to extract out the first 3~4 smallest eigen values. So I guess a general selection rule works better.

[yixuan]

Oh, I should have explained my question more clearly. If you only want the smallest k eigenvalues then it corresponds to the case of sigma = 0, and if you need eigenvalues closest to some given number s, then it means sigma = s.
In your case I can assume sigma = 0, which makes the implementation easier.

[yixuan]

You can try the example attached below. Note that in this case you have to make sure that the matrix is positive definite.

cg-example.zip

[ghost]

Thanks a lot!. This will surely help!

[yixuan]

@bodhi1991 Do you have any success on this?

[GYengera]

I am working with a sparse symmetric positive semidefinite matrix. Is the Sparse Symmetric Shift solver the most optimum solver for computing the minimum eigenvalue of such a matrix?

[ghost]

Hi Yixuan,
Apologies for the delay. I had been travelling. Working with the conjugate gradient and incomplete LUT did not give me any improvement at all, which is very strange. I have been trying to figure out what could be the reason behind this but did not get major breakthrough here. Maybe the Arnoldi is not suitable for large scale eigen problems.

[yixuan]

@bodhi1991 I think you can first test the CG solver in one iteration, and see if is faster than an LDLT solver. And then you may compare the number of iterations via the num_iterations() and num_operations() functions for different solvers. If you have a reproducible demo example to share, I can also help you to debug.

[yixuan]

@GYengera It's hard to say "optimum", but generally SymEigsShiftSolver is a good starting point. You can always benchmark different software packages to pick the best for you.

[y-yao]

Hi Yixuan,

Why are the dimension of the matrices indicated by a signed integer in the code? For example the m_n variable in SymEigsBase.h. This limits the size of the matrices that can be diagonalized using Spectra. For sparse matrices, it is not rare to go beyond 2 billion in dimension. Wouldn't size_t (unsigned long long) be a better type here?

Thank you.

[yixuan]

@y-yao I may change that to Eigen::Index, which is expanded to std::ptrdiff_t by default. Will let you know when ready. Thanks.