Generalized Eigenvalues

[yixuan]

Moved from yixuan/arpack-eigen#2

[skn123]

Hi Yixuan, any update on this?

[yixuan]

Ah, sorry for the slow progress on this. When I find some free time I'll catch up.

[Barry845]

Nice work.
Right now I'm using my own C++ interface to ARPACK-NG based on Arpack++ and in-house code. But your code looks quite light and easy to integrate (we are using Eigen for our project), the only thing missing is Generalized Eigenvalue support.
Any plans on doing this ?

[yixuan]

One thing I was struggling with is how to design the API of generalized solver: there are many different modes to compute the generalized eigenvalues, according to ARPACK manual. How to unify these interfaces is a problem. Maybe I should learn some design from ARPACK++, or do you have any suggestions? @skn123 @Barry845

[skn123]

I think a simple API would suffice. LHS matrix and a RHS matrix

[manusethi]

Hi @yixuan
In one of the comments on arpack-eigen you mentioned that we can construct our own operations class to handle the generalized eigenvalue problem. Could you give some direction on how that can be done?

[yixuan]

Hi @manusethi , please see the comments below.

[manusethi]

Thanks @yixuan. However, I do not want to invert the matrix. Both my matrices are sparse and inverting one of them will make it dense.

[yixuan]

Computing generalized eigenvalue does require some form of matrix inversion, either on the A matrix or on the B matrix. For sparse matrix there is a sparse Cholesky decomposition algorithm, which in Eigen is done by the SimplicialLLT solver.

Actually I have already written a generalized eigen solver using Cholesky decomposition, but the documentation is not ready yet.

[manusethi]

Thanks @yixuan, I appreciate it indeed! Could you help me with calling the eigs function for this class or provide a brief documentation? I am having some problem providing the cholesky decomposition.

[yixuan]

How do you store your sparse matrix? You may want to convert it to an Eigen sparse matrix, and then I could show you some sample code.

[manusethi]

Hi @yixuan I have already converted my matrix into a Eigen sparse matrix. I just don't know how to call your function. I get some errors. I tried a few times. Here is what I did:

Eigen::SparseMatrix<double> A  ;
Eigen::SparseMatrix<double> B ;

/*
* I have not specified above how I have allocated 
* values in A and B. But Both A and B are correct 
* and I have verified their values.
* After that I did the following:
*/

Spectra::SparseGenMatProd<double> op1(A) ;
Spectra::SparseGenMatProd<double> op2(B) ;

/* 
* After that I don't know how to proceed.
* I also tried the cholesky but don't know how to specify in your new function.
*/

Eigen::SimplicialLLT< Eigen::SparseMatrix<double> > lltofB(B) ;

Spectra::SymGEigsSolver<double, Spectra::SMALLEST_MAGN,
        Spectra::SparseGenMatProd<double>, Eigen::SimplicialLLT< Eigen::SparseMatrix<double>>, Spectra::GEIGS_CHOLESKY > eigs(&op1, &lltofB, nvec, nvec+6);

// I also tried the following:

Spectra::SymGEigsSolver<double, Spectra::SMALLEST_MAGN,
        Spectra::SparseGenMatProd<double>, Spectra::SparseGenMatProd<double>, Spectra::GEIGS_CHOLESKY > eigs(&op1, &op2, nvec, nvec+6);

I didn't really understand how to call so just tried a few things and for sure these wouldn't work. I would greatly appreciate a sample code.

Also, does Arpack also take the inverse of the B matrix while computing generalized eigenvalues?

[yixuan]

Hi all, just to let you know that finally I've added a generalized eigen solver to Spectra. The documentation and example code can be found here.

This implemented solver is the most basic one, and other more complicated algorithms are under development, but I think it can already solve a large number of symmetric generalized eigenvalue problems.

@manusethi, ARPACK also requires you to implement the operation inv(B) * v, and you need to program that by yourself. Now in Spectra this can be done via the wrapper class SparseCholesky. See the code in the documentation.

[manusethi]

Hi @yixuan. Thank you very much.

BTW regarding the inverse, I have been using ARPACK++. As per their documentation I didn't need to program the inverse of B matrix. But probably they are taking the inverse like what you have described.

[manusethi]

Hi @yixuan

I just tried your code and it worked great! Now I no longer have to deal with Arpack++ and other dependencies like SuperLU, Umfpack, BLAS etc which were creating a lot of problems for me when using with my MPI based code. Thank you so much!

[yixuan]

Glad to hear that!

[xchen-cs]

Hey @yixuan
Any plan to support the case when B is indefinite?

[yixuan]

@TerranRush
If I get it right, you are talking about the case where B is indefinite but is invertible, correct?

[xchen-cs]

Yes, just as the ARPACK buckling mode.

[yixuan]

@TerranRush
I plan to implement the "regular inverse mode" first, i.e., solving inv(B) * A * x = lambda * x. Does this suffice to your problem?

[xchen-cs]

@yixuan
Actually, I don't know the algorithm details. According to ARPACK document, I'm looking for this one:
c Mode 4: K_x = lambda_KG_x, K symmetric positive semi-definite,
c KG symmetric indefinite
c ===> OP = (inv[K - sigma_KG])*K and B = K.
c ===> Buckling mode
I tried the unsupported ARPACK module in Eigen, but it looks this case is not well handled.
So I found your library here :)

[yixuan]

It depends on whether you want to calculate the largest (or smallest) generalized eigenvalues, or the ones that are closest (or farthest) to some given number sigma. The first case will use regular inverse mode, and the latter corresponds to the buckling mode.

[oberbichler]

Hi!
Is it possible to combine SymGEigsShiftSolver with a ShiftSolve operator? I try to calculate the smallest generalized eigenvalues with large symetric sparse matrices (~60000x60000).

[KainFanfon]

I also need to calculate the smallest generalized eigenvalues with large
symetric sparse matrices. Did you (@oberbichler) find a way? Did you solve your problem?
There was a contribution attempt by Anna Araslanova (LOBPCG algorithm), is this
already available?

[oberbichler]

Finally we used another algorithm. Here you will find our implementation.

[yixuan]

@oberbichler Thanks for the pointer. I'm also trying to learn new algorithms that better fit such problems.

[lobpcg]

Finally we used another algorithm. Here you will find our implementation.

Anna Araslanova (LOBPCG algorithm) is already available in contributed and should likely be faster. You might want to give it a try.

[oberbichler]

Thanks for the hint! A lot of compiling errors... still WIP?

[yixuan]

Yeah, sorry about that...
I just experimentally merged the PR and didn't have time to write a test. Probably I can finish that in one or two weeks. The release of v0.7 of Spectra will be a signal.

[oberbichler]

👍 Would be nice to check the perfomance with large systems.

After fixing the compile time errors it drops a runtime error (invalid matrix product). I will wait for the next release 😛

[yixuan]

For people who are still interested, SymGEigsShiftSolver has been implemented in the 1.y.z branch. Example code is available in the testing program.

Note that this new branch contains several API-breaking changes that will be included in the next major release, so some (minor) code migration is needed.

[yixuan]

I think it is ready to close this issue now, as Spectra v1.0.0 has been released and generalized eigenvalues are fully supported.

It is also funny to see how long I took to close the very first issue of Spectra. 😆