Merge pull request #102 from ytdHuang/dev/tracedist

Support `svdvals`, generalize `norm`, and introduce `tracedist` for `Qobj`

docs/src/api.md

@@ -32,6 +32,7 @@ size
 eltype
 length
 LinearAlgebra.tr
+LinearAlgebra.svdvals
 LinearAlgebra.norm
 LinearAlgebra.kron
 tensor
@@ -57,6 +58,7 @@ expect
 wigner
 get_coherence
 n_th
+tracedist
 get_data
 mat2vec
 vec2mat

src/general_functions.jl

@@ -1,6 +1,6 @@
 export get_data, get_coherence, expect, ptrace
 export mat2vec, vec2mat
-export entropy_vn, entanglement
+export entropy_vn, entanglement, tracedist
 export gaussian, n_th
 
 export row_major_reshape, tidyup, tidyup!, meshgrid, sparse_to_dense, dense_to_sparse
@@ -328,9 +328,15 @@ function n_th(ω::Real, T::Real)::Float64
     return 1 / (exp(ω / T) - 1)
 end
 
+@doc raw"""
+    tracedist(ρ::QuantumObject, σ::QuantumObject)
 
+Calculates the [trace distance](https://en.wikipedia.org/wiki/Trace_distance) between two [`QuantumObject`](@ref):
+``T(\rho, \sigma) = frac{1}{2} \lVert \rho - \sigma \rVert_1``
 
-
+Note that `A` and `B` must be either [`Ket`](@ref) or [`Operator`](@ref).
+"""
+tracedist(ρ::QuantumObject{<:AbstractArray{T1},ObjType1}, σ::QuantumObject{<:AbstractArray{T2},ObjType2}) where {T1,T2,ObjType1<:Union{KetQuantumObject,OperatorQuantumObject},ObjType2<:Union{KetQuantumObject,OperatorQuantumObject}} = norm(ket2dm(ρ) - ket2dm(σ), 1) / 2
 
 function _ptrace_ket(QO::AbstractArray{T1}, dims::Vector{<:Integer}, sel::Vector{T2}) where
     {T1,T2<:Integer}

src/negativity.jl

@@ -47,7 +47,7 @@ function negativity(ρ::QuantumObject, subsys::Int; logarithmic::Bool=false)
     end
 
     ρ_pt = partial_transpose(ρ, mask)
-    tr_norm = sum(sqrt, eigvals(ρ_pt' * ρ_pt))
+    tr_norm = norm(ρ_pt, 1) # trace norm
 
     if logarithmic
         return log2(tr_norm)
@@ -89,7 +89,7 @@ function _partial_transpose(ρ::QuantumObject{<:AbstractArray, OperatorQuantumOb
         [pt_dims[n, 3 - mask2[n]] for n in 1:nsys]  # opposite value in mask2 (1 -> 2, and 2 -> 1)
     ]
     return QuantumObject(
-        reshape(PermutedDimsArray(reshape(ρ.data, (ρ.dims..., ρ.dims...)), pt_idx), size(ρ)),
+        reshape(permutedims(reshape(ρ.data, (ρ.dims..., ρ.dims...)), pt_idx), size(ρ)),
         Operator,
         ρ.dims
     )

src/quantum_object.jl

@@ -458,9 +458,21 @@ julia> tr(a' * a)
 LinearAlgebra.tr(A::QuantumObject{<:AbstractArray{T},OpType}) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} = tr(A.data)
 
 """
-    norm(A::QuantumObject)
+    svdvals(A::QuantumObject)
 
-Returns the norm of `A`.
+Return the singular values of a [`QuantumObject`](@ref) in descending order
+"""
+LinearAlgebra.svdvals(A::QuantumObject{<:AbstractVector}) = svdvals(A.data)
+LinearAlgebra.svdvals(A::QuantumObject{<:DenseMatrix}) = svdvals(A.data)
+LinearAlgebra.svdvals(A::QuantumObject{<:AbstractSparseMatrix}) = svdvals(sparse_to_dense(A.data))
+
+"""
+    norm(A::QuantumObject, p::Real=2)
+
+If `A` is either [`Ket`](@ref), [`Bra`](@ref), [`OperatorKet`](@ref), or [`OperatorBra`](@ref), returns the standard vector `p`-norm of `A`.
+If `A` is either [`Operator`](@ref) or [`SuperOperator`](@ref), returns [Schatten](https://en.wikipedia.org/wiki/Schatten_norm) `p`-norm of `A`.
+
+Note that the default value of `p=2`
 
 # Examples
 
@@ -483,7 +495,11 @@ julia> norm(ψ)
 1.0
 ```
 """
-LinearAlgebra.norm(A::QuantumObject{<:AbstractArray{T},OpType}) where {T,OpType<:QuantumObjectType} = norm(A.data)
+LinearAlgebra.norm(A::QuantumObject{<:AbstractArray{T},OpType}, p::Real=2) where {T,OpType<:Union{KetQuantumObject,BraQuantumObject,OperatorKetQuantumObject,OperatorBraQuantumObject}} = norm(A.data, p)
+function LinearAlgebra.norm(A::QuantumObject{<:AbstractArray{T},OpType}, p::Real=2) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
+    p == 2.0 && return norm(A.data, 2)
+    return norm(svdvals(A), p)
+end
 LinearAlgebra.normalize(A::QuantumObject{<:AbstractArray{T},OpType}) where {T,OpType<:QuantumObjectType} =
     QuantumObject(normalize(A.data), OpType(), A.dims)
 LinearAlgebra.normalize!(A::QuantumObject{<:AbstractArray{T},OpType}) where {T,OpType<:QuantumObjectType} =

test/quantum_objects.jl

@@ -211,10 +211,32 @@
     ψ2_size = size(ψ2)
     @test opstring == "Quantum Object:   type=OperatorBra   dims=$ψ2_dims   size=$ψ2_size\n$datastring"
 
+    # get coherence
     ψ = coherent(30, 3)
     α, δψ = get_coherence(ψ)
     @test isapprox(abs(α), 3, atol=1e-5) && abs2(δψ[1]) > 0.999
 
+    # svdvals, Schatten p-norm
+    vd = Qobj(  rand(ComplexF64, 10))
+    vs = Qobj(sprand(ComplexF64, 100, 0.1))
+    Md = Qobj(  rand(ComplexF64, 10, 10))
+    Ms = Qobj(sprand(ComplexF64, 10, 10, 0.5))
+    @test svdvals(vd)[1] ≈ √(vd' * vd)
+    @test svdvals(vs)[1] ≈ √(vs' * vs)
+    @test norm(Md, 1) ≈ sum(sqrt, abs.(eigenenergies(Md' * Md))) atol=1e-6
+    @test norm(Ms, 1) ≈ sum(sqrt, abs.(eigenenergies(Ms' * Ms))) atol=1e-6
+
+    # trace distance
+    ψz0 = basis(2, 0)
+    ψz1 = basis(2, 1)
+    ρz0 = dense_to_sparse(ket2dm(ψz0))
+    ρz1 = dense_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.
+    @test tracedist(ψz1, ρz0) ≈ 1.
+    @test tracedist(ρz0, ρz1) ≈ 1.
+
     # Broadcasting
     a = destroy(20)
     for op in ((+), (-), (*), (^))