3 ms·
For large and sparse matrices you can use the (restarted) Arnoldi method to compute a partial Schur decomposition AQ=QR where Q is tall and skinny and R has a f
by stabbles 3y ago
For large and sparse matrices you can use the (restarted) Arnoldi method to compute a partial Schur decomposition AQ=QR where Q is tall and skinny and R has a few dominant eigenvalues on the diagonal (i.e. eigenvalues on the boundary of the convex hull).
MATLAB uses ARPACK's implementation of this when you call `eigs`
I wrote my own implementation ArnoldiMethod.jl in julia, which unlike MATLAB/ARPACK supports arbitrary number types, and also should be more stable in general, and equally fast.
[1] https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl
- sfpotter 3y agoWhy should it be more stable?
- stabbles 3y agoAll Arnoldi methods basically repeatedly grow and shrink a search space until they contain good approximations to Schur vectors. Shrinking (=restart) is necessary to keep constant memory requirements, as well as performance. Restart is somewhat non-trivial because it has to preserve the Arnoldi relation. ARPACK uses Implicit Restart, which is a very beautiful trick based on "perfect shifts": a single iteration of the shifted QR algorithm where the shift is taken as an exact eigenvalue is a way to reorder basis vectors (moving vectors you want to remove to the back, so you can truncate them) while preserving the Arnoldi relation. However, perfect shifts in the QR algorithm are good on paper but can lead to arbitrary loss of precision in finite arithmetic, especially with nearly repeated eigenvalues. I think newer implementations (like SLEPc) are often based on Stewart's Krylov-Schur paper. It's not as pretty but certainly more stable. In short: give ARPACK a matrix with a cluster of almost repeated eigenvalues at the boundary of the convex hull, and it tends to fail.
- sfpotter 3y agoThanks, very helpful. I’ve used ARPACK extensively and have had no complaints with it, but I haven’t tried this particularly pathological case.
- ilayn 3y agoI guess it is a rite of passage to rewrite it. I'm doing it for SciPy too together with Propack in [1]. Somebody already mentioned your repo. Thank you for your efforts. [1]: https://github.com/scipy/scipy/issues/18566 https://github.com/scipy/scipy/issues/18566