From c12a1505bd9f780f422c1381f4bee5c0b717758a Mon Sep 17 00:00:00 2001 From: Akira Kyle Date: Fri, 19 Jul 2024 15:18:56 -0600 Subject: [PATCH] Implement KrausOperators, conversions, is_trace_preserving, is_valid_channel --- src/QuantumOpticsBase.jl | 2 +- src/superoperators.jl | 89 ++++++++++++++++++++++++++++++++++++---- 2 files changed, 81 insertions(+), 10 deletions(-) diff --git a/src/QuantumOpticsBase.jl b/src/QuantumOpticsBase.jl index 183e23a..747ab41 100644 --- a/src/QuantumOpticsBase.jl +++ b/src/QuantumOpticsBase.jl @@ -36,7 +36,7 @@ export Basis, GenericBasis, CompositeBasis, basis, #superoperators SuperOperator, DenseSuperOperator, DenseSuperOpType, SparseSuperOperator, SparseSuperOpType, - ChoiState, + ChoiState, KrausOperators, is_trace_preserving, is_valid_channel, spre, spost, sprepost, liouvillian, identitysuperoperator, #fock FockBasis, number, destroy, create, diff --git a/src/superoperators.jl b/src/superoperators.jl index 0bb5271..f964e99 100644 --- a/src/superoperators.jl +++ b/src/superoperators.jl @@ -11,8 +11,9 @@ mutable struct SuperOperator{B1,B2,T} <: AbstractSuperOperator{B1,B2} basis_r::B2 data::T function SuperOperator{BL,BR,T}(basis_l::BL, basis_r::BR, data::T) where {BL,BR,T} - if length(basis_l[1])*length(basis_l[2]) != size(data, 1) || - length(basis_r[1])*length(basis_r[2]) != size(data, 2) + if (length(basis_l) != 2 || length(basis_r) != 2 || + length(basis_l[1])*length(basis_l[2]) != size(data, 1) || + length(basis_r[1])*length(basis_r[2]) != size(data, 2)) throw(DimensionMismatch("Tried to assign data of size $(size(data)) to Hilbert spaces of sizes $(length.(basis_l)), $(length.(basis_r))")) end new(basis_l, basis_r, data) @@ -318,26 +319,72 @@ end end # TODO should all of PauliTransferMatrix, ChiMatrix, ChoiState, and KrausOperators subclass AbstractSuperOperator? + """ - Base class for the Choi representation of superoperators. + 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) + +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) +end + +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 + eigs = eigvals(Matrix(m)) + eigs[@. abs(eigs) < tol || eigs > 0] .= 0 + return iszero(eigs) +end + """ + ChoiState(B1, B2, data) +Superoperator represented as a choi state stored as a sparse or dense matrix. +""" mutable struct ChoiState{B1,B2,T} <: AbstractSuperOperator{B1,B2} basis_l::B1 basis_r::B2 data::T - function ChoiState{BL,BR,T}(basis_l::BL, basis_r::BR, data::T) where {BL,BR,T} + function ChoiState{BL,BR,T}(basis_l::BL, basis_r::BR, data::T; tol=1e-9) where {BL,BR,T} if (length(basis_l) != 2 || length(basis_r) != 2 || length(basis_l[1])*length(basis_l[2]) != size(data, 1) || length(basis_r[1])*length(basis_r[2]) != size(data, 2)) throw(DimensionMismatch("Tried to assign data of size $(size(data)) to Hilbert spaces of sizes $(length.(basis_l)), $(length.(basis_r))")) end - new(basis_l, basis_r, data) + if any(abs.(data - data') .> tol) + @warn "Trying to construct ChoiState from non-hermitian data" + end + new(basis_l, basis_r, Hermitian(data)) end end -ChoiState{BL,BR}(b1::BL,b2::BR,data::T) where {BL,BR,T} = ChoiState{BL,BR,T}(b1,b2,data) -ChoiState(b1::BL,b2::BR,data::T) where {BL,BR,T} = ChoiState{BL,BR,T}(b1,b2,data) -ChoiState(b,data) = ChoiState(b,b,data) +ChoiState{BL,BR}(b1::BL,b2::BR,data::T; tol=1e-9) where {BL,BR,T} = ChoiState{BL,BR,T}(b1,b2,data;tol=tol) +ChoiState(b1::BL,b2::BR,data::T; tol=1e-9) where {BL,BR,T} = ChoiState{BL,BR,T}(b1,b2,data; tol=tol) +ChoiState(b,data; tol=tol) = ChoiState(b,b,data; tol=tol) # TODO: document why we have super_to_choi return non-trace one density matrices. # https://forest-benchmarking.readthedocs.io/en/latest/superoperator_representations.html @@ -367,5 +414,29 @@ function _super_choi((r2, l2), (r1, l1), data::SparseMatrixCSC) return (l1, l2), (r1, r2), sparse(data) end -ChoiState(op::SuperOperator) = ChoiState(_super_choi(op.basis_l, op.basis_r, op.data)...) +ChoiState(op::SuperOperator; tol=1e-9) = ChoiState(_super_choi(op.basis_l, op.basis_r, op.data)...; 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) + SuperOperator(op::ChoiState) = SuperOperator(_super_choi(op.basis_l, op.basis_r, op.data)...) +ChoiState(kraus::KrausOperators) = ChoiState(SuperOperator(kraus)) + +function KrausOperators(choi::ChoiState; tol=1e-9) + if (!samebases(choi.basis_l[1], choi.basis_l[2]) || + !samebases(choi.basis_r[1], 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")) + end + bl, br = choi.basis_l[1], choi.basis_r[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)" + 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)