scipy.linalg.sparse.eigsh returns negative and non-consistent eigenvalues for positive semi-definite matrix

Question

I try to use scipy.linalg.sparse.eigsh (let's call it method 1 : M1) to compute the smallest eigenvalues of the Laplacian matrix of a real symmetric semi-definite matrix W.

As a benchmark, I ran the computation against scipy.linalg.eigh (method 2 : M2). It seems M1 returns different eigenvalues from M2, and moreover thoses eigenvalues seems to be wrong ones.

Here is the code :

# Initialization
N=10
np.random.seed(0) # for reproducibility 
Nvp=4

# Create a real symmetric semi-definite matrix W, precision float64
X=np.random.random((N,N))
W=np.dot(X, X.T) 
W=np.array(W, dtype=np.float64) 

# Compute its Laplacian matrix A, precision float64
d=np.array([sum(W[:][j]) for j in range(N)], dtype=np.float64)
D=np.diag(d)
A=D-W
for i in range(N):
    for j in range(N):
        A[i][j]/=np.sqrt(d[i], dtype=np.float64)*np.sqrt(d[j], dtype=np.float64)
A=np.array(A, dtype=np.float64)

Let's check A is correctly formated :

>>> A.dtype
dtype('float64')
>>> np.allclose(A, A.T)
True

Now let's run some tests :

## 1) Compute A's smallest eigenvalues by 2 different means
wA2, vA2 = la.eigh(A)
wA1, vA1 = sparse.eigsh(A, k=Nvp, sigma=0, which='LM')

## 2) Compute W's smallest eigenvalues by 2 different means
wW2, vW2 = la.eigh(W)
wW1, vW1 = sparse.eigsh(W, k=Nvp, sigma=0, which='LM')

# Output computed eigenvalues
print(wA2[:Nvp])
print(wA1[:Nvp])
print(wW2[:Nvp])
print(wW1[:Nvp])

Here is the output :

>>>[-1.88737914e-15  9.03999970e-01  9.23513143e-01  9.52678423e-01]

[-4.93242313e-01 -8.14381749e-17  9.22235466e-01  9.44848237e-01]

[0.00575077 0.04770667 0.08565863 0.16319798]

[0.00575077 0.04770667 0.08565863 0.16319798]

This first output shows M1 computes a negative eigenvalue for A, which is not mathematically possible. Moreover, if one checks the computation of the others, let's say the 3rd one :

>>> np.dot(A, vA2[:,2])-wA2[2]*vA2[:,2]
array([-0.01183104, -0.25059123,  0.47372558, -0.31358618, -0.2968036 ,
   -0.04832199,  0.40973325, -0.01369472,  0.33267859, -0.27122678])

It is not even close to zero. I must add the computed eigenvalues differ each time. To me it is due to the initialization vector. I would say scipy.linalg.sparse.eigsh does not iterate enough to get close enough to the real result, but setting maxiter=1000000 does not affect the strange results. Concerning the negative eigenvalue, I'm unfortunately clueless.

I'm running :

Python 3.7.0 (default, Jun 28 2018, 13:15:42)

[GCC 7.2.0] :: Anaconda, Inc. on linux

NumPy and SciPy have been built against Intel MKL.

Can anyone enlighten me ? Thank you in advance for your time.

Hm, looks strange. If no one answers here, you should consider posting this as an issue to scipy (github.com/scipy/scipy/issues) — TheWaveLad
"...compute the smallest eigenvalues..." Instead of sparse.eigsh(A, k=Nvp, sigma=0, which='LM'), have you tried sparse.eigsh(A, k=Nvp, which='SM')? That worked for me. — Warren Weckesser
Hi Warren, thank you for your answer. Indeed this solution works for me as well, but the very point I chose sparse.eigsh is to use the shift-inverted mode feature, which is way faster than other features, and which requires to specify sigma=0 and which='LM'. I refer you to SciPy documentation [docs.scipy.org/doc/scipy/reference/tutorial/arpack.html] for more details. — Nero

Warren Weckesser Warren Weckesser · Accepted Answer · 2018-11-07T14:50:32

Your matrix A is singular, i.e. 0 is an eigenvalue of A. The value of sigma in the shift-invert method must not itself be an eigenvalue. Take a look at the section "Shift and Invert Spectral Transformation Mode" of the ARPACK documentation.

The same formulas are shown in the SciPy tutorial. The formula for the shift-invert method is

inv(A - sigma*M) @ M @ x = nu*x

(Argh, I really wish stackoverflow had LaTeX markup!) When sigma is 0, A - sigma*M is just A, and inv(A) does not exist.

scipy.linalg.sparse.eigsh returns negative and non-consistent eigenvalues for positive semi-definite matrix

1 Answers