Non-eigenvectors from DGEEV

Open discussion regarding features, bugs, issues, vendors, etc.

Non-eigenvectors from DGEEV

Postby pav » Wed Apr 21, 2010 3:57 pm

Hi,

Any ideas why for the following matrix: http://projects.scipy.org/numpy/raw-att ... 54/mat.txt

DGEEV returns vectors for which the relative error |A v - lambda v|/||A||_F is of the order of 1e-2? Test program here: http://projects.scipy.org/numpy/raw-att ... 4/test.f90

And so, with LAPACK 3.2.1,

Code: Select all
$ gfortran -o test test.f90 -llapack -lblas
$ ./test < mat.txt
  4.09736095556244961E-026
  4.43089563288366399E-023
   3.4612399980026800      LARGE ERROR
   8.9078382202313087      LARGE ERROR
   8.8674244567312162      LARGE ERROR
   8.9066982533988046      LARGE ERROR
  3.83826475736065412E-032
...
pav
 
Posts: 5
Joined: Tue Sep 16, 2008 2:58 pm

Re: Non-eigenvectors from DGEEV

Postby Julien Langou » Wed Apr 21, 2010 8:14 pm

That's a really cool nasty little matrix that you got there! Congratulations ...
Numpy does not work: http://projects.scipy.org/numpy/ticket/1454
Octave does not work.
Matlab does not work.
LAPACK does not work.
Thanks for reporting this

Code: Select all
clear
A=load('~/Downloads/mat.txt', 'ascii');               
n = 28;
[V,D]=eig(A);                                         
[oo,I] = sort( diag(D) ); V = V(:,I); D = D(I,I);
relgap = zeros(n,1);
relgap(1) = (D(2,2)-D(1,1))/abs(D(1,1));
for i=2:27, relgap(i) = min( [ ((D(i,i)-D(i-1,i-1)))   (D(i+1,i+1)-D(i,i)) ])/abs(D(i,i)); end
relgap(n) = (D(n,n)-D(n-1,n-1))/abs(D(n,n));
fprintf('    eig                 relerr       relgap\n');
for i=1:28, fprintf('%2d %20.16f %e %e\n',i,D(i,i),norm(A*V(:,i)-V(:,i)*D(i,i))/norm(A),relgap(i));, end;


Code: Select all
    eig                 relerr       relgap
 1  -3.1751079149941352 2.033247e-02 4.195976e-16
 2  -3.1751079149941339 2.445700e-02 1.398659e-16
 3  -3.1751079149941335 2.161849e-02 1.398659e-16
 4  -3.1751079149941308 3.274260e-13 8.391952e-16
 5  -2.0498377961923029 1.524756e-02 1.299876e-15
 6  -2.0498377961923002 1.524127e-02 2.166460e-16
 7  -2.0498377961922998 1.468891e-11 2.166460e-16
 8  -2.0498377961922971 1.499134e-02 1.299876e-15
 9  -0.1887301587301591 3.198971e-18 1.176519e-15
10  -0.1887301587301589 1.525037e-18 0.000000e+00
11  -0.1887301587301589 1.664842e-18 0.000000e+00
12  -0.1887301587301586 9.324637e-19 1.617714e-15
13  -0.0495989036825447 1.161271e-03 2.937903e-15
14  -0.0495989036825445 1.165823e-03 1.678802e-15
15  -0.0495989036825444 1.161071e-03 1.678802e-15
16  -0.0495989036825441 1.176836e-03 6.155606e-15
17  -0.0204389846397406 1.378855e-03 3.394931e-15
18  -0.0204389846397405 1.380317e-03 2.376452e-15
19  -0.0204389846397405 1.376691e-03 2.376452e-15
20  -0.0204389846397403 1.377147e-03 7.978087e-15
21  -0.0103199818359591 1.461548e-03 2.151599e-14
22  -0.0103199818359589 1.459206e-03 1.008562e-14
23  -0.0103199818359588 1.466116e-03 1.008562e-14
24  -0.0103199818359586 1.462142e-03 1.512843e-14
25  -0.0077181030285824 1.488029e-03 2.034081e-14
26  -0.0077181030285823 1.480665e-03 9.327554e-15
27  -0.0077181030285822 1.469812e-03 9.327554e-15
28  -0.0077181030285821 1.485513e-03 1.562084e-14


So a lots of cluster of eigenvalues. But, still DGEEV should do its job.
I also looked at "spy(A)" and observed a great structure in the matrix.

Yes, there is a bug in DGEEV. Thanks for reporting it. We will investigate.

Best wishes,
Julien.
Last edited by Julien Langou on Thu Apr 22, 2010 1:06 pm, edited 1 time in total.
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Non-eigenvectors from DGEEV

Postby Julien Langou » Wed Apr 21, 2010 8:27 pm

In Matlab doing:
Code: Select all
% method 1
[Q,H]=hess(A);                                         
[V,D]=eig(H);
V = Q*V;

instead of
Code: Select all
% method 2
[V,D]=eig(A);

removes the problem.

(If Matlab does what I think it does ...)

Method 2 (eig) calls DGEEV directly. So that means: (1) balancing on A, (2) reduce to Hessenberg form, (3) DGEHRD, etc.
Method 1 does: (1) reduce to Hessenberg form, (2) balancing on H, (3) DGEHRD, etc.

Method 1 works, method 2 does not.

Julien
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Non-eigenvectors from DGEEV

Postby Julien Langou » Thu Apr 22, 2010 10:45 am

Following on the previous clue that it might be the balancing the guilty part, it appears that balancing can be turned off in Matlab.
[V,D]=eig(A,'nobalance');

works great.
Julien.
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Non-eigenvectors from DGEEV

Postby dmday » Thu Apr 22, 2010 12:05 pm

To whom it may concern,
Matlab worked for me. 7.9.0.529 (R2009b)
The matrix is, modulo absolute perturbations the size of the machine precision, a Kronecker product.
A = numpy54;
amax = max( max( abs(A))), % approximately 122
A11 = A(1:7,1:7);
I = eye(4,4);
norm( kron(I,A11) - A )/amax, % approximately eps
[V11,D11] = eig(A11);
V = kron(I,V11);
D = kron(I,D11); % [V,D] = eig(A);
Each eigenvalue has multiplicity 4. The condition number of V
is sensitive to the algorithm used. Balance reduces A11 to
Schur (lower triangular) form.
dmday
 
Posts: 2
Joined: Mon Dec 03, 2007 6:18 pm

Re: Non-eigenvectors from DGEEV

Postby pav » Thu Apr 22, 2010 12:09 pm

Julien Langou wrote:That's a really cool nasty little matrix that you got there! Congratulations ... Yes, there is a bug in DGEHRD. Thanks for reporting it. We will investigate.


Very good, thanks! (Though the honor of finding the matrix goes to the person who reported this bug for Numpy :)
pav
 
Posts: 5
Joined: Tue Sep 16, 2008 2:58 pm

Re: Non-eigenvectors from DGEEV

Postby Julien Langou » Thu Aug 15, 2013 10:26 am

This problem is referenced as BUG0057 in our bug list: http://www.netlib.org/lapack/bug_list.html.

The same bug was recently (March 28th, 2013) reported on the forum, see topic 4270: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=4270 and has been reported as BUG0106.

Both bugs (BUG0057 and BUG0106) have the same origin and we fixed the problem in our SVN repository. (See SVN revision 1413 on May 26th, 2013.) We will release the fix in LAPACK 3.5.0.

Cheers,
Julien.
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA


Return to User Discussion

Who is online

Users browsing this forum: Google [Bot] and 1 guest