Cleanup, TODOs, and tests

Author: akirakyle

src/superoperators.jl

@@ -319,62 +319,13 @@ end
 end
 
 # TODO should all of PauliTransferMatrix, ChiMatrix, ChoiState, and KrausOperators subclass AbstractSuperOperator?
-"""
-    KrausOperators(B1, B2, data)
-
-Superoperator represented as a list of Kraus operators.
-Note unlike the SuperOperator or ChoiState types where
-its possible to have basis_l[1] != basis_l[2] and basis_r[1] != basis_r[2]
-which allows representations of maps between general linear operators defined on H_A \to H_B,
-a quantum channel can only act on valid density operators which live in H_A \to H_A.
-Thus the Kraus representation is only defined for quantum channels which map
-(H_A \to H_A) \to (H_B \to H_B).
-    """
-mutable struct KrausOperators{B1,B2,T} <: AbstractSuperOperator{B1,B2}
-    basis_l::B1
-    basis_r::B2
-    data::T
-    function KrausOperators{BL,BR,T}(basis_l::BL, basis_r::BR, data::T) where {BL,BR,T}
-        if (any(!samebases(basis_r, M.basis_r) for M in data) ||
-            any(!samebases(basis_l, M.basis_l) for M in data))
-            throw(DimensionMismatch("Tried to assign data with incompatible bases"))
-        end
-
-        new(basis_l, basis_r, data)
-    end
-end
-KrausOperators{BL,BR}(b1::BL,b2::BR,data::T) where {BL,BR,T} = KrausOperators{BL,BR,T}(b1,b2,data)
-KrausOperators(b1::BL,b2::BR,data::T) where {BL,BR,T} = KrausOperators{BL,BR,T}(b1,b2,data)
-KrausOperators(b,data) = KrausOperators(b,b,data)
-tensor(a::KrausOperators, b::KrausOperators) =
-    KrausOperators(a.basis_l ⊗ b.basis_l, a.basis_r ⊗ b.basis_r,
-                   [A ⊗ B for A in a.data for B in b.data])
-dagger(a::KrausOperators) = KrausOperators(a.basis_r, a.basis_l, [dagger(op) for op in a.data])
-*(a::KrausOperators{B1,B2}, b::KrausOperators{B2,B3}) where {B1,B2,B3} =
-    KrausOperators(a.basis_l, b.basis_r, [A*B for A in a.data for B in b.data])
-*(a::KrausOperators, b::KrausOperators) = throw(IncompatibleBases())
-*(a::KrausOperators{BL,BR}, b::Operator{BR,BR}) where {BL,BR} =
-    sum(op*b*dagger(op) for op in a.data)
-
-function minimize_kraus_rank(kraus; tol=1e-9)
-    bl, br = kraus.basis_l, kraus.basis_r
-    dim = length(bl)*length(br)
-
-    A = stack(reshape(op.data, dim) for op in kraus.data; dims=1)
-    F = qr(A; tol=tol)
-    rank = maximum(findnz(F.R)[1]) # rank(F) doesn't work but should
-    #@assert (F.R ≈ (sparse(F.Q') * A[F.prow,F.pcol])[1:dim,:])
-    #@assert (all(iszero(F.R[rank+1:end,:])))
-
-    ops = [Operator(bl, br, copy(reshape( # copy materializes reshaped view
-        F.R[i,invperm(F.pcol)], (length(bl), length(br))))) for i=1:rank]
-    return KrausOperators(bl, br, ops)
-end
-
 """
     ChoiState <: AbstractSuperOperator
 
 Superoperator represented as a choi state.
+
+The convention is chosen such that the input operators live in `(basis_l[1], basis_r[1])` while
+the output operators live in `(basis_r[2], basis_r[2])`.
 """
 mutable struct ChoiState{B1,B2,T} <: AbstractSuperOperator{B1,B2}
     basis_l::B1
@@ -399,6 +350,9 @@ dagger(a::ChoiState) = ChoiState(dagger(SuperOperator(a)))
 *(a::ChoiState, b::Operator) = SuperOperator(a)*b
 ==(a::ChoiState, b::ChoiState) = (SuperOperator(a) == SuperOperator(b))
 
+# TOOD: decide whether to document and export this
+choi_to_operator(c::ChoiState) = Operator(c.basis_l[1]⊗c.basis_l[2], c.basis_r[1]⊗c.basis_r[2], c.data)
+
 # reshape swaps within systems due to colum major ordering
 # https://docs.qojulia.org/quantumobjects/operators/#tensor_order
 function _super_choi((l1, l2), (r1, r2), data)
@@ -426,50 +380,133 @@ end
 ChoiState(op::SuperOperator) = ChoiState(_super_choi(op.basis_l, op.basis_r, op.data)...)
 SuperOperator(op::ChoiState) = SuperOperator(_super_choi(op.basis_l, op.basis_r, op.data)...)
 
-*(a::ChoiState, b::SuperOperator) = SuperOperator(a)*b
-*(a::SuperOperator, b::ChoiState) = a*SuperOperator(b)
+
+"""
+    KrausOperators <: AbstractSuperOperator
+
+Superoperator represented as a list of Kraus operators.
+Note unlike the SuperOperator or ChoiState types where it is possible to have
+`basis_l[1] != basis_l[2]` and `basis_r[1] != basis_r[2]`
+which allows representations of maps between general linear operators defined on ``H_A \\to H_B``,
+a quantum channel can only act on valid density operators which live in ``H_A \\to H_A``.
+Thus the Kraus representation is only defined for quantum channels which map
+``(H_A \\to H_A) \\to (H_B \\to H_B)``.
+"""
+mutable struct KrausOperators{B1,B2,T} <: AbstractSuperOperator{B1,B2}
+    basis_l::B1
+    basis_r::B2
+    data::T
+    function KrausOperators{BL,BR,T}(basis_l::BL, basis_r::BR, data::T) where {BL,BR,T}
+        if (any(!samebases(basis_r, M.basis_r) for M in data) ||
+            any(!samebases(basis_l, M.basis_l) for M in data))
+            throw(DimensionMismatch("Tried to assign data with incompatible bases"))
+        end
+
+        new(basis_l, basis_r, data)
+    end
+end
+KrausOperators{BL,BR}(b1::BL,b2::BR,data::T) where {BL,BR,T} = KrausOperators{BL,BR,T}(b1,b2,data)
+KrausOperators(b1::BL,b2::BR,data::T) where {BL,BR,T} = KrausOperators{BL,BR,T}(b1,b2,data)
+
+tensor(a::KrausOperators, b::KrausOperators) =
+    KrausOperators(a.basis_l ⊗ b.basis_l, a.basis_r ⊗ b.basis_r,
+                   [A ⊗ B for A in a.data for B in b.data])
+dagger(a::KrausOperators) = KrausOperators(a.basis_r, a.basis_l, [dagger(op) for op in a.data])
+*(a::KrausOperators{B1,B2}, b::KrausOperators{B2,B3}) where {B1,B2,B3} =
+    KrausOperators(a.basis_l, b.basis_r, [A*B for A in a.data for B in b.data])
+*(a::KrausOperators, b::KrausOperators) = throw(IncompatibleBases())
+*(a::KrausOperators{BL,BR}, b::Operator{BR,BR}) where {BL,BR} = sum(op*b*dagger(op) for op in a.data)
+
+"""
+    canonicalize(kraus::KrausOperators; tol=1e-9)
+
+Canonicalize the set kraus operators by performing a qr decomposition.
+A quantum channel with kraus operators ``{A_k}`` is in cannonical form if and only if
+
+```math
+\\Tr A_i^\\dagger A_j \\sim \\delta_{i,j}
+```
+
+If the input dimension is d and output dimension is d' then the number of kraus
+operators returned is guaranteed to be no greater than dd' and will furthermore
+be equal the Kraus rank of the channel up to numerical imprecision controlled by `tol`.  
+"""
+function canonicalize(kraus::KrausOperators; tol=1e-9)
+    bl, br = kraus.basis_l, kraus.basis_r
+    dim = length(bl)*length(br)
+
+    A = stack(reshape(op.data, dim) for op in kraus.data; dims=1)
+    F = qr(A; tol=tol)
+    # rank(F) for some reason doesn't work but should
+    rank = maximum(findnz(F.R)[1]) 
+    # Sanity checks that help illustrate what qr() returns:
+    # @assert (F.R ≈ (sparse(F.Q') * A[F.prow,F.pcol])[1:dim,:])
+    # @assert (all(iszero(F.R[rank+1:end,:])))
+
+    ops = [Operator(bl, br, copy(reshape( # copy materializes reshaped view
+        F.R[i,invperm(F.pcol)], (length(bl), length(br))))) for i=1:rank]
+    return KrausOperators(bl, br, ops)
+end
+
+# TODO: check if canonicalize and orthogonalize are equivalent
+orthogonalize(kraus::KrausOperators; tol=1e-9) = KrausOperators(ChoiState(kraus); tol=tol)
 
 SuperOperator(kraus::KrausOperators) =
     SuperOperator((kraus.basis_l, kraus.basis_l), (kraus.basis_r, kraus.basis_r),
                   (sum(conj(op)⊗op for op in kraus.data)).data)
 
-ChoiState(kraus::KrausOperators; tol=1e-9) =
+ChoiState(kraus::KrausOperators) =
     ChoiState((kraus.basis_r, kraus.basis_l), (kraus.basis_r, kraus.basis_l),
               (sum((M=op.data; reshape(M, (length(M), 1))*reshape(M, (1, length(M))))
-                   for op in kraus.data)); tol=tol)
+                   for op in kraus.data)))
 
-function KrausOperators(choi::ChoiState; tol=1e-9)
+function KrausOperators(choi::ChoiState; tol=1e-9, warn=true)
     if (!samebases(choi.basis_l[1], choi.basis_r[1]) ||
         !samebases(choi.basis_l[2], choi.basis_r[2]))
-        throw(DimensionMismatch("Tried to convert choi state of something that isn't a quantum channel mapping density operators to density operators"))
+        throw(DimensionMismatch("Tried to convert Choi state of something that isn't a quantum channel mapping density operators to density operators"))
+    end
+    # TODO: consider using https://github.com/jlapeyre/IsApprox.jl
+    if !ishermitian(choi.data) || !isapprox(choi.data, choi.data', atol=tol)
+        throw(ArgumentError("Tried to convert nonhermitian Choi state"))
     end
     bl, br = choi.basis_l[2], choi.basis_l[1]
-    #ishermitian(choi.data) || @warn "ChoiState is not hermitian"
     # TODO: figure out how to do this with sparse matrices using e.g. Arpack.jl or ArnoldiMethod.jl
     vals, vecs = eigen(Hermitian(Matrix(choi.data)))
     for val in vals
-        (abs(val) > tol && val < 0) && @warn "eigval $(val) < 0 but abs(eigval) > tol=$(tol)"
+        if warn && (abs(val) > tol && val < 0)
+            @warn "eigval $(val) < 0 but abs(eigval) > tol=$(tol)"
+        end
     end
     ops = [Operator(bl, br, sqrt(val)*reshape(vecs[:,i], length(bl), length(br)))
            for (i, val) in enumerate(vals) if abs(val) > tol && val > 0]
     return KrausOperators(bl, br, ops)
 end
 
-KrausOperators(op::SuperOperator; tol=1e-9) = KrausOperators(ChoiState(op; tol=tol); tol=tol)
+KrausOperators(op::SuperOperator; tol=1e-9) = KrausOperators(ChoiState(op); tol=tol)
 
+# TODO: document superoperator representation precident: everything of mixed type returns SuperOperator
 *(a::ChoiState, b::SuperOperator) = SuperOperator(a)*b
 *(a::SuperOperator, b::ChoiState) = a*SuperOperator(b)
 *(a::KrausOperators, b::SuperOperator) = SuperOperator(a)*b
 *(a::SuperOperator, b::KrausOperators) = a*SuperOperator(b)
 *(a::KrausOperators, b::ChoiState) = SuperOperator(a)*SuperOperator(b)
 *(a::ChoiState, b::KrausOperators) = SuperOperator(a)*SuperOperator(b)
 
-function is_trace_preserving(kraus::KrausOperators; tol=1e-9)
-    m = I(length(kraus.basis_r)) - sum(dagger(M)*M for M in kraus.data).data
-    m[abs.(m) .< tol] .= 0
-    return iszero(m)
+# TODO: document this
+is_trace_preserving(kraus::KrausOperators; tol=1e-9) =
+    isapprox.(I(length(kraus.basis_r)) - sum(dagger(M)*M for M in kraus.data).data, atol=tol)
+
+function is_trace_preserving(choi::ChoiState; tol=1e-9)
+    if (!samebases(choi.basis_l[1], choi.basis_r[1]) ||
+        !samebases(choi.basis_l[2], choi.basis_r[2]))
+        throw(DimensionMismatch("Choi state is of something that isn't a quantum channel mapping density operators to density operators"))
+    end
+    bl, br = choi.basis_l[2], choi.basis_l[1]
+    return isapprox.(I(length(br)) - ptrace(choi_to_operator(choi), 2), atol=tol)
 end
 
+is_trace_preserving(super::SuperOperator; tol=1e-9) = is_trace_preserving(ChoiState(super); tol=tol)
+
 # this check seems suspect... since it fails while the below check on choi succeeeds
 #function is_valid_channel(kraus::KrausOperators; tol=1e-9)
 #    m = I(length(kraus.basis_r)) - sum(dagger(M)*M for M in kraus.data).data
@@ -484,8 +521,20 @@ end
 #    return iszero(eigs)
 #end
 
-function is_valid_channel(choi::ChoiState; tol=1e-9)
+# TODO: pull out check for valid quantum channel
+# TODO: document this
+function is_completely_positive(choi::ChoiState; tol=1e-9)
+    if (!samebases(choi.basis_l[1], choi.basis_r[1]) ||
+        !samebases(choi.basis_l[2], choi.basis_r[2]))
+        throw(DimensionMismatch("Choi state is of something that isn't a quantum channel mapping density operators to density operators"))
+    end
     any(abs.(choi.data - choi.data') .> tol) && return false
     eigs = eigvals(Hermitian(Matrix(choi.data)))
     return all(@. abs(eigs) < tol || eigs > 0)
 end
+
+is_completely_positive(kraus::KrausOperators; tol=1e-9) = is_completely_positive(ChoiState(kraus); tol=tol)
+is_completely_positive(super::SuperOperator; tol=1e-9) = is_completely_positive(ChoiState(super); tol=tol)
+
+# TODO: document this
+is_cptp(sop) = is_completly_positive(sop) && is_trace_preserving(sop)

test/test_superoperators.jl

@@ -247,20 +247,64 @@ N = exp(log(2) * sparse(L)) # 50% loss channel
 @test (0.5 - real(tr(ρ^2))) < 1e-10 # one photon state becomes maximally mixed
 @test tracedistance(ρ, normalize(dm(fockstate(b, 0)) + dm(fockstate(b, 1)))) < 1e-10
 
-# Testing 0-2-4 binomial code encoder
+# 0-2-4 binomial code encoder
 b_logical = SpinBasis(1//2)
 b_fock = FockBasis(5)
 z_l = normalize(fockstate(b_fock, 0) + fockstate(b_fock, 4))
 o_l = fockstate(b_fock, 2)
-encoder_kraus = z_l ⊗ dagger(spinup(b_logical)) + o_l ⊗ dagger(spindown(b_logical))
-    encoder_sup = sprepost(encoder_kraus, dagger(encoder_kraus))
-decoder_sup = sprepost(dagger(encoder_kraus), encoder_kraus)
-@test SuperOperator(ChoiState(encoder_sup)).data == encoder_sup.data
-@test decoder_sup == dagger(encoder_sup)
-@test ChoiState(decoder_sup) == dagger(ChoiState(encoder_sup))
-@test decoder_sup*encoder_sup ≈ dense(identitysuperoperator(b_logical))
-@test decoder_sup*ChoiState(encoder_sup) ≈ dense(identitysuperoperator(b_logical))
-@test ChoiState(decoder_sup)*encoder_sup ≈ dense(identitysuperoperator(b_logical))
-@test SuperOperator(ChoiState(decoder_sup)*ChoiState(encoder_sup)) ≈ dense(identitysuperoperator(b_logical))
+enc_proj = z_l ⊗ dagger(spinup(b_logical)) + o_l ⊗ dagger(spindown(b_logical))
+dec_proj = dagger(enc_proj)
+enc_sup = sprepost(enc_proj, dec_proj)
+dec_sup = sprepost(dec_proj, enc_proj)
+enc_kraus = KrausOperators(b_fock, b_logical, [enc_proj])
+dec_kraus = KrausOperators(b_logical, b_fock, [dec_proj])
+## testing conversions
+@test dec_sup == dagger(enc_sup)
+@test dec_kraus == dagger(enc_kraus)
+@test ChoiState(enc_sup) == ChoiState(enc_kraus)
+@test ChoiState(dec_sup) == dagger(ChoiState(enc_sup))
+@test ChoiState(dec_kraus) == dagger(ChoiState(enc_kraus))
+@test SuperOperator(ChoiState(enc_sup)) == enc_sup
+@test SuperOperator(KrausOperators(enc_sup)) == enc_sup
+@test KrausOperators(ChoiState(enc_kraus)) == enc_kraus
+@test KrausOperators(SuperOperator(enc_kraus)) == enc_kraus
+## testing multipication
+@test dec_sup*enc_sup ≈ dense(identitysuperoperator(b_logical))
+@test SuperOperator(dec_kraus*enc_kraus) ≈ dense(identitysuperoperator(b_logical))
+@test dec_sup*enc_kraus ≈ dense(identitysuperoperator(b_logical))
+@test dec_kraus*enc_sup ≈ dense(identitysuperoperator(b_logical))
+@test SuperOperator(dec_kraus*enc_kraus) ≈ dense(identitysuperoperator(b_logical))
+@test dec_sup*ChoiState(enc_sup) ≈ dense(identitysuperoperator(b_logical))
+@test ChoiState(dec_sup)*enc_sup ≈ dense(identitysuperoperator(b_logical))
+@test SuperOperator(ChoiState(dec_sup)*ChoiState(enc_sup)) ≈ dense(identitysuperoperator(b_logical))
+
+@test avg_gate_fidelity(dec_sup*enc_sup) ≈ 1
+@test avg_gate_fidelity(dec_kraus*enc_kraus) ≈ 1
+@test avg_gate_fidelity(ChoiState(dec_sup)*ChoiState(enc_sup)) ≈ 1
+
+# test amplitude damping channel
+function amplitude_damp_kraus_op(b, γ, l)
+    op = SparseOperator(b)
+    for n=l:(length(b)-1)
+        op.data[n-l+1, n+1] = sqrt(binomial(n,l) * (1-γ)^(n-l) * γ^l)
+    end
+    return op
+end
+
+function test_kraus_channel(N, γ, tol)
+    b = FockBasis(N)
+    super = exp(liouvillian(identityoperator(b), [destroy(b)]))
+    kraus = KrausOperators(b, b, [amplitude_damp_kraus_op(b, γ, i) for i=0:(N-1)])
+    @test SuperOperator(kraus) ≈ super
+    @test ChoiState(kraus) ≈ ChoiState(super)
+    @test kraus ≈ KrausOperators(super; tol=tol)
+    @test is_trace_preserving(kraus; tol=tol)
+    @test is_valid_channel(kraus; tol=tol)
+end
+
+test_kraus_channel(10, 0.1, 1e-8)
+test_kraus_channel(20, 0.1, 1e-8)
+test_kraus_channel(10, 0.01, 1e-8)
+test_kraus_channel(20, 0.01, 1e-8)
 
 end # testset