yixuan / arpack-eigen

A header-only C++ redesign of ARPACK using the Eigen library
23 stars 4 forks source link

Sparse Matrices #1

Closed skn123 closed 9 years ago

skn123 commented 9 years ago

Hi Yixuan Do you have any examples for using Sparse matrices?

yixuan commented 9 years ago

I haven't written one yet, but will add one later since it should be easy.

skn123 commented 9 years ago

Yixuan; I think a sparse matrix support would be the first step as ARPACK is used primarily with large very large sparse matrices. Let me know if I can help / contribute in anyway. Actually, I like the code that you have shared. Maybe you should contact ggael and add it to the Official Eigen repo?

yixuan commented 9 years ago

Hi @skn123 , I've just updated the code and there is now an example in the README file. Basically you just need to use the class SparseGenMatProd to wrap your sparse matrix, and all other part will be the same as before.

ARPACK-Eigen is still an immature project that has only a few months' development. When I think it's ready, I'll ask Gaël whether it deserves to be merged into Eigen. In either way, you can use this library in a neat way since it is header-only. :-)

skn123 commented 9 years ago

Hi yixuan. Indeed, I did see that. One more point. If we are working with symmetric matrices, do we have to pass the full matrix or an upper triangular view would suffice?

yixuan commented 9 years ago

At the current moment you need to pass the full matrix, if you use the built-in matrix wrappers. However, it's pretty straightforward to define your own matrix operation. See how easy it is to write the include/MatOp/SparseGenMatProd.h file.

Actually, you simply need to change the last line from y.noalias() = mat * x; to y.noalias() = mat.selfadjointView<Eigen::Lower>() * x; to force the matrix to be symmetric.

I'm considering how to make this procedure more intuitive and automatic.

skn123 commented 9 years ago

Actually, I figured that out when I studied the code. I think the input matrix is constant and does not change. Thanks for the code. This is indeed great. Would it be also possible to add a JDQR version of this also? I think the ML community frequently use this or JDQR.

yixuan commented 9 years ago

JDQR seems to be an interesting algorithm. If I understand correctly, it is another method to calculate partial eigen decomposition, similar to ARPACK. Is it correct?

skn123 commented 9 years ago

Yes it is. Apparently, it is more stable (rather more stable than legacy ARPACK). Example would be this https://github.com/jhaberstroh/cs294-davidson

yixuan commented 9 years ago

Great. I'll take a look once I find time. Since JDQR is another issue, I would close this one first.