Use ARPACK to implement eigs? (wrong eigenvalues generated by linalg.eigh/DSYEVR) (Trac #1159) #1686

Closed
scipy-gitbot opened this Issue Apr 25, 2013 · 13 comments

Projects

None yet

2 participants

@scipy-gitbot

Original ticket http://projects.scipy.org/scipy/ticket/1159 on 2010-04-21 by trac user ysr, assigned to unknown.

Wrong eigenvalues are generated by the following example. I have to attach the example matrix in order to reproduce the error.

import scipy.linalg
import scipy.io
L = scipy.io.loadmat('Lmat.mat')['L']
# let us compute the first (smallest) 5 eigenvectors
d,v= scipy.linalg.eigh(L,eigvals=(0,4))
# let us check if the eigenvectors are orthogonal
scipy.multiply(v[:,0],v[:,1]).sum()
# answer is zero (to a floating point error). Good!
# now let us compute the first (smallest) 2 eigenvectors
d,v= scipy.linalg.eigh(L,eigvals=(0,1))
# let us check if the eigenvectors are orthogonal
scipy.multiply(v[:,0],v[:,1]).sum()
# they are not orthogonal!! the second eigenvector was not calculated correctly!
@scipy-gitbot

Attachment added by trac user ysr on 2010-04-21: Lmat.mat

@scipy-gitbot

@pv wrote on 2010-04-21

It's a bug in the underlying LAPACK library and affects all programs depending on its DSYEVR implementation. This could/should be reported to LAPACK authors (see http://www.netlib.org/lapack/).

A Fortran test program reproducing the issue is attached.


The issues of concern for Scipy are:

  • should we provide an interface at all to this routine, if the underlying library is buggy
  • are there cheap ways to detect or work around this bug
@scipy-gitbot

Attachment added by @pv on 2010-04-21: test.f90

@scipy-gitbot

Attachment added by @pv on 2010-04-21: L.txt

@scipy-gitbot

@pv wrote on 2010-04-21

(Note that similar conclusions are very likely to apply for any bugs discovered in Scipy's any other linear algebra routines: the Python interface is very light, and so most bugs are likely to lie in the underlying library.)

@scipy-gitbot

@pv wrote on 2010-04-21

Ok, I forwarded it to LAPACK people.

Leaving this bug report open, though -- we might want to add some cautionary words into the docstring. I note that Octave always uses ARPACK for eigs, and it might be more stable in this respect.

@scipy-gitbot

trac user ysr wrote on 2010-04-22

Thanks for coding it up in fortran and forwarding it to the Lapack people!

Yes, I think arpack does not have this problem (since MATLAB works correctly on it and they use the Arnoldi Iterations). I will try it in Octave. There should be an eigs/eigh built on top of arpack in scipy. Can scipy 0.7.2 get this feature? Computing wrong eigenvalues is pretty horrific for a linear algebra library.

@scipy-gitbot

@pv wrote on 2010-04-22

ARPACK is available as scipy.sparse.linalg.eigen.arpack.eigen (in Scipy 0.7.1), but not used by the eigs-type routines at the moment. Perhaps it should be, at least optionally, since ARPACK is probably better tested as computing a few eigenvalues is more or less its main purpose.

(However, for some reason I often get exceptions when calling the ARPACK -- probably needs a separate bug report for that.)

I think it's unlikely at the moment that there will be any change in 0.7.2 this late in the release schedule. 0.8 is a more realistic target for this...

@scipy-gitbot

Title changed from Wrong eigenvalues generated by scipy.linalg.eigh to Use ARPACK to implement eigs? (wrong eigenvalues generated by linalg.eigh/DSYEVR) by @pv on 2010-04-22

@scipy-gitbot

trac user ysr wrote on 2010-04-22

Interestingly, scipy.sparse.linalg.eigen.arpack.eigen gives similar issues, but scipy.sparse.linalg.eigen.arpack.eigen_symmetric seems to give the correct answers. Octave's eigs also has issues for this matrix when we compute the smallest eigenvectors. MATLAB seems to give the correct answers.

@scipy-gitbot

@pv wrote on 2010-12-11

The same bug appears also in the DSYEVX routine, so switching to it does not help.

@tonysyu tonysyu pushed a commit to tonysyu/scipy that referenced this issue May 9, 2013
@rgommers rgommers BUG: skip unknown wavfile chunks when reading; still return the data.
Closes #1585, #1626 and #1686.  Thanks to richli for the patch.
c74b315
@pv
Member
pv commented Aug 18, 2013

The LAPACK bug bug0056 is still there.

@pv
Member
pv commented Dec 10, 2015

Fixed in lapack 3.6.0

@pv pv closed this Dec 10, 2015
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment