Search code examples
pythonnumpyscipylinear-algebra

Algebraic value sorting of eigsh()

I am using scipy.sparse.linalg.eigsh() to find the first 10 smallest eigenvalues of a large, sparse and Hermitian matrix H. I'm finding the eigenvalues following the documentation

For which = ‘LM’ or ‘SA’:

If return_eigenvectors is True, eigenvalues are sorted by algebraic value.

If return_eigenvectors is False, eigenvalues are sorted by absolute value.

The eigenvalues I get for H do not seem to be in any apparent order. For instance, if I am not returning eigenvectors

E1 = H.eigsh(k=10,which='SA',maxiter=1E4,return_eigenvectors=False)

gives

E1: [-5.8445426 -5.87106398 -5.87106398 -5.88635936 -6.37709649 -6.38166066 -6.39249619 -6.39249619 -6.39756682 -6.42135482]

If I am returning eigenvalues via E2,V = H.eigsh(k=10,which='SA',maxiter=1E4,return_eigenvectors=True) or equally E2, V = H.eigsh(k=NEvals, which='SA')

I get,

E2: [-6.42135482 -6.39756682 -6.39249619 -6.38166066 -6.37709649 -6.39249619 -5.88635936 -5.87106398 -5.8445426 -5.87106398]

Questions

  1. In E1 the eigenvalues are in ascending order. Why do I need to sort the eigenvalues again after I specified 'SA' to get the smallest eigenvalue?

  2. Why are the eigenvalues in E2 not in any order? Or better, what does it mean to sort by algebraic value? Do I need to sort the eigenvalues again or is there reason (perhaps to do with the complex eigenvectors) as to why they're in this order?

  3. Finally, where I can learn about how (if at all) the eigenvectors affect the order of the eigenvalues? Any explanation with a physics application, where finding the ground state corresponds to the finding the lowest eigenvalue, is also appreciated.

I looked for documentation of eigsh.

Solution

  • I remember discussing this issue on github. The fact that the return order depends on eigenvector is surprising, there is no obvious reason for that, but at least it is documented. As of why exactly, this is how ARPACK works internally, I did not dig into it.

    Now your case abides by the documentation. For case E1 with return_eigenvector=False, you get abs(E1[i]) < abs(E1[i+1]). For case E2 with return_eigenvector=True, you have E2[i]) < E2[i+1]). Both cases are sorted, but with a different sort depending on return_eigenvector. The eigenvectors are sorted according to the eigenvalues. You can always sort them with a different conventions that suits you better:

    E2, V = H.eigsh(k=NEvals, which='SA')
so = np.abs(E2).argsort()
E2 = E2[so]
V = V[:, so]

    This is not related to the properties of the input matrix, just a convention in the function.