Introduce some entropy related functions

CHANGELOG.md

@@ -11,6 +11,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Change save callbacks from `PresetTimeCallback` to `FunctionCallingCallback`. ([#410])
 - Align `eigenstates` and `eigenenergies` to QuTiP. ([#411])
 - Introduce `vector_to_operator` and `operator_to_vector`. ([#413])
+- Introduce some entropy related functions. ([#416])
+  - `entropy_linear`
+  - `entropy_mutual`
+  - `entropy_conditional`
+  - `entropy_relative`
 
 ## [v0.27.0]
 Release date: 2025-02-14
@@ -149,3 +154,4 @@ Release date: 2024-11-13
 [#410]: /~https://github.com/qutip/QuantumToolbox.jl/issues/410
 [#411]: /~https://github.com/qutip/QuantumToolbox.jl/issues/411
 [#413]: /~https://github.com/qutip/QuantumToolbox.jl/issues/413
+[#416]: /~https://github.com/qutip/QuantumToolbox.jl/issues/416

docs/src/resources/api.md

@@ -250,10 +250,14 @@ ExponentialSeries
 PseudoInverse
 ```
 
-## [Metrics](@id doc-API:Metrics)
+## [Entropy and Metrics](@id doc-API:Entropy-and-Metrics)
 
 ```@docs
 entropy_vn
+entropy_relative
+entropy_linear
+entropy_mutual
+entropy_conditional
 entanglement
 tracedist
 fidelity

src/QuantumToolbox.jl

@@ -117,6 +117,7 @@ include("spectrum.jl")
 include("wigner.jl")
 include("spin_lattice.jl")
 include("arnoldi.jl")
+include("entropy.jl")
 include("metrics.jl")
 include("negativity.jl")
 include("steadystate.jl")

src/entropy.jl

@@ -0,0 +1,211 @@
+#=
+Entropy related functions.
+=#
+
+export entropy_vn, entropy_relative, entropy_linear, entropy_mutual, entropy_conditional
+export entanglement
+
+@doc raw"""
+    entropy_vn(ρ::QuantumObject; base::Int=0, tol::Real=1e-15)
+
+Calculates the [Von Neumann entropy](https://en.wikipedia.org/wiki/Von_Neumann_entropy) ``S = - \textrm{Tr} \left[ \hat{\rho} \log \left( \hat{\rho} \right) \right]``, where ``\hat{\rho}`` is the density matrix of the system.
+
+# Notes
+
+- `ρ` is the quantum state, can be either a [`Ket`](@ref) or an [`Operator`](@ref).
+- `base` specifies the base of the logarithm to use, and when using the default value `0`, the natural logarithm is used.
+- `tol` describes the absolute tolerance for detecting the zero-valued eigenvalues of the density matrix ``\hat{\rho}``.
+
+# Examples
+
+Pure state:
+```jldoctest
+julia> ψ = fock(2,0)
+
+Quantum Object:   type=Ket   dims=[2]   size=(2,)
+2-element Vector{ComplexF64}:
+ 1.0 + 0.0im
+ 0.0 + 0.0im
+
+julia> entropy_vn(ψ, base=2)
+-0.0
+```
+
+Mixed state:
+```jldoctest
+julia> ρ = maximally_mixed_dm(2)
+
+Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
+2×2 Diagonal{ComplexF64, Vector{ComplexF64}}:
+ 0.5-0.0im      ⋅    
+     ⋅      0.5-0.0im
+
+julia> entropy_vn(ρ, base=2)
+1.0
+```
+"""
+function entropy_vn(
+    ρ::QuantumObject{ObjType};
+    base::Int = 0,
+    tol::Real = 1e-15,
+) where {ObjType<:Union{KetQuantumObject,OperatorQuantumObject}}
+    T = eltype(ρ)
+    vals = eigenenergies(ket2dm(ρ))
+    indexes = findall(x -> abs(x) > tol, vals)
+    length(indexes) == 0 && return zero(real(T))
+    nzvals = vals[indexes]
+    logvals = base != 0 ? log.(base, Complex.(nzvals)) : log.(Complex.(nzvals))
+    return -real(mapreduce(*, +, nzvals, logvals))
+end
+
+@doc raw"""
+    entropy_relative(ρ::QuantumObject, σ::QuantumObject; base::Int=0, tol::Real=1e-15)
+
+Calculates the [quantum relative entropy](https://en.wikipedia.org/wiki/Quantum_relative_entropy) of ``\hat{\rho}`` with respect to ``\hat{\sigma}``: ``D(\hat{\rho}||\hat{\sigma}) = \textrm{Tr} \left[ \hat{\rho} \log \left( \hat{\rho} \right) \right] - \textrm{Tr} \left[ \hat{\rho} \log \left( \hat{\sigma} \right) \right]``.
+
+# Notes
+
+- `ρ` is a quantum state, can be either a [`Ket`](@ref) or an [`Operator`](@ref).
+- `σ` is a quantum state, can be either a [`Ket`](@ref) or an [`Operator`](@ref).
+- `base` specifies the base of the logarithm to use, and when using the default value `0`, the natural logarithm is used.
+- `tol` describes the absolute tolerance for detecting the zero-valued eigenvalues of the density matrix ``\hat{\rho}``.
+
+# References
+- [Nielsen-Chuang2011; section 11.3.1, page 511](@citet)
+"""
+function entropy_relative(
+    ρ::QuantumObject{ObjType1},
+    σ::QuantumObject{ObjType2};
+    base::Int = 0,
+    tol::Real = 1e-15,
+) where {
+    ObjType1<:Union{KetQuantumObject,OperatorQuantumObject},
+    ObjType2<:Union{KetQuantumObject,OperatorQuantumObject},
+}
+    check_dimensions(ρ, σ)
+
+    # the logic of this code follows the detail given in the reference of the docstring
+    # consider the eigen decompositions:
+    #   ρ = Σ_i p_i |i⟩⟨i|
+    #   σ = Σ_j q_j |j⟩⟨j|
+    ρ_result = eigenstates(ket2dm(ρ))
+    σ_result = eigenstates(ket2dm(σ))
+
+    # make sure all p_i and q_j are real
+    any(p_i -> imag(p_i) >= tol, ρ_result.values) && error("Input `ρ` has non-real eigenvalues.")
+    any(q_j -> imag(q_j) >= tol, σ_result.values) && error("Input `σ` has non-real eigenvalues.")
+    p = real(ρ_result.values)
+    q = real(σ_result.values)
+    Uρ = ρ_result.vectors
+    Uσ = σ_result.vectors
+
+    # create P_ij matrix (all elements should be real)
+    P = abs2.(Uρ' * Uσ) # this equals to ⟨i|j⟩⟨j|i⟩
+
+    # return +∞ if kernel of σ overlaps with support of ρ, i.e., supp(p) ⊆ supp(q)
+    # That is, if σ is not full rank, S(ρ||σ) = +∞
+    # note that, one special case is that S(ρ||σ) = 0 (if ρ == σ)
+    ((transpose(p .>= tol) * (P .>= tol) * (q .< tol)) == 0) || return Inf
+
+    # Avoid -∞ from log(0), these terms will be multiplied by zero later anyway
+    replace!(q_j -> abs(q_j) < tol ? 1 : q_j, q)
+    p_vals = filter(p_i -> abs(p_i) >= tol, p)
+
+    if base == 0
+        log_p = log.(p_vals)
+        log_q = log.(q)
+    else
+        log_p = log.(base, p_vals)
+        log_q = log.(base, q)
+    end
+
+    # the relative entropy is guaranteed to be ≥ 0
+    # so we calculate the value to 0 to avoid small violations of the lower bound.
+    return max(0.0, dot(p_vals, log_p) - dot(p, P, log_q))
+end
+
+@doc raw"""
+    entropy_linear(ρ::QuantumObject)
+
+Calculates the quantum linear entropy ``S_L = 1 - \textrm{Tr} \left[ \hat{\rho}^2 \right]``, where ``\hat{\rho}`` is the density matrix of the system.
+
+Note that `ρ` can be either a [`Ket`](@ref) or an [`Operator`](@ref).
+"""
+entropy_linear(ρ::QuantumObject{ObjType}) where {ObjType<:Union{KetQuantumObject,OperatorQuantumObject}} =
+    1.0 - purity(ρ) # use 1.0 to make sure it always return value in Float-type
+
+@doc raw"""
+    entropy_mutual(ρAB::QuantumObject, selA, selB; kwargs...)
+
+Calculates the [quantum mutual information](https://en.wikipedia.org/wiki/Quantum_mutual_information) ``I(A:B) = S(\hat{\rho}_A) + S(\hat{\rho}_B) - S(\hat{\rho}_{AB})`` between subsystems ``A`` and ``B``.
+
+Here, ``S`` is the [Von Neumann entropy](https://en.wikipedia.org/wiki/Von_Neumann_entropy), ``\hat{\rho}_{AB}`` is the density matrix of the entire system, ``\hat{\rho}_A = \textrm{Tr}_B \left[ \hat{\rho}_{AB} \right]``, ``\hat{\rho}_B = \textrm{Tr}_A \left[ \hat{\rho}_{AB} \right]``.
+
+# Notes
+
+- `ρAB` can be either a [`Ket`](@ref) or an [`Operator`](@ref).
+- `selA` specifies the indices of the sub-system `A` in `ρAB.dimensions`. See also [`ptrace`](@ref).
+- `selB` specifies the indices of the sub-system `B` in `ρAB.dimensions`. See also [`ptrace`](@ref).
+- `kwargs` are the keyword arguments for calculating Von Neumann entropy. See also [`entropy_vn`](@ref).
+"""
+function entropy_mutual(
+    ρAB::QuantumObject{ObjType,<:AbstractDimensions{N}},
+    selA::Union{Int,AbstractVector{Int},Tuple},
+    selB::Union{Int,AbstractVector{Int},Tuple};
+    kwargs...,
+) where {ObjType<:Union{KetQuantumObject,OperatorQuantumObject},N}
+    # check if selA and selB matches the dimensions of ρAB
+    sel_A_B = (selA..., selB...)
+    (length(sel_A_B) != N) && throw(
+        ArgumentError(
+            "The indices in `selA = $(selA)` and `selB = $(selB)` does not match the given QuantumObject which has $N sub-systems",
+        ),
+    )
+    allunique(sel_A_B) || throw(ArgumentError("Duplicate selection indices in `selA = $(selA)` and `selB = $(selB)`"))
+
+    ρA = ptrace(ρAB, selA)
+    ρB = ptrace(ρAB, selB)
+    return entropy_vn(ρA; kwargs...) + entropy_vn(ρB; kwargs...) - entropy_vn(ρAB; kwargs...)
+end
+
+@doc raw"""
+    entropy_conditional(ρAB::QuantumObject, selB; kwargs...)
+
+Calculates the [conditional quantum entropy](https://en.wikipedia.org/wiki/Conditional_quantum_entropy) with respect to sub-system ``B``: ``S(A|B) = S(\hat{\rho}_{AB}) - S(\hat{\rho}_{B})``.
+
+Here, ``S`` is the [Von Neumann entropy](https://en.wikipedia.org/wiki/Von_Neumann_entropy), ``\hat{\rho}_{AB}`` is the density matrix of the entire system, and ``\hat{\rho}_B = \textrm{Tr}_A \left[ \hat{\rho}_{AB} \right]``.
+
+# Notes
+
+- `ρAB` can be either a [`Ket`](@ref) or an [`Operator`](@ref).
+- `selB` specifies the indices of the sub-system `B` in `ρAB.dimensions`. See also [`ptrace`](@ref).
+- `kwargs` are the keyword arguments for calculating Von Neumann entropy. See also [`entropy_vn`](@ref).
+"""
+entropy_conditional(
+    ρAB::QuantumObject{ObjType,<:AbstractDimensions{N}},
+    selB::Union{Int,AbstractVector{Int},Tuple};
+    kwargs...,
+) where {ObjType<:Union{KetQuantumObject,OperatorQuantumObject},N} =
+    entropy_vn(ρAB; kwargs...) - entropy_vn(ptrace(ρAB, selB); kwargs...)
+
+@doc raw"""
+    entanglement(ρ::QuantumObject, sel; kwargs...)
+
+Calculates the [entanglement entropy](https://en.wikipedia.org/wiki/Entropy_of_entanglement) by doing the partial trace of `ρ`, selecting only the dimensions with the indices contained in the `sel` vector, and then use the Von Neumann entropy [`entropy_vn`](@ref).
+
+# Notes
+
+- `ρ` can be either a [`Ket`](@ref) or an [`Operator`](@ref).
+- `sel` specifies the indices of the remaining sub-system. See also [`ptrace`](@ref).
+- `kwargs` are the keyword arguments for calculating Von Neumann entropy. See also [`entropy_vn`](@ref).
+"""
+function entanglement(
+    ρ::QuantumObject{OpType},
+    sel::Union{Int,AbstractVector{Int},Tuple},
+    kwargs...,
+) where {OpType<:Union{KetQuantumObject,OperatorQuantumObject}}
+    _ρ = normalize(ρ)
+    ρ_tr = ptrace(_ρ, sel)
+    val = entropy_vn(ρ_tr; kwargs...)
+    return (val > 0) * val
+end

src/metrics.jl

@@ -2,81 +2,7 @@
 Functions for calculating metrics (distance measures) between states and operators.
 =#
 
-export entropy_vn, entanglement, tracedist, fidelity
-
-@doc raw"""
-    entropy_vn(ρ::QuantumObject; base::Int=0, tol::Real=1e-15)
-
-Calculates the [Von Neumann entropy](https://en.wikipedia.org/wiki/Von_Neumann_entropy)
-``S = - \textrm{Tr} \left[ \hat{\rho} \log \left( \hat{\rho} \right) \right]`` where ``\hat{\rho}``
-is the density matrix of the system.
-
-The `base` parameter specifies the base of the logarithm to use, and when using the default value 0,
-the natural logarithm is used. The `tol` parameter
-describes the absolute tolerance for detecting the zero-valued eigenvalues of the density
-matrix ``\hat{\rho}``.
-
-# Examples
-
-Pure state:
-```jldoctest
-julia> ψ = fock(2,0)
-
-Quantum Object:   type=Ket   dims=[2]   size=(2,)
-2-element Vector{ComplexF64}:
- 1.0 + 0.0im
- 0.0 + 0.0im
-
-julia> ρ = ket2dm(ψ)
-
-Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
-2×2 Matrix{ComplexF64}:
- 1.0+0.0im  0.0+0.0im
- 0.0+0.0im  0.0+0.0im
-
-julia> entropy_vn(ρ, base=2)
--0.0
-```
-
-Mixed state:
-```jldoctest
-julia> ρ = maximally_mixed_dm(2)
-
-Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
-2×2 Diagonal{ComplexF64, Vector{ComplexF64}}:
- 0.5-0.0im      ⋅    
-     ⋅      0.5-0.0im
-
-julia> entropy_vn(ρ, base=2)
-1.0
-```
-"""
-function entropy_vn(ρ::QuantumObject{OperatorQuantumObject}; base::Int = 0, tol::Real = 1e-15)
-    T = eltype(ρ)
-    vals = eigenenergies(ρ)
-    indexes = findall(x -> abs(x) > tol, vals)
-    length(indexes) == 0 && return zero(real(T))
-    nzvals = vals[indexes]
-    logvals = base != 0 ? log.(base, Complex.(nzvals)) : log.(Complex.(nzvals))
-    return -real(mapreduce(*, +, nzvals, logvals))
-end
-
-"""
-    entanglement(QO::QuantumObject, sel::Union{Int,AbstractVector{Int},Tuple})
-
-Calculates the entanglement by doing the partial trace of `QO`, selecting only the dimensions
-with the indices contained in the `sel` vector, and then using the Von Neumann entropy [`entropy_vn`](@ref).
-"""
-function entanglement(
-    QO::QuantumObject{OpType},
-    sel::Union{AbstractVector{Int},Tuple},
-) where {OpType<:Union{BraQuantumObject,KetQuantumObject,OperatorQuantumObject}}
-    ψ = normalize(QO)
-    ρ_tr = ptrace(ψ, sel)
-    entropy = entropy_vn(ρ_tr)
-    return (entropy > 0) * entropy
-end
-entanglement(QO::QuantumObject, sel::Int) = entanglement(QO, (sel,))
+export tracedist, fidelity
 
 @doc raw"""
     tracedist(ρ::QuantumObject, σ::QuantumObject)

src/qobj/arithmetic_and_attributes.jl

@@ -568,16 +568,15 @@ Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
 function ptrace(QO::QuantumObject{KetQuantumObject}, sel::Union{AbstractVector{Int},Tuple})
     _non_static_array_warning("sel", sel)
 
-    n_s = length(sel)
-    if n_s == 0 # return full trace for empty sel
+    if length(sel) == 0 # return full trace for empty sel
         return tr(ket2dm(QO))
     else
         n_d = length(QO.dimensions)
 
         (any(>(n_d), sel) || any(<(1), sel)) && throw(
             ArgumentError("Invalid indices in `sel`: $(sel), the given QuantumObject only have $(n_d) sub-systems"),
         )
-        (n_s != length(unique(sel))) && throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
+        allunique(sel) || throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
         (n_d == 1) && return ket2dm(QO) # ptrace should always return Operator
     end
 
@@ -596,16 +595,15 @@ function ptrace(QO::QuantumObject{OperatorQuantumObject}, sel::Union{AbstractVec
 
     _non_static_array_warning("sel", sel)
 
-    n_s = length(sel)
-    if n_s == 0 # return full trace for empty sel
+    if length(sel) == 0 # return full trace for empty sel
         return tr(QO)
     else
         n_d = length(QO.dimensions)
 
         (any(>(n_d), sel) || any(<(1), sel)) && throw(
             ArgumentError("Invalid indices in `sel`: $(sel), the given QuantumObject only have $(n_d) sub-systems"),
         )
-        (n_s != length(unique(sel))) && throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
+        allunique(sel) || throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
         (n_d == 1) && return QO
     end
 

src/qobj/functions.jl

@@ -125,7 +125,7 @@ to_dense(::Type{T}, A::AbstractSparseArray) where {T<:Number} = Array{T}(A)
 to_dense(::Type{T1}, A::AbstractArray{T2}) where {T1<:Number,T2<:Number} = Array{T1}(A)
 to_dense(::Type{T}, A::AbstractArray{T}) where {T<:Number} = A
 
-function to_dense(::Type{M}) where {M<:SparseMatrixCSC}
+function to_dense(::Type{M}) where {M<:Union{Diagonal,SparseMatrixCSC}}
     T = M
     par = T.parameters
     npar = length(par)