Re: [eigen] Band Lanczos Algorith

On Tue, Jan 13, 2015 at 1:54 PM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:

Hi,

sorry for the late reply. I haven't found time to look at the fork yet, but let me try reply some of your questions anyways.

On Mon, Dec 29, 2014 at 7:01 AM, Micah Chambers <micahc.vt@xxxxxxxxx> wrote:
Ok, so I have finished a first cut of the class and a matching SVD solver that is faster for large matrices than JacobiSVD.

For large matrices you should rather compare to the BDCSVD class implementing a much faster divide and conquer algorithm.

I also have added 8 tests that demonstrate them.

in test/svd_common.h we have a set of common routines to verify the reliability of SVD solvers. Also look at test/*svd.cpp files to see how to use 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.

If it is working on a BandMatrix or SparseMatrix, then you won't find a common base class to inherit from.

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.

Yes we have one as numext::conj(Scalar), and many others as numext::real(Scalar), numext::imag(Scalar), etc.

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.

Actually, templating on a matrix type is just a convenient way to specify many attributes at once: scalar type, compile time sizes, maximal compile-time sizes, expected input storage (dense/sparse/etc.). Then the implementer is free to use whatever is best suited for the internal matrices and vectors. As a start you can only consider the provided Scalar type and neglect the other information.

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)

will do....

thank you for sharing your efforts!

gael

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

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