Bugzilla – Bug 230

ARPACK support

Last modified: 2014-08-27 10:21:07 UTC

One way to get eigen decomposition on sparse matrices is to wrap ARPACK++, here is a start: http://forum.kde.org/viewtopic.php?f=74&t=84238&p=192268#p192198 and here is a copy of schizofrenikh's message from the forum: #include <arpack++/arssym.h> #include <vector> #include <eigen2/Eigen/Eigen> class TMatrixForArpackpp { private: Eigen::SparseMatrix< double > &M; public: 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 prob.FindEigenvectors(); //prob.FindEigenvalues(); int m = prob.ConvergedEigenvalues(); eigenValues.resize(m); eigenVectors.resize(m); 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++

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.

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: STATIC_ASSERT(false,NUMERIC_TYPE_MUST_BE_REAL); 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.

Thanks for the info. 2) Small issue, but just calling EIGEN_STATIC_ASSERT(false, NUMERIC_TYPE_MUST_BE_REAL); 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.

"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!

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.

Created attachment 288 [details] Patch for initial version of ArpackSelfAdjointEigenSolver

Created attachment 289 [details] Test driver for ArpackSelfAdjointEigenSolver, -larpack needed to build.

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?

Created attachment 298 [details] ARPACK support

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. Thanks.

Created attachment 350 [details] New patch fixing some tiny bugs including EIGEN_STATIC_ASSERTs

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.

https://bitbucket.org/eigen/eigen/commits/c5057358673a/ Changeset: c5057358673a User: Kolja Brix Date: 2013-06-18 09:44:40 Summary: Bug 230, fix compilation issues and wrong static assertions

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?

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.

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.

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.

I've been working on some improvements to Eigen's ARPACK support at <https://bitbucket.org/garrison/eigen-arpack>. I made some minor cleanups to ArpackSelfAdjointEigenSolver.hpp, and I also have created ArpackEigenSolver.hpp. Currently ArpackEigenSolver works only for the type std::complex<double>. It should be easy to also extend to work with std::complex<float>. Working with regular float and double will be a bit tricky and involve a bit of refactoring, since ARPACK takes a slightly different set of arguments when called for real types. There are a few suboptimal things with the current version, which can be noticed by doing a diff against ArpackSelfAdjointEigenSolver.hpp. I tried to keep the differences to a minimum, since much of the code has been copy-pasted from there and tweaked. Also, the patch series contains a patch similar to the one provided for bug #858, a patch which I made before noticing that report.

I had a quick look to your changesets, and I don't really understand why the sigma variable has to be of type Scalar while it is really real. Regarding the new class, would it make sense to introduce a common base class to reduce duplicated code?

Gael, you are right about sigma. I tried to make it consistent in both versions but failed. In ArpackSelfAdjointEigenSolver sigma should be type RealScalar, but in ArpackEigenSolver it should be type ComplexScalar. I am very unfamiliar with hg so I am unsure how to modify my revision history to fix this. (Actually, I remember hearing once that hg discourages such modifications, so I'm unsure what to do instead.) For background, there are three classes problems as far as ARPACK is concerned (see http://www.caam.rice.edu/software/ARPACK/UG/node32.html): (i) real symmetric, (ii) real non-symmetric, and (iii) complex. Each is supported for both single and double precision. Notice that for complex types, the user is supposed to call the same routine for both Hermitian and non-symmetric problems. There's really two possible solutions as far as Eigen's API is concerned: make them aware of this difference, and tell them to use ArpackEigenSolver any time the matrix is complex (even if it is Hermitian); or, have ArpackSelfAdjointEigensolver abstract over this difference by accepting complex matrices, taking a triangular view, and then in the end converting the eigenvalues to real before returning. If the latter option is chosen, we would also need to abstract over the "which" field with an enum, because e.g. if one wants the smallest algebraic eigenvalue using ArpackSelfAdjointEigensolver, the "which" string would need to be "SA" for real problems and "SR" for complex problems.

FWIW, scipy chooses the latter option I mentioned: if eigsh (sparse hermitian) is called with a complex (Hermitian) matrix, it calls znaupd and simply discards the (small) imaginary part of each returned eigenvalue.

yes I'd also go the scipy way. Regarding code factorization, I was thinking about something like our CholmodSupport module (https://bitbucket.org/eigen/eigen/src/a1f5fed0d734c2993474868950af76a6889ba9f2/Eigen/src/CholmodSupport/CholmodSupport.h?at=default) where a CholmodBase class implements all the common code to interface cholmod while the derived classes (CholmodSimplicialLLT, CholmodSimplicialLDLT, CholmodDecomposition) simply provide the end-user specific interface. That being said, I don't known if there is much to factorize or not as I did not look carefully to your wrapper code.