Bug in Matlab's Eig and lapack's dgeev

Post here if you want to report a bug to the LAPACK team

Bug in Matlab's Eig and lapack's dgeev

Postby RonMorgan » Thu Mar 28, 2013 8:40 pm

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
RonMorgan
 
Posts: 1
Joined: Thu Mar 28, 2013 7:14 pm

Re: Bug in Matlab's Eig and lapack's dgeev

Postby rodney » Fri Mar 29, 2013 2:03 pm

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
rodney
 
Posts: 48
Joined: Thu Feb 10, 2011 8:20 pm
Location: Colorado College

Re: Bug in Matlab's Eig and lapack's dgeev

Postby eem2314 » Sun Apr 21, 2013 11:46 pm

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?
eem2314
 
Posts: 1
Joined: Sun Apr 21, 2013 11:23 pm

Re: Bug in Matlab's Eig and lapack's dgeev

Postby Julien Langou » Sat Feb 15, 2014 10:01 pm

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.
Julien Langou
 
Posts: 734
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA


Return to Bug report

Who is online

Users browsing this forum: No registered users and 0 guests