Re: [eigen] Band Lanczos Algorith

Ok, so I have finished a first cut of the class and a matching SVD solver that is faster for large matrices than JacobiSVD. I also have added 8 tests that demonstrate them. There are a couple of issues.

1) Right now it doesn't inherent from any of the solvers, I might need you guys' help figuring out how different stopping conditions should match up.

2) It seems like someone should have already created a complex conjugate function that is similar to C++11's (ie works for double type), but I couldn't find it so I added my own.

3) I templated over Scalar rather than Matrix, it seemed like it would be a lot of work to template over Matrix type so I thought I should just get a decent working version first. I looked at the SelfAdjointEigenSolver for example but it looked ... confusing. If templating over MatrixType is preferring it would be helpful if someone could give me tips.

Please take a look at the fork I made, its here: https://bitbucket.org/micahc/eigen. (see unsupported/Eigen/src/IterativeSolvers/BandLanczosSelfAdjointEigenSolver.h and unsupported/Eigen/src/IterativeSolvers/TruncatedLanczosSVD.h)

Thanks!

-Micah

On Wed, Dec 24, 2014 at 7:35 AM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:

On Tue, Dec 23, 2014 at 9:55 PM, Micah Chambers <micahc.vt@xxxxxxxxx> wrote:
Hey everyone. As part of my PhD work I needed an implementation of the Band Lanczos Algorithm, so I implemented it from the "Templates for the Solution of Algebraic Eigenvalue Problems" book (link) using eigen, but its not really a full blown template class. I figured it might be worth expanding and adding to the Unsupported set of libraries. I have a couple of questions about contributing, since it will still be a bit of work to make as generic as the rest of the decomposition classes.

1) Does this already exist somewhere?

AFAIK, nothing like this yet.

2) How does Eigen deal with Hermitian/non-Hermitian naming, since usually a very similar algorithm exists for both. Should I make a separate class or just add a flag for whether the matrix is hermitian?

So far we used to have different classes because the underlying algorithms where really different. Implementation-wise, if the code is the same, you can factorize it through a common internal free function, and if the API is also common, then a common base class will do it. We use the word "SelfAdjoint" instead of Hermitian though.

3) In these sorts of problems its sometimes possible that the Matrix being decomposed can't actually be formulated, so a function needs to passed that takes a vector and returns a vector. Are there any classes that already do this so I can make sure my API looks the same?

What does this function do? A matrix vector product? In general, we request that the abstract matrix class follows the part of the Eigen API required by the algorithm, like implementing rows(), cols(), and operator*(const Eigen::MatrixBase<Derived> &) for iterative solvers as BiCGSTAB.

4) Which branch should I check out to add stuff to?

The default one.

5) Can I just put the tests in unsupported/test?

That's the right place ! You might look at the test/bicgstab.cpp and test/conjugate_gradient.cpp files for some examples which makes use of a common set of testing functions for sparse solvers (in test/sparse_solver.h)

Feel free to ask more questions as the internal mechanisms of Eigen are not yet documented.

cheers,
Gael

Thanks

-Micah

--
Micah Chambers, M.S.
Ahmanson-Lovelace Brain Mapping Center
University of California Los Angeles
635 Charles E. Young Drive SW, Suite 225
Los Angeles, CA 90095-7334
310-206-2101 (Office)
310-206-5518 (Fax)
703-307-7692 (Cell)
email: micahcc@xxxxxxxx

Micah Chambers, M.S.

Ahmanson-Lovelace Brain Mapping Center

University of California Los Angeles

635 Charles E. Young Drive SW, Suite 225

Los Angeles, CA 90095-7334

310-206-2101 (Office)
310-206-5518 (Fax)

703-307-7692 (Cell)

email: micahcc@xxxxxxxx