Resurrect Dhillon's algorithm for symmetric eigen-decomposition?

[Darío de México <da...@gmail.com>]

Hi everyone,

>From version 2.0 => 2.1, we seem to have replaced the algorithm for
calculating the eigen-factorization of symmetric matrices.
You guys were previously using this specialized algorithm from Mr. Dhillon
(plus other stuff):

"This implementation is based on Inderjit Singh Dhillon thesis A New O(n2)
Algorithm for the Symmetric Tridiagonal Eigenvalue/Eigenvector Problem, on
Beresford N. Parlett and Osni A. Marques paper An Implementation of the
dqds Algorithm (Positive Case) and on the corresponding LAPACK routines
(DLARRE, DLASQ2, DLAZQ3, DLAZQ4, DLASQ5 and DLASQ6)."

Javadoc:
https://commons.apache.org/proper/commons-math/javadocs/api-2.0/org/apache/commons/math/linear/EigenDecompositionImpl.html

but mysteriously, you changed the algorithm on version 2.1; and that seems
to have started the trend up to current version 3.6:

"This implementation is based on the paper by A. Drubrulle, R.S. Martin and
J.H. Wilkinson 'The Implicit QL Algorithm' in Wilksinson and Reinsch (1971)
Handbook for automatic computation, vol. 2, Linear algebra,
Springer-Verlag, New-York"

https://commons.apache.org/proper/commons-math/javadocs/api-2.1/org/apache/commons/math/linear/EigenDecompositionImpl.html
https://commons.apache.org/proper/commons-math/javadocs/api-3.6/org/apache/commons/math3/linear/EigenDecomposition.html

I have tested the version 3.6 and 2.0 (with some manual patches you
published for 2.0), and the difference is quite significant. For a
symmetric matrix of 867x867, the times are as follows on my laptop:

3.6: around 14 secs
2.0: around 3 secs!

I could send a test-case that shows the concrete time differences in a
reproducible way, and probably some measures of how different/equal is the
result produced by each algorithm.  But before going there, wanted to hear
your thoughts: would it make sense to bring back the 2.0 algorithm the for
symmetric case? I could contribute with the patch, but wanted to understand
on the first place the reasons for the algorithm change; perhaps there were
good reasons for doing that, which I am not aware of.

Thanks in advance.

PS: I initially posted the suggestion on JIRA, MATH-1334, but they
recommended me to bring the discussion here.