Non-eigenvectors from DGEEV

[pav]

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
...


[Julien Langou]

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.

[Julien Langou]

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]

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.

[dmday]

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.

[pav]

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 :)