Sparse Eigensolver issue.

[Ericgig]

Bug Description

The sparse eigensolver can fail when eigenvalues are purely imaginary.

Code to Reproduce the Bug

import qutip 
(qutip.num(4) * 1j).eigenenergies(sparse=True)

Code Output

array([-3.33066907e-16+3.j, -1.06901857e-16+1.j,  1.03300843e-16+1.j,
        4.44089210e-16+3.j])

Expected Behaviour

All unique eigenstates should be returned.

Your Environment

Present in both `master` and `dev.major`.

Additional Context

No response

[pschindler]

I see a similar issue with just calculating the eigenstates of a simple $\sigma_x \otimes \sigma_x$ operator

import qutip as q
import scipy

op = q.tensor(q.sigmax(), q.sigmax())
eigv_sp, eigs_sp = op.eigenstates(sparse=True)
eigv_full, eigs_full = op.eigenstates(sparse=False)
eigv_scipy, eigs_scipy = scipy.linalg.eigh(op.full())
print(eigv_sp, eigs_sp)
print(eigv_full, eigs_full)
print(eigv_scipy, eigs_scipy)

The sparse=True gives eigenstates

[Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket
 Qobj data =
 [[ 0.36754009-0.27845703j]
  [-0.49166296+0.21364338j]
  [ 0.49166296-0.21364338j]
  [-0.36754009+0.27845703j]]

whereas sparse=False and the scipy.linalg.eigh give the correct states

Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket
 Qobj data =
 [[-0.70710678]
  [ 0.        ]
  [ 0.        ]
  [-0.70710678]]

[AGaliciaMartinez]

Hi @pschindler, thanks for reporting this! Although the eigenenstates are different for both cases (sparse=True and sparse=False) I believe that they are both equally valid since the operator has degenerated eigenvalues. The approach for the computation of the eigenstates for both sparse=True and sparse=False is different so that may explain why you see different results.

For the case that @Ericgig presented though, the eigenvalues are non-degenerated. I tracked down the issue and it is related to how we compute these eigenvalues. When we ask for all the eigenvalues what we actually do is to compute the smallest half in with one method (scipy.sparse.linalg(which='LR')) and the largest half with a different method (scipy.sparse.linalg(which='SR')). Smallest half and largest half here means relative to the real part of the eigenvalue. This makes complete sense for the cases where the eigenvalues are real but unfortunately not so much sense when we only care about the imaginary part. This is what the which is doing, is selecting the values with smallest/largest real part.

If we include a small real component the bug completely disappears:

In [56]: (qutip.num(4)*(1j+0.00001)).eigenenergies(sparse=True)
Out[56]:
array([1.23048985e-16-2.21177095e-16j, 1.00000000e-05+1.00000000e+00j,
       2.00000000e-05+2.00000000e+00j, 3.00000000e-05+3.00000000e+00j])

Because now there is a real part to use in the ordering.

The behavior of eigenstates is definitely wrong and one way of solving it would be to not use both LR and SR options when all the eigenvalues are required and instead chose internally just one. This would still return results ordered by their real part but at least all the eigenvalues would be present (only once).

We may want to extend the behavior of eigenstates to not only accept sort=low and sort=high and instead have a similar behavior to the one provided by scipy and allow sorting by imaginary part or absolute value.

[pschindler]

Hi @Ericgig

Thank you for the clarification. You are right for the $\sigma_x \otimes \sigma_x$ operator, the eigenvectors span indeed the same space.

I stumbled upon weird behavior of eigenstates() when calculating the eigenstates for following operator (shift operator on a qudit):

import qutip as q
import scipy

max_n = 6
op = 0
for i in range(max_n - 1):
  op += q.projection(max_n, i, i+1)
op += q.projection(max_n, max_n-1, 0)

eigv_sp, eigs_sp = op.eigenstates(sparse=True)
eigv_full, eigs_full = op.eigenstates(sparse=False)
eigv_scipy, eigs_scipy = scipy.linalg.eigh(op.full())
print(f"Sparse eigenvalues: {eigv_sp}")
print(f"Full eigenvalues: {eigv_full}")
print(f"Scipy eigenvalues: {eigv_scipy}") 

Here I get differing eigenvalues for the sparse=True and the other two methods:

Sparse eigenvalues: [-1.  -0.5 -0.5  0.5  0.5  1. ]
Full eigenvalues: [-1.  0.  0.  0.  0.  1.]
Scipy eigenvalues: [-1.  0.  0.  0.  0.  1.]

[AGaliciaMartinez]

Hi @pschindler,

Indeed, the behavior you show is a bug in qutip 4.7.0. Thanks for reporting this! However, it is actually a bug in the sparse=True case. The scipy case gives a different result compared to qutip sparse because of the method used. eigh assumes the input is hermitian but op is not. You can see the correct eigenvalues using:

In [7]: scipy.linalg.eigvals(op.full())
Out[7]:
array([-1. -2.22044605e-16j, -0.5+8.66025404e-01j, -0.5-8.66025404e-01j,
        0.5+8.66025404e-01j,  0.5-8.66025404e-01j,  1. +0.00000000e+00j])

I am genuinely surprised that the sparse=False is bugged in the latest version of qutip. Fortunately, I tested this in the master branch of qutip (this contains the qutip 5.0.0 alpha version) and the bug seems to be solved there. This operator is indeed subjected to the same bug Eric shows (it is non-hermitian). However, due to the real part of the eigenvalues, it actually shows the correct answer (similar to the example I showed in my previous comment).

I will try to find the source of the error but it seems very likely that we incorrectly label the operator as hermitian and hence uses the wrong method to find eigenvalues.

[pschindler]

Hi @AGaliciaMartinez . Good to hear from you again ;)

It seems that the hermitian check for the operator op evaluates as True

qutip.cy.spmath.zcsr_isherm(op.data) == True

[pschindler]

Also, if I add op.isherm = False to manual set isherm it seems to work:

import qutip as q
import scipy

max_n = 6
op = 0
for i in range(max_n - 1):
  op += q.projection(max_n, i, i+1)
op += q.projection(max_n, max_n-1, 0)
op.isherm = False

eigv_sp, eigs_sp = op.eigenstates(sparse=True)
eigv_full, eigs_full = op.eigenstates(sparse=False)
eigv_scipy, eigs_scipy = scipy.linalg.eigh(op.full())
print(f"Sparse eigenvalues: {eigv_sp}")
print(f"Full eigenvalues: {eigv_full}")
print(f"Scipy eigenvalues: {eigv_scipy}") 

[AGaliciaMartinez]

Hi @AGaliciaMartinez . Good to hear from you again ;)

Hi! 😄

It seems that the hermitian check for the operator op evaluates as True

Yes, this is definitely the source of the error. Nice catch! I will try to compare the 4.7.1 version of qutip with the new pre lease of version 5.0.0 and see if we can provide a quick fix for it.

[hodgestar]

Fixed in #2082, I believe.

[Ericgig]

Nope, the original issue is not fixed.
We could close it as a won't fix, and refuse to compute all eigenvalues when using the sparse eigen solver as scipy does.
But there is place for improvement that I think we should implement before closing.

[gadhvirushiraj]

Hello @AGaliciaMartinez,
I would like to work on this issue. Can you give some details on where to start ?

[Ericgig]

I do not think it's a good first issue...
There are multiple possible solutions and it is not clear which is the best one:

1. Ensure all eigen values and eigen vectors are computed.
   1.1 This could be done by finding an alternative to scipy to solve for sparse matrix decomposition, but it would add extra requirement and it must be easy to install with pip.
   1.2 Use a third call to scipy.sparse.linalg.eigs/ scipy.sparse.linalg.eigsh using sigma to find all the pseudo degenerate values. It need to be done while being careful about special cases. This would make the sparse eigen solver even slower than it already is.
2. Allow to sort by other orders, such as the complex value or magnitude. This would not solve the issue, but give some options to go around it. However the dense eigen solver would be expected to behave the same and the dense eigen solver does not have this options, so we need to implement it ourselves. This sorting option would also be required for any other data type we want to support. To go around that, the Qobj methods could be responsible for the sorting, the sort argument would only be passed to the dispatched function when sparse=True.
3. Refuse to use the sparse solver for the full decomposition as scipy does. The sparse eigensolver is both slow and hard to use in this case. Also the eigen vectors computed by the sparse solver are stored in a dense array, so it is not that useful to save memory...
4. Implement our own sparse eigen decomposition...
5. Only document the issue without fixing. The default is the dense solver, which is both faster and safer in almost all cases.
6. Something else...

I would personally go for either 1.2 of 3.
Adding different sorting options in Qobj methods (2) and removing them from the dispatched functions, could be nice, but it's not a fix.

The place to start would be to debate on the solution. If you have some opinion on the matter, please let us know.
Once we've made up our mind, you can implement it, but only 5 would be an easy fix in line with other good first issue.

[AGaliciaMartinez]

5 seems also fine for me. We default to the dense solver when all the eigenvalues are required, even if sparse=True and we then raise a Warning explaining why we do this and how to stop the warning. IIf I understood it correctly, 3 would be raising an Error instead of a warning, which may be a little bit more annoying for someone that just wants the code to work (?). I am fine with both though and I think they could be suitable issues for a "medium" level "good_first_issue" (if that makes sense 😅 ).

1 and 4 seem indeed too much work right now (but perhaps they can be tackled in the future). 2 may actually be separate an enhacement, that could be useful for some specific scenarios (?). But maybe it is best to wait for when this scenario is present to actually motivate the enhancement (?).

[Ericgig]

By (3), I meant removing the csr specialisation from the dispatched function and call the function directly when not all values are desired. When Qobj.eigenstates is called, the sparse operation being not available, it would defer to the dense code, no warning nor error.
We can't fully trust eigs_csr as it is, so I would like to make it fully working (1) or remove it from the dispatch (3).

Yes (2) can be seen as separate enhancement.

For me (5) was purely documentation, a warning is fine, but when sparse is specified, it would use eigh_csr and too bad if wrong results are obtained...

(3) would also mean removing the

if isinstance(L.data, _data.CSR) and not sparse:
    L = L.to(_data.Dense)
evals, evecs = _data.eigs(L.data)

since they are no longer used (spectrum.py, floquet.py, qobj.py) and reviewing the tests cases.