sparse implementation of `eigenenergies` not sorting according to `which`

[eunjongkim]

Bug Description

Hi, I am trying to find the lowest eigenenergies of a large Hamiltonian by using sparse implementation. However, I believe just specifying which=:SR as a keyword argument doesn't work. Could you help me with this?

Best,
Eunjong

Code to Reproduce the Bug

using QuantumToolbox

# Hilbert Space Dimensions
N_t = 16
N_r = 300

# annihilation operators
a_op = destroy(N_r)
b_op = destroy(N_t)

E_C = 0.22 * 2π # in units of GHz
E_J = 110 * E_C # in units of GHz
n_g = 0.0

φ_zpf_t = (2 * E_C / E_J) ^ (1/4)
n_zpf_t = 1 / (2 * φ_zpf_t)

φ_t = φ_zpf_t * (b_op + b_op')
n_t = -im * n_zpf_t * (b_op - b_op')

H_t = 4 * E_C * (n_t - n_g) ^ 2 - E_J * cos(φ_t)
eig_E_t = eigenenergies(H_t)

ω10 = real(eig_E_t[2] - eig_E_t[1])
ω21 = real(eig_E_t[3] - eig_E_t[2])
α = ω21 - ω10

println("Qubit Transition Frequency $(ω10/(2π)) GHz (Dense eigval)")
println("Anharmonicity $(α/(2π)) GHz (Dense eigval)")


eig_E_t_sp = eigenenergies(H_t; sparse=true, eigvals=5, which=:SR)

ω10_sp = real(eig_E_t_sp[2] - eig_E_t_sp[1])
ω21_sp = real(eig_E_t_sp[3] - eig_E_t_sp[2])
α_sp = ω21_sp - ω10_sp

println("Qubit Transition Frequency $(ω10_sp/(2π)) GHz (Sparse eigval)")
println("Anharmonicity $(α_sp/(2π)) GHz (Sparse eigval)")

Code Output

Expected Behaviour

I would expect dense and sparse to return the same results.

Your Environment

QuantumToolbox.jl: Quantum Toolbox in Julia
≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡
Copyright © QuTiP team 2022 and later.
Current admin team:
    Alberto Mercurio and Yi-Te Huang

Package information:
====================================
Julia              Ver. 1.11.4
QuantumToolbox     Ver. 0.33.0
SciMLOperators     Ver. 0.4.0
LinearSolve        Ver. 3.23.0
OrdinaryDiffEqCore Ver. 1.26.2

System information:
====================================
OS       : macOS (x86_64-apple-darwin24.0.0)
CPU      : 16 × Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz
Memory   : 16.0 GB
WORD_SIZE: 64
LIBM     : libopenlibm
LLVM     : libLLVM-16.0.6 (ORCJIT, skylake)
BLAS     : libopenblas64_.0.3.27.dylib (ilp64)
Threads  : 12 (on 16 virtual cores)

Additional Context

No response

[albertomercurio]

Hi @eunjongkim.

The sparse solver is efficient to find the largest magnitude eigenvalues. To find the smallest one, one trick is performing the shift-inverse method. This is automatically implemented if you set the sigma keyword argument. This should be a value smaller than the supposed ground state eigenvalue. For your transmon case I think the ground state is around -Ej, so set sigma below this value.

Also note that:

• the which argument doesn't currently work for the sparse method, so we should keep this issue open until this is fixed
• You don't need it because it is very inefficient. Use sigma instead.

[ytdHuang]

Also note that, the function eigsolve in our package doesn't call KrylovKit.eigsolve.

That's why the keyword argument which = :SR does not work.