Bug 230 - ARPACK support
ARPACK support
Status: NEW
Product: Eigen
Classification: Unclassified
Component: Sparse
All All
: --- enhancement
Assigned To: Nobody
Depends on:
  Show dependency treegraph
Reported: 2011-03-22 00:16 UTC by Gael Guennebaud
Modified: 2014-07-09 14:17 UTC (History)
3 users (show)

Patch for initial version of ArpackSelfAdjointEigenSolver (20.01 KB, application/octet-stream)
2012-07-24 06:28 UTC, David Harmon
no flags Details
Test driver for ArpackSelfAdjointEigenSolver, -larpack needed to build. (2.24 KB, application/octet-stream)
2012-07-24 06:29 UTC, David Harmon
no flags Details
ARPACK support (31.30 KB, application/octet-stream)
2012-10-07 01:21 UTC, David Harmon
no flags Details
New patch fixing some tiny bugs including EIGEN_STATIC_ASSERTs (3.68 KB, patch)
2013-06-18 08:06 UTC, Kolja Brix
no flags Details | Diff

Description Gael Guennebaud 2011-03-22 00:16:47 UTC
One way to get eigen decomposition on sparse matrices is to wrap ARPACK++, here is a start:

and here is a copy of schizofrenikh's message from the forum:

    #include <arpack++/arssym.h>
    #include <vector>
    #include <eigen2/Eigen/Eigen>

    class TMatrixForArpackpp {
        Eigen::SparseMatrix< double > &M;
        TMatrixForArpackpp( Eigen::SparseMatrix< double > &_M ) :
            M(_M) {
        void MultMv( double *v, double *w ) {
            Eigen::Map< Eigen::VectorXd > v_map = Eigen::VectorXd::Map( v, M.rows() );
            Eigen::Map< Eigen::VectorXd > w_map = Eigen::VectorXd::Map( w, M.cols() );
            w_map = M.marked<Eigen::SelfAdjoint|Eigen::LowerTriangular>()*v_map;

    void Lanczos_symm( Eigen::SparseMatrix< double > &M, std::vector< double > &eigenValues, std::vector< Eigen::VectorXd > &eigenVectors ) {
        int n = M.cols();
        TMatrixForArpackpp M_arp( M );
        ARSymStdEig< double, TMatrixForArpackpp > prob( n, n-1, &M_arp, &TMatrixForArpackpp::MultMv, "SM" ); //"SM" to get smallest eigenvalues
        int m = prob.ConvergedEigenvalues();
        eigenValues.assign( prob.RawEigenvalues(), prob.RawEigenvalues()+m );
        for( int i=0; i<m; i++ )
            eigenVectors[i] = Eigen::VectorXd::Map( prob.RawEigenvector(i), n);
        /*eigenValues = *( prob.StlEigenvalues() );
        eigenVectors = *( prob.StlEigenvectors() );*/

To compile with arpack on my debian system I had to add  -lblas -lgfortran -lgfortranbegin -lnsl -larpack -larpack++
Comment 1 David Harmon 2012-07-21 19:42:13 UTC
ARPACK++ is difficult to compile on many platforms, and even if you get there, does not support 64-bit mode, which can be a deal-breaker for a lot of people using other external libraries.

Anyways, I have a simple wrapper for ARPACK that I propose submitting after I clean it up and add some more functionality, since ARPACK packages are readily available on pretty much all *nix platforms (including Mac OS X).

I copied the dense eigen decomposition interface for my wrapper, but I have a few questions on the "Eigen way" of doing things before I submit:

1) Storage type for eigenvectors. In the dense eigen-code storage is MatrixType, which is fine since it is solving for all eigenvectors. However, in this case, if we only want m << n vectors, then the storage should be n x m.

I currently use Matrix<Rows, Dynamic, Scalar>. Is this ok? (Rows is a typedef of the nbr of rows, and can be Dynamic as well)

2) If an unsupported Scalar type is used, is it better to have a runtime assertion tripped or compile-time error? The advantage I see of a runtime assertion is that I can explicitly say "That type is not supported", but a compile-time error will be more cryptic when it doesn't find the necessary template implementation.

3) How do I call the Eigen SimplicialLLT code from a header file? I need a linear solver for generalized eigen problems. I tried to use SimplicialLLT, but calling from a header file gave a bunch of errors. Are there known ways around this?

4) Is there an interface for cleanly allowing the user to choose between linear solver types that I can copy? I would like to give the option of built-in solver, CHOLMOD, UMFPACK, etc. I'm sure there's other code that does this, and I can just copy what it does if someone points me to it.

5) Do we trust what the user gives us? For example, in the generalized problem Av=\lambda B v, different choices must be made based on if B is PD, PSD, etc. Do we trust what the user says or do some error checking?

The choices I see are a) always trust what the user says, and if we get garbage, then so be it. b) never trust the user, check the rank of B. c) check the rank of B if NDEBUG is not defined.

That's all my questions for now.
Comment 2 Gael Guennebaud 2012-07-21 22:22:27 UTC
1) It should simply be Matrix<Scalar,Dynamic,Dynamic>, the size of the input matrix is assumed to be large, so Dynamic.

2) Compile-time error. To this end you can use a static assertion in the general not supported template instantiation:

If NUMERIC_TYPE_MUST_BE_REAL is not appropriate, we can add one in Core/util/StaticAssert.h

3) I don't really understand your issue. You must include Eigen/SparseCholesky, and then you can use SimplicialLDLT. I guess that ARPACK requires a function pointer. ou have to write it as a wrapper around SimplicialLDLT. Ideally, this wrapper function could be a template function supporting any solver, and the solver class should be specified by the user with a default to SimplicialLDLT.

4) Half answered in 3). You could make it a template parameter of ArpackEigenSolver:

template<typename MatrixType, typename Solver = SimplicialLDLT<MatrixType> >
class ArpackEigenSolver {

5) We trust what the user say, but we need to report to the user when ARPACK failed. Our standard is to report errors via the .info() method which returns Success or NumericalIssue, NoConvergence, InvalidInput.
Comment 3 David Harmon 2012-07-22 19:02:50 UTC
Thanks for the info.

2) Small issue, but just calling


  gives me a "NUMERIC_TYPE_MUST_BE_REAL is not member of static_assertion<false>" error. I explicitly put in static_assertion<true>::NUMERIC_TYPE_MUST_BE_REAL and the error went away. I didn't seen anywhere else in Eigen doing this, so I must be doing something wrong.

3) Dumb mistake on my part. I wasn't following the Eigen convention of letting the caller be responsible for including necessary header files. if I put include Eigen/SparseCholesky in my ArpackSelfAdjointEigenSolver header file, crazy things happen, but if my caller simply says

#include <Eigen/SparseCholesky>
#include <ArpackSelfAdjointEigenSolver>

  and uses the arpack solver, then everything is fine.

4) No wrapper is needed. Arpack computes stuff, returns it to you, you perform the linear solve, then call again, repeat until convergence.

I seem to recall the old linear solver module using a template to choose the solver (Cholmod, Umfpack, etc.). I'll follow that format.

I'll have the symmetric and generalized symmetric code ready in the next week or two. It's a lot of repeated code, but since the dense solver does it, I can create two separate classes for SelfAdjoint problems and GeneralizedSelfAdjoint problems.
Comment 4 Gael Guennebaud 2012-07-23 20:58:28 UTC
"lot of repeated code"-> there must be something wrong. Perhaps you could share your code even if it is not ready at all. Better to get early feedbacks to avoid useless work!
Comment 5 David Harmon 2012-07-24 06:25:55 UTC
I'm attaching a patch of the class ArpackSelfAdjointEigenSolver that solves for the specified number of eigenvalues/eigenvectors for both dense and sparse matrices, with either float or double Scalar type. Error checking is included as well. I've put the code in unsupported/, mirroring the dense setup, except naming the module "ArpackSupport" instead of "Eigenvalues".

I'm also attaching a test driver that generates random dense/sparse symmetric matrices and computes the eigenvalues/vectors, and then prints out || A v - \lambda v || for each one.
Comment 6 David Harmon 2012-07-24 06:28:47 UTC
Created attachment 288 [details]
Patch for initial version of ArpackSelfAdjointEigenSolver
Comment 7 David Harmon 2012-07-24 06:29:43 UTC
Created attachment 289 [details]
Test driver for ArpackSelfAdjointEigenSolver, -larpack needed to build.
Comment 8 David Harmon 2012-10-07 01:19:58 UTC
I've had a chance to work on this some more and have added support for the generalized eigenvalue problem. I've tested it for a variety of problems and have been using it in my own research. I'd like it to go "live" as an unsupported module so others can start using/testing it.

Currently it only works with self-adjoint problems, and does not support the arpack functionality of finding all eigenvalues around a given value (still trying to get it to work).

I created a new file in unsupported/Eigen/ArpackSupport. The source is in unsupported/Eigen/src/Eigenvalues/. The code is pretty well commented (I think).

The patch with the new files will follow, or should I commit to the repository directly?
Comment 9 David Harmon 2012-10-07 01:21:37 UTC
Created attachment 298 [details]
ARPACK support
Comment 10 David Harmon 2012-12-04 02:08:41 UTC
I've had some people contact me after applying the patch manually themselves and using this extension. Could somebody with write access to the repository please just check it in so people can start using it properly?

If there is some other way I should submit the patch, please let me know.

Comment 11 Kolja Brix 2013-06-18 08:06:41 UTC
Created attachment 350 [details]
New patch fixing some tiny bugs including EIGEN_STATIC_ASSERTs
Comment 12 Kolja Brix 2013-06-18 08:09:02 UTC
Unfortunately, there were still some tiny bugs in the ARPACK wrapper, i.e. wrong a include file name. These issues are fixed with the new patch including the EIGEN_STATIC_ASSERT issue.

Thanks to David and Gael for proving the ARPACK wrapper, special thanks to Gael for helping with fixing the EIGEN_STATIC_ASSERT issue.
Comment 13 Gael Guennebaud 2013-06-18 09:45:28 UTC
Changeset:   c5057358673a
User:        Kolja Brix
Date:        2013-06-18 09:44:40
Summary:     Bug 230, fix compilation issues and wrong static assertions
Comment 14 Jim Garrison 2013-08-16 02:11:12 UTC
The code in ArpackSelfAdjointEigenSolver.h says "by default, the upper triangular part is used, but can be changed through the template parameter."  However, the documentation for the regular (dense) SelfAdjointEigenSolver says "Only the lower triangular part of the input matrix is referenced."  To keep consistency, should ArpackSelfAdjointEigenSolver be changed to instead look at the lower triangular part of the matrix by default?
Comment 15 Gael Guennebaud 2013-08-20 13:26:23 UTC
I agree with you that the default should be to read the lower part, but that's unclear how to do it, and the default linear solver is SimplicialLLT<Scalar> which already reads the lower part. So, it might be that the doc is wrong.
Comment 16 Jim Garrison 2014-07-05 08:51:31 UTC
Also, I have found that ArpackGeneralizedSelfAdjointEigenSolver fails to compile if used with a dense matrix type -- it only seems to work with sparse matrices.  It would be nice if it worked with any type of matrix.

Actually, even more important: it would be nice if we did not have to store the entire sparse matrix in memory at all if it can be easily computed on the fly.  For instance, I have a function that takes a row index, then calls a callback for each nonzero column entry in that row.  I would love to be able to use this directly with Eigen::VectorXd.
Comment 17 Gael Guennebaud 2014-07-09 14:17:54 UTC
You can already write a small class inheriting from SparseMatrixBase and implementing an InnerIterator to loop through the non-zeros. A good example is SparseView which iterate over the non-zeros of a dense matrix.

Note You need to log in before you can comment on or make changes to this bug.