From 3612c7b5d1569979fc9294f38c244b7a42d8ab72 Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Mon, 17 Feb 2025 21:26:45 +0100 Subject: [PATCH 1/2] [no ci] Align `eigenstates` and `eigenenergies` to QuTiP --- benchmarks/eigenvalues.jl | 3 +- src/qobj/eigsolve.jl | 43 +++++++++++++-------- src/steadystate.jl | 2 +- test/core-test/eigenvalues_and_operators.jl | 16 ++++---- 4 files changed, 37 insertions(+), 27 deletions(-) diff --git a/benchmarks/eigenvalues.jl b/benchmarks/eigenvalues.jl index 3fc4287a0..3fab5302d 100644 --- a/benchmarks/eigenvalues.jl +++ b/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 diff --git a/src/qobj/eigsolve.jl b/src/qobj/eigsolve.jl index 73db30c7b..f180c0272 100644 --- a/src/qobj/eigsolve.jl +++ b/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) diff --git a/src/steadystate.jl b/src/steadystate.jl index 94cbd3222..e0296954a 100644 --- a/src/steadystate.jl +++ b/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) diff --git a/test/core-test/eigenvalues_and_operators.jl b/test/core-test/eigenvalues_and_operators.jl index 481aedce0..3faebf0dc 100644 --- a/test/core-test/eigenvalues_and_operators.jl +++ b/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 From 5c4dd042ff1194fd0886613a90da14770fe17f3f Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Mon, 17 Feb 2025 21:28:31 +0100 Subject: [PATCH 2/2] Make changelog --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index ce0837a02..253449b48 100644 --- a/CHANGELOG.md +++ b/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