add keyword argument assume_hermitian

Author: ytdHuang

CHANGELOG.md

@@ -8,7 +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)
 
 - Change default solver detection in `eigensolve` when using `sigma` keyword argument (shift-inverse algorithm). If the operator is a `SparseMatrixCSC`, the default solver is `UMFPACKFactorization`, otherwise it is automatically chosen by LinearSolve.jl, depending on the type of the operator. ([#580]) 
-- Generalize the definition of `liouvillian`. It no longer expects the Hamiltonian to be Hermitian. ([#581])
+- Add keyword argument `assume_hermitian` to `liouvillian`. This allows users to disable the assumption that the Hamiltonian is Hermitian. ([#581])
 
 ## [v0.38.1]
 Release date: 2025-10-27

src/QuantumToolbox.jl

@@ -134,13 +134,4 @@ include("visualization.jl")
 # deprecated functions
 include("deprecated.jl")
 
-# module initialization
-function __init__()
-    # this will be called immediately after the module is loaded (e.g., by import or using) at runtime for the first time.
-    # useful for announcement
-    @warn "The definition of `liouvillian` L for a given Hamiltonian H is changed to L[⋅] = -i( H[⋅] - [⋅]H' ), with ' representing complex conjugation. " *
-          "QuantumToolbox.jl no longer expects H to be Hermitian."
-    return nothing
-end
-
 end

src/qobj/superoperators.jl

@@ -39,18 +39,38 @@ function _spost(B::AbstractSciMLOperator, Id::AbstractMatrix)
     return kron(B_T, Id)
 end
 
-## intrinsic liouvillian 
-function _liouvillian(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
+## intrinsic liouvillian
+_liouvillian(
+    H::MT,
+    Id::AbstractMatrix,
+    assume_hermitian::Val{true},
+) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} = _liouvillian_Herm(H, Id)
+_liouvillian(
+    H::MT,
+    Id::AbstractMatrix,
+    assume_hermitian::Val{false},
+) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} = _liouvillian_NonHerm(H, Id)
+
+## intrinsic liouvillian (assume H is Hermitian)
+_liouvillian_Herm(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} =
+    -1im * (_spre(H, Id) - _spost(H, Id))
+_liouvillian_Herm(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian_Herm(H.A, Id))
+_liouvillian_Herm(H::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(H.λ, _liouvillian_Herm(H.L, Id))
+_liouvillian_Herm(H::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _liouvillian_Herm(op, Id), +, H.ops) # use mapreduce to sum (+) all ops
+
+## intrinsic liouvillian (assume H is Non-Hermitian)
+function _liouvillian_NonHerm(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
     CFType = _complex_float_type(H)
     return CFType(-1.0im) * _spre(H, Id) + CFType(1.0im) * _spost(H', Id)
 end
-_liouvillian(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian(H.A, Id))
-function _liouvillian(H::ScaledOperator, Id::AbstractMatrix)
+_liouvillian_NonHerm(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian_NonHerm(H.A, Id))
+function _liouvillian_NonHerm(H::ScaledOperator, Id::AbstractMatrix)
     CFType = _complex_float_type(H)
     return CFType(-1.0im) * ScaledOperator(H.λ, _spre(H.L, Id)) +
            CFType(1.0im) * ScaledOperator(conj(H.λ), _spost(H.L', Id))
 end
-_liouvillian(H::AddedOperator, Id::AbstractMatrix) = AddedOperator(map(op -> _liouvillian(op, Id), H.ops))
+_liouvillian_NonHerm(H::AddedOperator, Id::AbstractMatrix) =
+    mapreduce(op -> _liouvillian_NonHerm(op, Id), +, H.ops) # use mapreduce to sum (+) all ops
 
 # intrinsic lindblad_dissipator
 function _lindblad_dissipator(O::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
@@ -145,12 +165,25 @@ lindblad_dissipator(O::AbstractQuantumObject{Operator}, Id_cache = I(size(O, 1))
 lindblad_dissipator(O::AbstractQuantumObject{SuperOperator}, Id_cache = nothing) = O
 
 @doc raw"""
-    liouvillian(H::AbstractQuantumObject, c_ops::Union{Nothing,AbstractVector,Tuple}=nothing, Id_cache=I(prod(H.dimensions)))
+    liouvillian(
+        H::AbstractQuantumObject,
+        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing,
+        Id_cache=I(prod(H.dimensions));
+        assume_hermitian::Union{Bool,Val} = Val(true),
+    )
 
-Construct the Liouvillian [`SuperOperator`](@ref) for a system Hamiltonian ``\hat{H}`` and a set of collapse operators ``\{\hat{C}_n\}_n``:
+Construct the Liouvillian [`SuperOperator`](@ref) for a system Hamiltonian ``\hat{H}`` and a set of collapse operators ``\{\hat{C}_n\}_n``.
+
+By default, when the Hamiltonian `H` is assumed to be Hermitian [`assume_hermitian = Val(true)` or `true`], the Liouvillian [`SuperOperator`](@ref) is defined as :
+
+```math
+\mathcal{L} [\cdot] = -i\left(\hat{H}[\cdot] - [\cdot]\hat{H}\right) + \sum_n \mathcal{D}(\hat{C}_n) [\cdot],
+```
+
+otherwise [`assume_hermitian = Val(false)` or `false`],
 
 ```math
-\mathcal{L} [\cdot] = -i\left(\hat{H}[\cdot] - [\cdot]\hat{H}^\dagger\right) + \sum_n \mathcal{D}(\hat{C}_n) [\cdot]
+\mathcal{L} [\cdot] = -i\left(\hat{H}[\cdot] - [\cdot]\hat{H}^\dagger\right) + \sum_n \mathcal{D}(\hat{C}_n) [\cdot],
 ```
 
 where 
@@ -162,27 +195,42 @@ where
 The optional argument `Id_cache` can be used to pass a precomputed identity matrix. This can be useful when the same function is applied multiple times with a known Hilbert space dimension.
 
 See also [`spre`](@ref), [`spost`](@ref), and [`lindblad_dissipator`](@ref).
+
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `assume_hermitian = Val(true)` instead of `assume_hermitian = true`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 function liouvillian(
     H::AbstractQuantumObject{OpType},
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    Id_cache::Diagonal = I(prod(H.dimensions)),
+    Id_cache::Diagonal = I(prod(H.dimensions));
+    assume_hermitian::Union{Bool,Val} = Val(true),
 ) where {OpType<:Union{Operator,SuperOperator}}
-    L = liouvillian(H, Id_cache)
+    L = liouvillian(H, Id_cache; assume_hermitian = assume_hermitian)
     if !(c_ops isa Nothing)
         L += _sum_lindblad_dissipators(c_ops, Id_cache)
     end
     return L
 end
 
-liouvillian(H::Nothing, c_ops::Union{AbstractVector,Tuple}, Id_cache::Diagonal = I(prod(c_ops[1].dims))) =
-    _sum_lindblad_dissipators(c_ops, Id_cache)
-
-liouvillian(H::Nothing, c_ops::Nothing) = 0
-
-liouvillian(H::AbstractQuantumObject{Operator}, Id_cache::Diagonal = I(prod(H.dimensions))) =
-    get_typename_wrapper(H)(_liouvillian(H.data, Id_cache), SuperOperator(), H.dimensions)
-
-liouvillian(H::AbstractQuantumObject{SuperOperator}, Id_cache::Diagonal) = H
+liouvillian(
+    H::Nothing,
+    c_ops::Union{AbstractVector,Tuple},
+    Id_cache::Diagonal = I(prod(c_ops[1].dims));
+    assume_hermitian::Union{Bool,Val} = Val(true),
+) = _sum_lindblad_dissipators(c_ops, Id_cache)
+
+liouvillian(H::Nothing, c_ops::Nothing; assume_hermitian::Union{Bool,Val} = Val(true)) = 0
+
+liouvillian(
+    H::AbstractQuantumObject{Operator},
+    Id_cache::Diagonal = I(prod(H.dimensions));
+    assume_hermitian::Union{Bool,Val} = Val(true),
+) = get_typename_wrapper(H)(_liouvillian(H.data, Id_cache, makeVal(assume_hermitian)), SuperOperator(), H.dimensions)
+
+liouvillian(
+    H::AbstractQuantumObject{SuperOperator},
+    Id_cache::Diagonal;
+    assume_hermitian::Union{Bool,Val} = Val(true),
+) = H
 
 _sum_lindblad_dissipators(c_ops, Id_cache::Diagonal) = sum(op -> lindblad_dissipator(op, Id_cache), c_ops)

src/settings.jl

@@ -15,7 +15,7 @@ function Base.show(io::IO, s::Settings)
 end
 
 @doc raw"""
-    QuantumToolbox.settings 
+    QuantumToolbox.settings
 
 Contains all the default global settings of QuantumToolbox.jl.
 

test/core-test/quantum_objects_evo.jl

@@ -182,15 +182,15 @@
     @testset "Time Dependent Operators and SuperOperators" begin
         N = 10
         a = destroy(N)
-        coef1(p, t) = exp(-1im * p.ω1 * t)
-        coef2(p, t) = sin(p.ω2 * t)
-        coef3(p, t) = sin(p.ω3 * t)
+        coef1(p, t) = exp(-p.γ * t)
+        coef2(p, t) = cos(p.ω1 * t)
+        coef3(p, t) = sin(p.ω2 * t)
         t = rand()
-        p = (ω1 = rand(), ω2 = rand(), ω3 = rand())
+        p = (γ = rand(), ω1 = rand(), ω2 = rand())
 
         # Operator
-        H_td = QobjEvo(((a, coef1), a' * a, (a', coef2)))
-        H_ti = coef1(p, t) * a + a' * a + coef2(p, t) * a'
+        H_td = QobjEvo(((a + a', coef1), a' * a, (-1im * (a - a'), coef2))) # a Hermitian Hamiltonian
+        H_ti = coef1(p, t) * (a + a') + a' * a - 1im * coef2(p, t) * (a - a')
         ψ = rand_ket(N)
         @test H_td(p, t) ≈ H_ti
         @test iscached(H_td) == true
@@ -225,6 +225,13 @@
         @test isconstant(L_td) == false
         @test issuper(L_td) == true
 
+        # test number of lazy operators and assume_hermitian = Val(false)
+        L_td_assume_herm = liouvillian(H_td)
+        L_td_assume_not_herm = liouvillian(H_td, assume_hermitian = Val(false))
+        @test length(L_td_assume_herm.data.ops) == 3 # 1-time-indep. + 2-time-dep.
+        @test length(L_td_assume_not_herm.data.ops) == 5 # 1-time-indep. + 2 x 2-time-dep.
+        @test L_td_assume_herm(p, t) ≈ L_td_assume_not_herm(p, t) # check the matrix since H_td is itself Hermitian
+
         coef_wrong1(t) = nothing
         coef_wrong2(p, t::ComplexF64) = nothing
         @test_logs (:warn,) (:warn,) liouvillian(H_td * H_td) # warnings from lazy tensor
@@ -245,6 +252,7 @@
             @inferred liouvillian(H_td, c_ops2)
             @inferred liouvillian(H_td2, c_ops1)
             @inferred liouvillian(H_td2, c_ops2)
+            @inferred liouvillian(H_td2, c_ops2, assume_hermitian = Val(false))
         end
     end
 end

test/core-test/time_evolution.jl

@@ -1046,6 +1046,11 @@ end
     @test sol_me.expect ≈ sol_me_td2.expect atol = 1e-6 * length(tlist)
     @test sol_mc.expect ≈ sol_mc_td2.expect atol = 1e-2 * length(tlist)
     # @test sol_sse.expect ≈ sol_sse_td2.expect atol = 1e-2 * length(tlist)
+
+    # test assume_hermitian = Val(false)
+    L_td_non_herm = liouvillian(H_td2, c_ops, assume_hermitian = Val(false))
+    sol_me_td3 = mesolve(L_td_non_herm, ψ0, tlist, e_ops = e_ops, progress_bar = Val(false), params = p)
+    @test sol_me.expect ≈ sol_me_td3.expect atol = 1e-6 * length(tlist)
 end
 
 @testitem "mcsolve, ssesolve and smesolve reproducibility" setup=[TESetup] begin

test/ext-test/cpu/autodiff/autodiff.jl

@@ -53,6 +53,7 @@ coef_γ(p, t) = sqrt(p[3])
 H = QobjEvo(a' * a, coef_Δ) + QobjEvo(a + a', coef_F)
 c_ops = [QobjEvo(a, coef_γ)]
 const L = liouvillian(H, c_ops)
+const L_assume_non_herm = liouvillian(H, c_ops, assume_hermitian = Val(false))
 
 function my_f_mesolve(p)
     sol = mesolve(
@@ -67,6 +68,19 @@ function my_f_mesolve(p)
     return real(expect(a' * a, sol.states[end]))
 end
 
+function my_f_mesolve_assume_non_herm(p)
+    sol = mesolve(
+        L_assume_non_herm,
+        ψ0_mesolve,
+        tlist_mesolve,
+        progress_bar = Val(false),
+        params = p,
+        sensealg = BacksolveAdjoint(autojacvec = EnzymeVJP()),
+    )
+
+    return real(expect(a' * a, sol.states[end]))
+end
+
 # Analytical solution
 n_ss(Δ, F, γ) = abs2(F / (Δ + 1im * γ / 2))
 
@@ -113,6 +127,7 @@ n_ss(Δ, F, γ) = abs2(F / (Δ + 1im * γ / 2))
 
         my_f_mesolve_direct(params)
         my_f_mesolve(params)
+        my_f_mesolve_assume_non_herm(params)
 
         grad_exact = Zygote.gradient((p) -> n_ss(p[1], p[2], p[3]), params)[1]
 
@@ -122,20 +137,31 @@ n_ss(Δ, F, γ) = abs2(F / (Δ + 1im * γ / 2))
         end
 
         @testset "Zygote.jl" begin
-            grad_qt = Zygote.gradient(my_f_mesolve, params)[1]
-            @test grad_qt ≈ grad_exact atol=1e-6
+            grad_qt1 = Zygote.gradient(my_f_mesolve, params)[1]
+            grad_qt2 = Zygote.gradient(my_f_mesolve_assume_non_herm, params)[1]
+            @test grad_qt1 ≈ grad_exact atol=1e-6
+            @test grad_qt2 ≈ grad_exact atol=1e-6
         end
 
         @testset "Enzyme.jl" begin
-            dparams = Enzyme.make_zero(params)
+            dparams1 = Enzyme.make_zero(params)
             Enzyme.autodiff(
                 Enzyme.set_runtime_activity(Enzyme.Reverse),
                 my_f_mesolve,
                 Active,
-                Duplicated(params, dparams),
+                Duplicated(params, dparams1),
             )[1]
 
-            @test dparams ≈ grad_exact atol=1e-6
+            dparams2 = Enzyme.make_zero(params)
+            Enzyme.autodiff(
+                Enzyme.set_runtime_activity(Enzyme.Reverse),
+                my_f_mesolve_assume_non_herm,
+                Active,
+                Duplicated(params, dparams2),
+            )[1]
+
+            @test dparams1 ≈ grad_exact atol=1e-6
+            @test dparams2 ≈ grad_exact atol=1e-6
         end
     end
 end