Open dgrnbrg opened 7 years ago
Hey @dgrnbrg, thanks for this example! I'm actually reasonably close to rewriting the eigenvalue/eigenvector routines for rulinalg
, so this is a useful example to check against. I've been pretty busy with real life the last few weeks, but I hope I should have some more time to look at it again in not too long!
HI! I also face this problem. I plan to port this code: https://github.com/fiji/Jama/blob/master/src/main/java/Jama/EigenvalueDecomposition.java into Rust tomorrow, because I have not found reliable eigen solver for rust which will work stably.
@stiv-yakovenko: thanks for your interest, and your willingness to attempt a port!
I see now that it's been 3 months ago since I said that I would hopefully be able to look at it again soon. Unfortunately that didn't really hold up, so it's good if someone wants to take a look at it.
I took a brief look at the Jama implementation (not the details at this point). I don't agree with their overall API design, though I'm sure their implementation might be very good. In particular, an "Eigenvalue Decomposition" only makes sense for symmetric matrices. For non-symmetic matrices, one rather talks about the "Real Schur Decomposition". The idea I had for the design of eigenvalue algorithms for rulinalg
is as follows:
RealSchurDecomposition
, which furthermore allows one to compute (complex) eigenvalues and (complex) eigenvectors (note that these are different from the Schur vectors)SymmetricEigenDecomposition
(possibly a shorter name...), which furthermore allows computation of (real) eigenvalues and eigenvectors.For complex eigenvalues, we'd have to introduce "HermitianEigenvalueDecomposition" for hermitian matrices. I haven't given much thought to the general complex case yet.
It seems that Jama tries to detect whether the matrix is (exactly) symmetric, and then compute either the symmetric eigenvalue decomposition or the real schur decomposition. I think detecting exact symmetry is a fragile approach: indeed, due to roundoff errors a matrix that is supposed to be symmetric may well have very slightly different elements in the corresponding element in the transpose.
Of course, .eigenvalues()
would be a convenient function to access this functionality without explicitly constructing the decompositions.
I'm certainly open for a discussion on how an eigenvalue API could look like in Rulinalg!
HI! I've finished this dirty work of transpiling and debugging the code. And for simple real-valued assymetric matrix the code works correctly. I propose to do some more testing, it may fail in some other cases. Now you can improve the code (I am rust beginner) and orchestrate it the way you like: https://gist.github.com/stiv-yakovenko/2c59c9588b615556f20d4a2537f6cd8d
Thank you for the effort, @stiv-yakovenko! Unfortunately I can't promise that I'll get to it anytime soon (though I hope so, but I'm currently spending my limited time on a different, but tangentially related effort at the moment). If I don't, perhaps someone else is willing to pick it up in the meantime though!
Hey @stiv-yakovenko and @Andlon: if it makes sense, I can provide code review & suggestions for the algorithm.
Well, for my own purposes the algorithm works good enough, I don't have time to integrate it into rulinalg. If, as far as I understand, this project has no funding, then sorry, I'll leave integration task for the community. BTW: algorithm also finds imaginary parts of eigen numbers.
@stiv-yakovenko - thanks for sharing your code. This gives a good start point for somebody looking to integrate into rulinalg.
@dgrnbrg please feel free to take the lead on that work if you want to. I'm sure it will be quite tricky but I am happy to try and provide feedback along the way.
Unfortunately I am in the same boat as @Andlon . I have a lot of other things going on and no time to dedicate for any serious attention to rulinalg for now...
By accident, my program which is using
eigenvalues()
has become very good at constructing matrices which never converge during the search for their eigenvalues. Below is an example of such a matrix--I can actually generate an infinite number of these. I've noticed that this only occurs for matrices that are around 50x50--my smaller matrices generated through the same procedure never have this issue: