Why does SciPy eigsh() produce erroneous eigenvalues in case of harmonic oscillator?

Question

I'm developing some larger code to perform eigenvalue computations of huge sparse matrices, in the context of computational physics. I test my routines against the simple harmonic oscillator in one dimension, since the eigenvalues are well known analytically. Doing so and comparing my own routines to SciPy's inbuilt solvers, I have come across the oddity displayed in the plot below. Here you can see the first 100 numerically computed eigenvalues $\lambda_{num}$ and analytical eigenvalues $\lambda_{ana}$

Around eigenvalue number 40, the numercial results start to diverge from the analytical ones. This doesn't surprise me (I won't go into why here, unless it comes up in the discussion). However, what is surprising to me is that eigsh() produces degenerate eigenvalues (around eigenvalue number 80). Why does eigsh() behave like that for even such a small number of eigenvalues?

import numpy as np
from scipy.sparse.linalg import eigsh
import myFunctions as myFunc
import matplotlib.pyplot as plt

#discretize x-axis
N = 100
xmin = -10.
xmax = 10.
accuracy = 1e-5
#stepsize
h = (xmax - xmin) / (N + 1.)
#exclude first and last points since we force wave function to be zero there
x = np.linspace(-10. + h,10. - h,N)
#create potential
V = x**2

def fivePoint(N,h,V):
    C0 = (np.ones(N))*30. / (12. * h * h) + V
    C1 = (np.ones(N)) * (-16.) / (12. * h * h)
    C2 = (np.ones(N)) / (12. * h * h)
    H = sp.spdiags([C2, C1, C0, C1, C2],[-2, -1, 0, 1, 2],N,N)
    return H

H = myFunc.fivePoint(N,h,V)
eigval,eigvec = eigsh(H, k=N-1, which='SM', tol=accuracy)

#comparison analytical and numerical eigenvalues
xAxes = np.linspace(0,len(eigval)-1,len(eigval))
analyticalEigval = 2. * (xAxes + 0.5)
plt.figure()
plt.plot(xAxes,eigval, '+', label=r"$\lambda_{num}$")
plt.plot(xAxes,analyticalEigval, label=r"$\lambda_{ana}$")
plt.xlabel("Number of Eigenvalue")
plt.ylabel("Eigenvalue")
plt.legend(loc=4)
plt.title("eigsh()-method: Comparison of $\lambda_{num}$ and $\lambda_{ana}$")
plt.show()

Answer 1

The degeneracy of some eigenvalues looks to me like the hallmark of the breakdown of the Lanczos algorithm. The Lanczos algorithm is one of the more commonly used methods to approximate the eigenvalues and eigenvectors of Hermitian matrices; it's what scipy.eigsh() uses, through a call to the ARPACK library.

In exact arithmetic, the Lanczos algorithm produces a set of orthogonal vectors, but in floating point arithmetic these can fail to be orthogonal and even become linearly dependent. The really annoying thing is that this loss of orthogonality happens precisely when one of the approximate eigenvalues has converged to one of the real eigenvalues -- the algorithm sabotages itself, so to speak. The result is that you'll get some spurious pairs of nearby eigenvalues. There are various fixes for this, for example using Gram-Schmidt to force any converged eigenvectors to be orthogonal at every step.

Nonetheless, no method is perfect, especially if you're trying to compute the entire spectrum of your matrix. So if you're trying to get the 50 smallest eigenvalues, you may be better off approximating the wave function by a vector with 100 elements and only asking eigsh() for the first 50 energy levels, rather than using a vector with 50 points and asking for all of the eigenvalues.

If you want to read more, look at Yousef Saad's Numerical Methods for Large Eigenvalue Problems.