Bug in Matlab's Eig and lapack's dgeev

[RonMorgan]

I have found that in some of my restarted Krylov programs a bug in the Matlab Eig command is causing significant problems. With help from Mark Embree, Dan Sorensen and Charles Puelz, it appears that the bug comes from LAPACK's dgeev. Apparently from posts on the topic "Non-eigenvectors from DGEEV", something is already known of this bug, but I am getting problems routinely, instead of for a "nasty little matrix". As the program proceeds, some eigenvectors become very accurate. This causes very small entries in the lower triangular portion of the H matrix that Eig is applied to. Then some of the eigenvectors start losing accuracy and can become very inaccurate. I have been dealing with the problem by "locking in" eigenvectors (zeroing part of H), but this can cause loss of accuracy for eigenvalues still converging if applied too soon.
Mark Embree has simplified it to this matrix which shows the problem well:
A = [1 1 0 0; 0 2 1 0; 0 0 3 1; eps^2 0 0 4], where eps is machine epsilon.
Eig gives one inaccurate eigenvector for this matrix. Also, if eps^2 is replaced with eps^4, there are two inaccurate eigenvectors. So very small perturbations cause very bad results.
In Matlab, using 'nobalance' seems to fix the problem, as does using schur before eig. Using hess does not always work. However, even though there are options to fix it, it seems that Matlab and LAPACK would want this problem eliminated.
Thanks,
Ron Morgan

[rodney]

Ron,

What I think you are seeing is one of the rare cases where balancing seriously degrades the accuracy of the computed eigenvectors. For a reference, see: David S. Watkins, A case where balancing is harmful, Electron. Trans. Numer. Anal., 23 (2006), pp. 1-4. which is available on his website at http://www.math.wsu.edu/faculty/watkins/pdfiles/balbad.pdf. Watkins writes in his conclusion:

"Balancing sometimes seriously degrades accuracy. In particular, one should not balance a matrix after it has been transformed to Hessenberg form. However, we must emphasize that balancing is usually not harmful and often very beneficial. When in doubt, balance."

The solution is as you discovered to use the 'nobalance' option in Matlab. In LAPACK, you can use dgeevx instead of dgeev to optionally turn off balancing. I tried your matrix using a modified version of dgeev with balancing turned off and it worked fine.

--Rodney

[eem2314]

It seems to me the real problem is that *geev calls *gebal with 'B' instead of 'P' as is done in *ggev for the generalized problem.
This means the matrix is permuted and scaled in *geev instead of of just permuted. No mention is made in the doc about
scaling - that is left to *geevx so is this not a bug?

[Julien Langou]

Hi Ron, thanks for the bug report. We (Rodney, Brad and I) wrote a technical report about this. See: http://arxiv.org/abs/1401.5766. We developed a modified balancing routine that we released in LAPACK 3.5.0. (See arxiv report for some explanation.) LAPACK 3.5.0 (see: http://www.netlib.org/lapack/lapack-3.5.0.html) is just fine with your matrix now. Cheers, Julien.