Wsrtjtyk

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.

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.