diff --git a/src/superoperators.jl b/src/superoperators.jl index 205e5c5..15b376a 100644 --- a/src/superoperators.jl +++ b/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,37 +380,111 @@ 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 @@ -464,12 +492,21 @@ KrausOperators(op::SuperOperator; tol=1e-9) = KrausOperators(ChoiState(op; tol=t *(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) diff --git a/test/test_superoperators.jl b/test/test_superoperators.jl index 803a085..30832e2 100644 --- a/test/test_superoperators.jl +++ b/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