add tests for entropy functions

Author: ytdHuang

src/entropy.jl

@@ -100,10 +100,12 @@ function entropy_relative(
     Uσ = σ_result.vectors
 
     # create P_ij matrix (all elements should be real)
-    P = abs2(Uρ' * Uσ) # this equals to ⟨i|j⟩⟨j|i⟩
+    P = abs2.(Uρ' * Uσ) # this equals to ⟨i|j⟩⟨j|i⟩
 
-    # return +∞ if kernel of σ overlaps with support of ρ
-    dot((p .>= tol), (P .>= tol), (q .< tol)) && return Inf
+    # 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)
@@ -119,7 +121,7 @@ function entropy_relative(
 
     # 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, dot(p_vals, log_p) - dot(p, P, log_q))
+    return max(0.0, dot(p_vals, log_p) - dot(p, P, log_q))
 end
 
 @doc raw"""
@@ -158,9 +160,11 @@ function entropy_mutual(
     BType<:Union{Int,AbstractVector{Int},Tuple},
 }
     # check if selA and selB matches the dimensions of ρAB
-    sel_A_B = vcat(selA, selB)
-    (length(sel_A_B) != N) && ArgumentError(
-        "The indices in `selA = $(selA)` and `selB = $(selB)` does not match the given QuantumObject which has $N sub-systems",
+    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)`"))
 
@@ -196,11 +200,10 @@ Calculates the entanglement by doing the partial trace of `QO`, selecting only t
 """
 function entanglement(
     QO::QuantumObject{OpType},
-    sel::Union{AbstractVector{Int},Tuple},
+    sel::Union{Int,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,))

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)

test/core-test/entanglement.jl

@@ -0,0 +1,101 @@
+@testset "entropy" begin
+    base = 2
+    λ = rand()
+    ψ = rand_ket(10)
+    ρ1 = rand_dm(10)
+    ρ2 = rand_dm(10)
+    σ1 = rand_dm(10)
+    σ2 = rand_dm(10)
+
+    dims = (2, 3)
+    ρAB = rand_dm((dims..., dims...))
+    selA = (1, 2)
+    selB = (3, 4)
+    ρA = ptrace(ρAB, selA)
+    ρB = ptrace(ρAB, selB)
+    nA = nB = prod(dims)
+    IA = qeye(nA, dims = dims)
+    IB = qeye(nB, dims = dims)
+
+    # quantum relative entropy
+    @test entropy_relative(ρ1, ψ) == Inf
+    @test entropy_relative(ρ1, rand_dm(10, rank = 9)) == Inf
+    @test entropy_relative(ψ, ψ) ≈ 0
+    @test entropy_relative(λ * ρ1 + (1 - λ) * ρ2, λ * σ1 + (1 - λ) * σ2) <=
+          λ * entropy_relative(ρ1, σ1) + (1 - λ) * entropy_relative(ρ2, σ2) # joint convexity
+
+    # relations between different entropies
+    @test entropy_relative(ρA, IA / nA) ≈ log(nA) - entropy_vn(ρA)
+    @test entropy_relative(ρB, IB / nB; base = base) ≈ log(base, nB) - entropy_vn(ρB; base = base)
+    @test entropy_relative(ρAB, tensor(ρA, ρB)) ≈ entropy_mutual(ρAB, selA, selB)
+    @test entropy_relative(ρAB, tensor(ρA, ρB)) ≈ entropy_mutual(ρAB, selA, selB)
+    @test entropy_relative(ρAB, tensor(ρA, IB / nB)) ≈ log(nB) - entropy_conditional(ρAB, selA)
+    @test entropy_linear(ρ1) == 1 - purity(ρ1)
+
+    ρ_wrong = Qobj(rand(ComplexF64, 10, 10))
+    @test_throws ErrorException entropy_relative(ρ1, ρ_wrong)
+    @test_throws ErrorException entropy_relative(ρ_wrong, ρ1)
+    @test_throws ArgumentError entropy_mutual(ρAB, 1, 3)
+    @test_throws ArgumentError entropy_mutual(ρAB, 1, (3, 4))
+    @test_throws ArgumentError entropy_mutual(ρAB, (1, 2), 3)
+    @test_throws ArgumentError entropy_mutual(ρAB, (1, 2), (1, 3))
+
+    @testset "Type Stability (entropy)" begin
+        @inferred entropy_vn(ρ1)
+        @inferred entropy_vn(ρ1, base = base)
+        @inferred entropy_relative(ρ1, ψ)
+        @inferred entropy_relative(ρ1, σ1, base = base)
+        @inferred entropy_linear(ρ1)
+        @inferred entropy_mutual(ρAB, selA, selB)
+        @inferred entropy_conditional(ρAB, selA)
+    end
+end
+
+@testset "Entanglement" begin
+    g = fock(2, 1)
+    e = fock(2, 0)
+    state = normalize(kron(g, e) + kron(e, g))
+    rho = state * state'
+    @test entanglement(state, 1) / log(2) ≈ 1
+    @test entanglement(rho, 1) / log(2) ≈ 1
+
+    @testset "Type Stability (entanglement)" begin
+        @inferred entanglement(state, 1)
+        @inferred entanglement(rho, 1)
+    end
+end
+
+@testset "trace distance" begin
+    ψz0 = basis(2, 0)
+    ψz1 = basis(2, 1)
+    ρz0 = to_sparse(ket2dm(ψz0))
+    ρz1 = to_sparse(ket2dm(ψz1))
+    ψx0 = sqrt(0.5) * (basis(2, 0) + basis(2, 1))
+    @test tracedist(ψz0, ψx0) ≈ sqrt(0.5)
+    @test tracedist(ρz0, ψz1) ≈ 1.0
+    @test tracedist(ψz1, ρz0) ≈ 1.0
+    @test tracedist(ρz0, ρz1) ≈ 1.0
+
+    @testset "Type Inference (trace distance)" begin
+        @inferred tracedist(ψz0, ψx0)
+        @inferred tracedist(ρz0, ψz1)
+        @inferred tracedist(ψz1, ρz0)
+        @inferred tracedist(ρz0, ρz1)
+    end
+end
+
+@testset "fidelity" begin
+    M = sprand(ComplexF64, 5, 5, 0.5)
+    M0 = Qobj(M * M')
+    ψ1 = Qobj(rand(ComplexF64, 5))
+    ψ2 = Qobj(rand(ComplexF64, 5))
+    M1 = ψ1 * ψ1'
+    @test isapprox(fidelity(M0, M1), fidelity(ψ1, M0); atol = 1e-6)
+    @test isapprox(fidelity(ψ1, ψ2), fidelity(ket2dm(ψ1), ket2dm(ψ2)); atol = 1e-6)
+
+    @testset "Type Inference (fidelity)" begin
+        @inferred fidelity(M0, M1)
+        @inferred fidelity(ψ1, M0)
+        @inferred fidelity(ψ1, ψ2)
+    end
+end

test/core-test/entropy_and_metric.jl

@@ -566,54 +566,20 @@
         end
     end
 
-    @testset "trace distance" begin
-        ψz0 = basis(2, 0)
-        ψz1 = basis(2, 1)
-        ρz0 = to_sparse(ket2dm(ψz0))
-        ρz1 = to_sparse(ket2dm(ψz1))
-        ψx0 = sqrt(0.5) * (basis(2, 0) + basis(2, 1))
-        @test tracedist(ψz0, ψx0) ≈ sqrt(0.5)
-        @test tracedist(ρz0, ψz1) ≈ 1.0
-        @test tracedist(ψz1, ρz0) ≈ 1.0
-        @test tracedist(ρz0, ρz1) ≈ 1.0
-
-        @testset "Type Inference (trace distance)" begin
-            @inferred tracedist(ψz0, ψx0)
-            @inferred tracedist(ρz0, ψz1)
-            @inferred tracedist(ψz1, ρz0)
-            @inferred tracedist(ρz0, ρz1)
-        end
-    end
-
-    @testset "sqrt and fidelity" begin
-        M = sprand(ComplexF64, 5, 5, 0.5)
-        M0 = Qobj(M * M')
-        ψ1 = Qobj(rand(ComplexF64, 5))
-        ψ2 = Qobj(rand(ComplexF64, 5))
-        M1 = ψ1 * ψ1'
-        @test sqrtm(M0) ≈ sqrtm(to_dense(M0))
-        @test isapprox(fidelity(M0, M1), fidelity(ψ1, M0); atol = 1e-6)
-        @test isapprox(fidelity(ψ1, ψ2), fidelity(ket2dm(ψ1), ket2dm(ψ2)); atol = 1e-6)
-
-        @testset "Type Inference (sqrt and fidelity)" begin
-            @inferred sqrtm(M0)
-            @inferred fidelity(M0, M1)
-            @inferred fidelity(ψ1, M0)
-            @inferred fidelity(ψ1, ψ2)
-        end
-    end
-
-    @testset "log, exp, sinm, cosm" begin
+    @testset "sqrt, log, exp, sinm, cosm" begin
         M0 = rand(ComplexF64, 4, 4)
         Md = Qobj(M0 * M0')
         Ms = to_sparse(Md)
         e_p = expm(1im * Md)
         e_m = expm(-1im * Md)
+        @test sqrtm(Md) ≈ sqrtm(Ms)
         @test logm(expm(Ms)) ≈ expm(logm(Md))
         @test cosm(Ms) ≈ (e_p + e_m) / 2
         @test sinm(Ms) ≈ (e_p - e_m) / 2im
 
         @testset "Type Inference" begin
+            @inferred sqrtm(Md)
+            @inferred sqrtm(Ms)
             @inferred expm(Md)
             @inferred expm(Ms)
             @inferred logm(Md)

test/core-test/quantum_objects.jl

@@ -12,7 +12,7 @@ core_tests = [
     "dynamical_fock_dimension_mesolve.jl",
     "dynamical-shifted-fock.jl",
     "eigenvalues_and_operators.jl",
-    "entanglement.jl",
+    "entropy_and_metric.jl",
     "generalized_master_equation.jl",
     "low_rank_dynamics.jl",
     "negativity_and_partial_transpose.jl",