I am thinking about adding an OpenCL implementation for symmetric eigendecomposition into Math. Now I have the dilemma of which algorithm to choose.
To quickly summarize there are usually two steps to the process. The first one is tridiagonalization, for which the obvious choice is block Housholder algorithm. For the second step there are multiple options. The slowest and most numerically accurate is QR iteration (this one is used in Eigen). Divide & conquer is almost as accurate, but faster. Multiple relatively robust representations (MRRR) is the fastest by a large margin, but the least accurate.
Once I already implemented MRRR and attempted to put it in Math (for anyone interested in history: https://github.com/stan-dev/math/pull/1159). There were some concerns about MRRR’s accuracy, so I am willing to implement another algorithm if people agree that would be better. Alternatively I could reimplement MRRR with increased precision (for example using double double scalars), which would probably still be faster than other algorithms, but I am not 100% sure if it would solve accuracy issues in all edge cases.
What about separate functions like eigenvalues_sym_mrrr and eigenvectors_sym_mrrr?
Either that or making eigenvalues_sym default to the fast MRRR so that by default people get the fast thing. Then we just expose eigenvalues_sym_qr or whatever as separate functions which use the slower Eigen algorithms.
I don’t have a sense of where eigensolvers are used directly in models. Is there a particular thing where this came up or is this generically improving performance in Stan Math? I did a model with an eigensolve once before and ended up using a separate library.
Yes, I think we still have no way to return two matrices, so we need separate functions for eigenvalues and eigenvectors.
Either that or making eigenvalues_sym default to the fast MRRR so that by default people get the fast thing. Then we just expose eigenvalues_sym_qr or whatever as separate functions which use the slower Eigen algorithms.
I don’t really want to implement two different solvers … They are not exactly simple.
Is there a particular thing where this came up or is this generically improving performance in Stan Math?
No, I don’t have a particular example. I am just trying to improve performance in general.