Audit use of eig versus eigh (and eigvals versus eigvalsh)

[rileyjmurray]

While reviewing a PR by Juan I came across a strange comment explaining deliberate use of eig instead of eigh, even though the matrix was Hermitian:
#480 (comment)

We should audit the codebase for these instances and review if those uses are sensible. I'm opening as a separate issue from Juan's PR in order to not slow him down.

[rileyjmurray]

Scipy's eig and eigh dispatch LAPACK's xGEEV and xHEEV/xSYEV respectively.

Question I sent to the LAPACK maintainers:

I have a vague recollection that zHEEV could be less reliable than constructing a Hermtian matrix in "general" format and then running zGEEV. Am I making this up, or is there something to that? If I always had the matrix in general format anyway, what are the pros/cons of these functions?

Response from LAPACK maintainer 1

Responding to
I have a vague recollection that zHEEV could be less reliable than constructing a Hermtian matrix in "general" format and then running zGEEV.

I am confused on what context someone should say this. The reverse would make sense. But this statement makes no sense. So maybe you simply got things reversed when you wrote it. If you know more about this claim, I would be interested the context, a reference, etc. This is certainly not my opinion and surprising to hear.

Asking how bad it is to use ZGEEV on a Hermitian matrix instead of ZHEEV is a little bit like asking how bad it is to use GETRF on Hermitian Positive Definite matrix instead of POTRF.

I think calling ZGEEV instead of ZHEEV for a Hermitian matrix might qualify as sin #1 in using LAPACK, and then calling GETRF for a Hermitian Positive Definite matrix instead of POTRF might be second to worse.

Pros of using ZHEEV instead of ZGEEV (for a Hermitian matrix):

ZHEEV knows the Hermitian structure of A so it will return some real eigenvalues, (ZGEEV will return complex eigenvalues)

ZGEEV will go to Schur form, and then it will have to compute eigenvectors. ZHEEV directly reads eigenvectors from the diagonalization.

ZHEEV will return eigenvectors that are orthogonal. Eigenvectors of close eigenvalues will not be orthogonal with ZGEEV. In particular, if you have a multiple eigenvalues, then that’s an issue. ZGEEV has not idea that the eigenvectors should be orthogonal and the fact that A is Hermitian in input will not help you (when there is a double eigenvalue).

Algorithms to go to reduce form is 4/3 n^3 (ZHEEV) instead of 10/3n^3 (ZGEEV) and this is not only FLOPS, you can also do more optimized algorithm such as successive band reduction.

ZHEEV reduces the complex Hermitian problem to a real arithmetic problem. So a lot of the computation in ZHEEV are done using real arithmetic. ZGEEV will be full on complex arithmetic.

Many good algorithms for ZHEEV (MRRR, D&C, etc.) that are not available for ZGEEV

If you run ZGEEV on a Hermitian matrix, you end up doing non-Hermitian perturbation to the problem, so you lose the Hermitian-ity of the problem.

ZHEEV does some perturbation on the matrix, but these are Hermitian matrix so you do not lose the nice properties of your problem.

Pros of using ZGEEV (for a Hermitian matrix) instead of ZHEEV:

None

When you use ZGEEV for a Hermitian matrix, you do inherit some of the good properties. Some you get their benefices “for free”. For example, cubic convergence of the QR algorithm. That’s nice. And ZGEEV on a Hermitian matrix will get this naturally. However there are a lots of nice properties that (1) you could use to go faster, and (2) returns properties expected by users.

ZGEEV for a Hermitian matrix should return eigenvalues that are more or less real, eigenvectors that are more or less orthogonal (if not multiple eigenvalues). But that’s really not the way to do it.

Response from LAPACK maintainer 2

I agree with [the first peron's] statement. Occasionally in Matlab, I've computed eig of a supposedly Hermitian matrix, and gotten complex eigenvalues out. That prompted me to look more closely to see if the matrix was exactly Hermitian, and found the upper and lower triangle were ever so slightly different due to rounding errors, so Matlab classified it as "general" (geev) instead of "Hermitian" (heev[d]). The end result was complex eigenvalues, which were clearly wrong for a Hermitian matrix.

There are many versions of the Hermitian eig (heev, heevd, heevr, heevx, etc.), which have varying accuracy properties, but all should be better and faster than geev.

My suggested course of action, endorsed by LAPACK maintainers

Based on all your responses, I'm leaning toward going through the project's codebase and changing eig[vals] calls to eig[vals]h calls for matrices that should be numerically symmetric. If I run into issues after doing that, I'll report back. Sound reasonable?

Other thoughts

Iterative eigenvalue solvers, like eigs and eigsh in scipy.linalg.sparse, can produce horribly incorrect output (usually with a warning). Conflating reliability/behavior of eigs and eigsh with the reliability and behavior of eig and eigh might be how eigh got a bad rap.