Bugzilla – Bug 230

ARPACK support

Last modified: 2014-07-09 14:17:54 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.