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.

So â€¦ any thoughts? @Erik_Strumbelj @rok_cesnovar @Bob_Carpenter @stevebronder @bbbales2 @bgoodri