cupyx.scipy.sparse.linalg.eigsh yields wrong results

Created on 4 Apr 2021 · 3Comments · Source: cupy/cupy

Conditions
- CuPy version: 9.0.0b3
- OS/Platform: Windows-10
- CUDA version: 11010
Code to reproduce

from cupyx.scipy.sparse import diags
from cupyx.scipy.sparse.linalg import eigsh
from cupyx.scipy.sparse import eye
from cupyx.scipy.sparse import kron

N = 50
T =  diags([-2. , 1. , 1. ], [0,-1, 1] , shape=(N, N))
I = eye(N)
H = kron(T,I) + kron(I,T)

eigenvalues, eigenvectors = eigsh(H, k = 20, which='LA' )
print(eigenvalues)

Out:

[2125.83508304 2126.78382114 2127.81013589 2128.15849898 2128.60714843
 2132.69078101 2133.04216459 2141.02943555 2143.24419361 2143.41556547
 2147.26750965 2187.31610205 2202.85676903 2204.06317034 2204.36118675
 2206.04539525 2208.65665711 2208.90268838 2223.28317277 2282.71630961]

These aren't the correct eigenvalues. These values also vary each time the function is executed. The ones computed using scipy are the correct:

from scipy.sparse import diags
from scipy.sparse.linalg import eigsh
from scipy.sparse import eye
from scipy.sparse import kron

N = 50
T =  diags([-2. , 1. , 1. ], [0,-1, 1] , shape=(N, N))
I = eye(N)
H = kron(T,I) + kron(I,T)

eigenvalues, eigenvectors = eigsh(H, k = 20, which='LA' )
print(eigenvalues)

Out: (correct eigenvalues)

[-0.12081226 -0.10927498 -0.10927498 -0.09790934 -0.09790934 -0.09445993
 -0.09445993 -0.07556511 -0.07556511 -0.0681076  -0.06419947 -0.06419947
 -0.04921278 -0.04921278 -0.03784714 -0.03784714 -0.03031796 -0.01895232
 -0.01895232 -0.00758669]

One additional note:
Replacing;
T = diags([-2. , 1. , 1. ], [0,-1, 1] , shape=(N, N))
by
T = diags([2. , -1. , -1. ], [0,-1, 1] , shape=(N, N))
yields always a correct result:

[7.87918774 7.89072502 7.89072502 7.90209066 7.90209066 7.90554007
 7.90554007 7.92443489 7.92443489 7.9318924  7.93580053 7.93580053
 7.95078722 7.95078722 7.96215286 7.96215286 7.96968204 7.98104768
 7.98104768 7.99241331]

but of course, it isn't the same matrix.

bug high

Source