Diagonalization of ill-conditioned matrix and imposibility to compute eigenvectors. Different results with numpy/scipy

Question

Im am dealing with a sparse Matrix with very small elements. Consider a vector:

vec=[-1.e-76 -1.e-72 -1.e-68 -1.e-64 -1.e-60 -1.e-56 -1.e-52 -1.e-48 -1.e-44
-1.e-40 -1.e-36 -1.e-32 -1.e-28 -1.e-24 -1.e-20 -1.e-16 -1.e-12 -1.e-08
-1.e-04 -1.e-02 -1.e-04 -1.e-08 -1.e-12 -1.e-16 -1.e-20 -1.e-24 -1.e-28
-1.e-32 -1.e-36 -1.e-40 -1.e-44 -1.e-48 -1.e-52 -1.e-56 -1.e-60 -1.e-64
-1.e-68 -1.e-72 -1.e-76]

For those interested, those numbers represent the hopping amplitudes of a 1D system. They are not zero. The hamiltonian is given by a sparse matrix:

H=sps.diags([vec,vec],[-1,1],dtype='f8')

I am interested on the eigenvalues, but even more on the eigenvectors . As far as I know, there are two ways of deal with the diagonalization: scipy.linalg and numpy.linalg and the former is better.

 denseHam=H.toarray()

The correct eigenvalue spectrum is given by all of these functions:

import numpy as np
import scipy.linalg as la
s1= la.eigvalsh(denseHam)
s2= np.linalg.eigvalsh(denseHam)
s3= np.linalg.eigvals(denseHam) #I did not expect that!

The correct spectrum is:

spectrum=[-3.16230928e-03 -3.16227766e-08 -3.16227766e-13 -3.16227766e-18
-3.16227766e-23 -3.16227766e-28 -3.16227766e-33 -3.16227766e-38
-3.16227766e-43 -3.16227766e-48 -3.16227766e-53 -3.16227766e-58
-3.16224604e-63  3.16224604e-63  3.16227766e-58  3.16227766e-53
 3.16227766e-48  3.16227766e-43  3.16227766e-38  3.16227766e-33
 3.16227766e-28  3.16227766e-23  3.16227766e-18  3.16227766e-13
 3.16227766e-08  3.16230928e-03]

Nevertheless, the other functions (which involve the computation of the eigenvectors also) fail, and I can't go on because I need the eigenvectors.

I have to say that C++ is able to compute correctly also the eigenvectors.

So I have two questions:

1. Why the function np.linalg.eigh(denseHam) gives different spectrum than np.linalg.eigvalsh(denseHam)?
2. Is there any way to compute correctly the eigenvectors with python?

Thank you very much in advance!

--- UPDATE------ I paste here a minimal complete example. Note the exporeous degeneracy of the numpy.linalg.eigh:

import numpy as np
import scipy.sparse as sps

vec=np.array([-1.e-76, -1.e-72, -1.e-68, -1.e-64, -1.e-60, -1.e-56, -1.e-52,
       -1.e-48, -1.e-44, -1.e-40, -1.e-36, -1.e-32, -1.e-28, -1.e-24,
       -1.e-20, -1.e-16, -1.e-12, -1.e-08, -1.e-04, -1.e-02, -1.e-04,
       -1.e-08, -1.e-12, -1.e-16, -1.e-20, -1.e-24, -1.e-28, -1.e-32,
       -1.e-36, -1.e-40, -1.e-44, -1.e-48, -1.e-52, -1.e-56, -1.e-60,
       -1.e-64, -1.e-68, -1.e-72, -1.e-76])
H=sps.diags([vec,vec],[-1,1],dtype='f8')
denseHam=H.toarray()

s1=np.linalg.eigvalsh(denseHam)
(s2,basis)=np.linalg.eigh(denseHam)

print("Note the difference between the eigenvalues computed with eigvalsh (1stcolumn) and eigh (2nd column)")
for elem in range(len(s1)):
    print (s1[elem],"     ",s2[elem])

Answer 1

It is somewhat suprising that they give different results, since I expect both numpy and scipy to call into LAPACK, but this has also generally been my experience.

Note that the scipy bindings offer more arguments to play with; and numpy is likely using different defaults. Some experimentation seems like it is called for; you problem isnt just having very small elements (that could be solved by a simple scaling if it leads to underflow), but your problem is also very 'stiff', with eigenvalues spanning more than 70 orders of magnitude. C++ might give you eigenvectors but I would not be surprised if they were polluted with numerical noise to the point of being useless.

It sounds like the kind of problem where it would be much more reliable to solve it in some kind of transformed/preconditioned space. The docstring doesnt say if the LAPACK functions can process 128-bit float; last time I tried they didnt, but if they do now be sure to use that instead of 64bit at least.

Answer 2

Short answer : try LAPACK's scipy.linalg.lapack.dsyev() !

# LAPACK's dsyev instead of LAPACK's dsyevd
(s3,basis3,error)=scipy.linalg.lapack.dsyev(denseHam)

Functions np.linalg.eigvalsh() and np.linalg.eigh() both call LAPCK's DSYEVD as stated in their documentations. Nevertheless, they deliver different eigenvalues: eigh() fails. The reason is likely engraved in the source code of DSYEVD. Indeed, if the eigenvectors are not needed, the routine scales the matrix, reduces the matrix to trigonal form (DSYTRD) and finally calls DSTERF. This last routine computes all eigenvalues of a symmetric tridiagonal matrix using the Pal-Walker-Kahan variant of the QL or QR algorithm. On the contrary, if eigenvectors are needed, DSYEVD switches to DSTEDC after scaling and reduction to trigonal shape. DSTEDC computes all eigenvalues and eigenvectors of a symmetric tridiagonal matrix using the divide and conquer method.

1. The small norm of the input matrix does not matter, as scaling is likely applied in such a case. Since the real symmetric matrix features eigenvalues of very different magnitude (from 3.16224604e-63 to 3.16230928e-03), it is ill-conditionned. The accuracy of most linear algebra procedures, including eigenvalue computation, is significantly affected by this property. The provided matrix is already in trigonal form.

2. np.linalg.eigh() fails. It is likely related to the use of a divide and conquer strategy.

3. np.linalg.eigvalsh() seems to work. Therefore, it looks like DSTERF worked. Nevertheless, this routine only provides the eigenvalues.

As LAPACK provides different algorithms to compute the eigenvalues and eigenvectors, your C++ code likely uses another one, such as dsyev(). After scaling and reducing the matrix to trigonal form, this routine first call DORGTR to generate the orthogonal matrix, then call DSTEQR to get the eigenvectors. Hopefully, this function can be called through scipy's Low-level LAPACK functions (scipy.linalg.lapack)

I added a few lines to your code to call this function. scipy.linalg.lapack.dsyev() computes eigenvalues similar to those of np.linalg.eigvalsh() for this ill conditionned matrix.

import numpy as np
import scipy.sparse as sps
import scipy.linalg.lapack

vec=np.array([-1.e-76, -1.e-72, -1.e-68, -1.e-64, -1.e-60, -1.e-56, -1.e-52,
       -1.e-48, -1.e-44, -1.e-40, -1.e-36, -1.e-32, -1.e-28, -1.e-24,
       -1.e-20, -1.e-16, -1.e-12, -1.e-08, -1.e-04, -1.e-02, -1.e-04,
       -1.e-08, -1.e-12, -1.e-16, -1.e-20, -1.e-24, -1.e-28, -1.e-32,
       -1.e-36, -1.e-40, -1.e-44, -1.e-48, -1.e-52, -1.e-56, -1.e-60,
       -1.e-64, -1.e-68, -1.e-72, -1.e-76])
H=sps.diags([vec,vec],[-1,1],dtype='f8')
denseHam=H.toarray()

# eigenvalues only, boils down to LAPACK's DSTERF
s1=np.linalg.eigvalsh(denseHam)
# LAPACK's dsyevd, divide and conquer
(s2,basis)=np.linalg.eigh(denseHam)
# LAPACK's dsyev instead of LAPACK's dsyevd
(s3,basis3,error)=scipy.linalg.lapack.dsyev(denseHam)

print("Note the difference between the eigenvalues computed with eigvalsh (1stcolumn) and eigh (2nd column)")
for elem in range(len(s1)):
    print (s1[elem],"     ",s2[elem], s3[elem])

You can also avoid the reduction to trigonal form as your matrix is already trigonal. Scaling is likely required to avoid numerical problems. Then DORGTR and DSTEQR can be called directly through LAPACK functions for Cython. But it requires more work and a compilation step.