Align eigenstates and eigenenergies to QuTiP (#411)

Author: albertomercurio

CHANGELOG.md

@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
 - Support for single `AbstractQuantumObject` in `sc_ops` for faster specific method in `ssesolve` and `smesolve`. ([#408])
+- Align `eigenstates` and `eigenenergies` to QuTiP. ([#411])
 
 ## [v0.27.0]
 Release date: 2025-02-14
@@ -143,3 +144,4 @@ Release date: 2024-11-13
 [#404]: https://github.com/qutip/QuantumToolbox.jl/issues/404
 [#405]: https://github.com/qutip/QuantumToolbox.jl/issues/405
 [#408]: https://github.com/qutip/QuantumToolbox.jl/issues/408
+[#411]: https://github.com/qutip/QuantumToolbox.jl/issues/411

benchmarks/eigenvalues.jl

@@ -16,7 +16,8 @@ function benchmark_eigenvalues!(SUITE)
     L = liouvillian(H, c_ops)
 
     SUITE["Eigenvalues"]["eigenstates"]["dense"] = @benchmarkable eigenstates($L)
-    SUITE["Eigenvalues"]["eigenstates"]["sparse"] = @benchmarkable eigenstates($L, sparse = true, sigma = 0.01, k = 5)
+    SUITE["Eigenvalues"]["eigenstates"]["sparse"] =
+        @benchmarkable eigenstates($L, sparse = true, sigma = 0.01, eigvals = 5)
 
     return nothing
 end

src/qobj/eigsolve.jl

@@ -257,7 +257,7 @@ end
     eigsolve(A::QuantumObject; 
         v0::Union{Nothing,AbstractVector}=nothing, 
         sigma::Union{Nothing, Real}=nothing,
-        k::Int = 1,
+        eigvals::Int = 1,
         krylovdim::Int = max(20, 2*k+1),
         tol::Real = 1e-8,
         maxiter::Int = 200,
@@ -266,6 +266,17 @@ end
 
 Solve for the eigenvalues and eigenvectors of a matrix `A` using the Arnoldi method.
 
+# Arguments
+- `A::QuantumObject`: the [`QuantumObject`](@ref) to solve eigenvalues and eigenvectors.
+- `v0::Union{Nothing,AbstractVector}`: the initial vector for the Arnoldi method. Default is a random vector.
+- `sigma::Union{Nothing, Real}`: the shift for the eigenvalue problem. Default is `nothing`.
+- `eigvals::Int`: the number of eigenvalues to compute. Default is `1`.
+- `krylovdim::Int`: the dimension of the Krylov subspace. Default is `max(20, 2*k+1)`.
+- `tol::Real`: the tolerance for the Arnoldi method. Default is `1e-8`.
+- `maxiter::Int`: the maximum number of iterations for the Arnoldi method. Default is `200`.
+- `solver::Union{Nothing, SciMLLinearSolveAlgorithm}`: the linear solver algorithm. Default is `nothing`.
+- `kwargs`: Additional keyword arguments passed to the solver.
+
 # Notes
 - For more details about `solver` and extra `kwargs`, please refer to [`LinearSolve.jl`](https://docs.sciml.ai/LinearSolve/stable/)
 
@@ -276,8 +287,8 @@ function eigsolve(
     A::QuantumObject;
     v0::Union{Nothing,AbstractVector} = nothing,
     sigma::Union{Nothing,Real} = nothing,
-    k::Int = 1,
-    krylovdim::Int = max(20, 2 * k + 1),
+    eigvals::Int = 1,
+    krylovdim::Int = max(20, 2 * eigvals + 1),
     tol::Real = 1e-8,
     maxiter::Int = 200,
     solver::Union{Nothing,SciMLLinearSolveAlgorithm} = nothing,
@@ -289,7 +300,7 @@ function eigsolve(
         type = A.type,
         dimensions = A.dimensions,
         sigma = sigma,
-        k = k,
+        eigvals = eigvals,
         krylovdim = krylovdim,
         tol = tol,
         maxiter = maxiter,
@@ -304,8 +315,8 @@ function eigsolve(
     type::Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject} = nothing,
     dimensions = nothing,
     sigma::Union{Nothing,Real} = nothing,
-    k::Int = 1,
-    krylovdim::Int = max(20, 2 * k + 1),
+    eigvals::Int = 1,
+    krylovdim::Int = max(20, 2 * eigvals + 1),
     tol::Real = 1e-8,
     maxiter::Int = 200,
     solver::Union{Nothing,SciMLLinearSolveAlgorithm} = nothing,
@@ -316,7 +327,7 @@ function eigsolve(
     v0 === nothing && (v0 = normalize!(rand(T, size(A, 1))))
 
     if sigma === nothing
-        res = _eigsolve(A, v0, type, dimensions, k, krylovdim, tol = tol, maxiter = maxiter)
+        res = _eigsolve(A, v0, type, dimensions, eigvals, krylovdim, tol = tol, maxiter = maxiter)
         vals = res.values
     else
         Aₛ = A - sigma * I
@@ -334,7 +345,7 @@ function eigsolve(
 
         Amap = EigsolveInverseMap(T, size(A), linsolve)
 
-        res = _eigsolve(Amap, v0, type, dimensions, k, krylovdim, tol = tol, maxiter = maxiter)
+        res = _eigsolve(Amap, v0, type, dimensions, eigvals, krylovdim, tol = tol, maxiter = maxiter)
         vals = @. (1 + sigma * res.values) / res.values
     end
 
@@ -349,7 +360,7 @@ end
         alg::OrdinaryDiffEqAlgorithm = Tsit5(),
         params::NamedTuple = NamedTuple(),
         ρ0::AbstractMatrix = rand_dm(prod(H.dimensions)).data,
-        k::Int = 1,
+        eigvals::Int = 1,
         krylovdim::Int = min(10, size(H, 1)),
         maxiter::Int = 200,
         eigstol::Real = 1e-6,
@@ -365,7 +376,7 @@ Solve the eigenvalue problem for a Liouvillian superoperator `L` using the Arnol
 - `alg`: The differential equation solver algorithm. Default is `Tsit5()`.
 - `params`: A `NamedTuple` containing the parameters of the system.
 - `ρ0`: The initial density matrix. If not specified, a random density matrix is used.
-- `k`: The number of eigenvalues to compute.
+- `eigvals`: The number of eigenvalues to compute.
 - `krylovdim`: The dimension of the Krylov subspace.
 - `maxiter`: The maximum number of iterations for the eigsolver.
 - `eigstol`: The tolerance for the eigsolver.
@@ -388,7 +399,7 @@ function eigsolve_al(
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
     params::NamedTuple = NamedTuple(),
     ρ0::AbstractMatrix = rand_dm(prod(H.dimensions)).data,
-    k::Int = 1,
+    eigvals::Int = 1,
     krylovdim::Int = min(10, size(H, 1)),
     maxiter::Int = 200,
     eigstol::Real = 1e-6,
@@ -406,12 +417,10 @@ function eigsolve_al(
         ).prob
     integrator = init(prob, alg)
 
-    # prog = ProgressUnknown(desc="Applications:", showspeed = true, enabled=progress)
-
     Lmap = ArnoldiLindbladIntegratorMap(eltype(H), size(L_evo), integrator)
 
-    res = _eigsolve(Lmap, mat2vec(ρ0), L_evo.type, L_evo.dimensions, k, krylovdim, maxiter = maxiter, tol = eigstol)
-    # finish!(prog)
+    res =
+        _eigsolve(Lmap, mat2vec(ρ0), L_evo.type, L_evo.dimensions, eigvals, krylovdim, maxiter = maxiter, tol = eigstol)
 
     vals = similar(res.values)
     vecs = similar(res.vectors)
@@ -488,7 +497,7 @@ Calculate the eigenenergies
 # Arguments
 - `A::QuantumObject`: the [`QuantumObject`](@ref) to solve eigenvalues
 - `sparse::Bool`: if `false` call [`eigvals(A::QuantumObject; kwargs...)`](@ref), otherwise call [`eigsolve`](@ref). Default to `false`.
-- `kwargs`: Additional keyword arguments passed to the solver
+- `kwargs`: Additional keyword arguments passed to the solver. If `sparse=true`, the keyword arguments are passed to [`eigsolve`](@ref), otherwise to [`eigen`](@ref).
 
 # Returns
 - `::Vector{<:Number}`: a list of eigenvalues
@@ -513,7 +522,7 @@ Calculate the eigenvalues and corresponding eigenvectors
 # Arguments
 - `A::QuantumObject`: the [`QuantumObject`](@ref) to solve eigenvalues and eigenvectors
 - `sparse::Bool`: if `false` call [`eigen(A::QuantumObject; kwargs...)`](@ref), otherwise call [`eigsolve`](@ref). Default to `false`.
-- `kwargs`: Additional keyword arguments passed to the solver
+- `kwargs`: Additional keyword arguments passed to the solver. If `sparse=true`, the keyword arguments are passed to [`eigsolve`](@ref), otherwise to [`eigen`](@ref).
 
 # Returns
 - `::EigsolveResult`: containing the eigenvalues, the eigenvectors, and some information from the solver. see also [`EigsolveResult`](@ref)

src/steadystate.jl

@@ -141,7 +141,7 @@ end
 function _steadystate(L::QuantumObject{SuperOperatorQuantumObject}, solver::SteadyStateEigenSolver; kwargs...)
     N = prod(L.dimensions)
 
-    kwargs = merge((sigma = 1e-8, k = 1), (; kwargs...))
+    kwargs = merge((sigma = 1e-8, eigvals = 1), (; kwargs...))
 
     ρss_vec = eigsolve(L; kwargs...).vectors[:, 1]
     ρss = reshape(ρss_vec, N, N)

test/core-test/eigenvalues_and_operators.jl

@@ -5,15 +5,15 @@
     resstring = sprint((t, s) -> show(t, "text/plain", s), result)
     valstring = sprint((t, s) -> show(t, "text/plain", s), result.values)
     vecsstring = sprint((t, s) -> show(t, "text/plain", s), result.vectors)
-    λs, ψs, Ts = eigenstates(σx, sparse = true, k = 2)
-    λs1, ψs1, Ts1 = eigenstates(σx, sparse = true, k = 1)
+    λs, ψs, Ts = eigenstates(σx, sparse = true, eigvals = 2)
+    λs1, ψs1, Ts1 = eigenstates(σx, sparse = true, eigvals = 1)
 
     @test all([ψ.type isa KetQuantumObject for ψ in ψd])
     @test typeof(Td) <: AbstractMatrix
     @test typeof(Ts) <: AbstractMatrix
     @test typeof(Ts1) <: AbstractMatrix
     @test all(abs.(eigenenergies(σx, sparse = false)) .≈ abs.(λd))
-    @test all(abs.(eigenenergies(σx, sparse = true, k = 2)) .≈ abs.(λs))
+    @test all(abs.(eigenenergies(σx, sparse = true, eigvals = 2)) .≈ abs.(λs))
     @test resstring ==
           "EigsolveResult:   type=$(Operator)   dims=$(result.dims)\nvalues:\n$(valstring)\nvectors:\n$vecsstring"
 
@@ -33,7 +33,7 @@
 
     vals_d, vecs_d, mat_d = eigenstates(H_d)
     vals_c, vecs_c, mat_c = eigenstates(H_c)
-    vals2, vecs2, mat2 = eigenstates(H_d, sparse = true, sigma = -0.9, k = 10, krylovdim = 30)
+    vals2, vecs2, mat2 = eigenstates(H_d, sparse = true, sigma = -0.9, eigvals = 10, krylovdim = 30)
     sort!(vals_c, by = real)
     sort!(vals2, by = real)
 
@@ -57,9 +57,9 @@
     L = liouvillian(H, c_ops)
 
     # eigen solve for general matrices
-    vals, _, vecs = eigsolve(L.data, sigma = 0.01, k = 10, krylovdim = 50)
+    vals, _, vecs = eigsolve(L.data, sigma = 0.01, eigvals = 10, krylovdim = 50)
     vals2, vecs2 = eigen(to_dense(L.data))
-    vals3, state3, vecs3 = eigsolve_al(L, 1 \ (40 * κ), k = 10, krylovdim = 50)
+    vals3, state3, vecs3 = eigsolve_al(L, 1 \ (40 * κ), eigvals = 10, krylovdim = 50)
     idxs = sortperm(vals2, by = abs)
     vals2 = vals2[idxs][1:10]
     vecs2 = vecs2[:, idxs][:, 1:10]
@@ -70,7 +70,7 @@
     @test isapprox(vec2mat(vecs[:, 1]) * exp(-1im * angle(vecs[1, 1])), vec2mat(vecs3[:, 1]), atol = 1e-5)
 
     # eigen solve for QuantumObject
-    result = eigenstates(L, sparse = true, sigma = 0.01, k = 10, krylovdim = 50)
+    result = eigenstates(L, sparse = true, sigma = 0.01, eigvals = 10, krylovdim = 50)
     vals, vecs = result
     resstring = sprint((t, s) -> show(t, "text/plain", s), result)
     valstring = sprint((t, s) -> show(t, "text/plain", s), result.values)
@@ -112,6 +112,6 @@
         @inferred eigenstates(H, sparse = false)
         @inferred eigenstates(H, sparse = true)
         @inferred eigenstates(L, sparse = true)
-        @inferred eigsolve_al(L, 1 \ (40 * κ), k = 10)
+        @inferred eigsolve_al(L, 1 \ (40 * κ), eigvals = 10)
     end
 end