diff --git a/Project.toml b/Project.toml index f2f793b..4bd8172 100644 --- a/Project.toml +++ b/Project.toml @@ -4,12 +4,10 @@ version = "0.5.4" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FastExpm = "7868e603-8603-432e-a1a1-694bd70b01f2" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" -KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" QuantumInterface = "5717a53b-5d69-4fa3-b976-0bf2f97ca1e5" @@ -22,12 +20,10 @@ UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6" [compat] Adapt = "1, 2, 3.3, 4" -ArnoldiMethod = "0.4.0" FFTW = "1.2" FastExpm = "1.1.0" FastGaussQuadrature = "0.5, 1" FillArrays = "0.13, 1" -KrylovKit = "0.8.1" LRUCache = "1" LinearAlgebra = "1" QuantumInterface = "0.3.3" diff --git a/src/superoperators.jl b/src/superoperators.jl index f24b098..94cb26e 100644 --- a/src/superoperators.jl +++ b/src/superoperators.jl @@ -1,8 +1,6 @@ import Base: isapprox import QuantumInterface: AbstractSuperOperator import FastExpm: fastExpm -import KrylovKit: eigsolve -using ArnoldiMethod # TODO: this should belong in QuantumInterface.jl abstract type OperatorBasis{BL<:Basis,BR<:Basis} end @@ -471,7 +469,7 @@ tensor(a::KrausOperators, b::KrausOperators) = [A ⊗ B for A in a.data for B in b.data]) """ - orthogonalize(kraus::KrausOperators; tol=1e-12) + orthogonalize(kraus::KrausOperators; tol=√eps) Orthogonalize the set kraus operators by performing a qr decomposition on their vec'd operators. Note that this is different than `canonicalize` which returns a kraus decomposition such @@ -481,9 +479,10 @@ 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', however it may be greater than the Kraus rank. -`orthogonalize` should always be much faster than canonicalize as it avoids an explicit eigendecomposition. +`orthogonalize` should always be much faster than canonicalize as it avoids an explicit eigendecomposition +and thus also preserves sparsity if the kraus operators are sparse. """ -function orthogonalize(kraus::KrausOperators; tol=1e-12) +function orthogonalize(kraus::KrausOperators; tol=_get_tol(kraus)) bl, br = kraus.basis_l, kraus.basis_r dim = length(bl)*length(br) @@ -501,7 +500,7 @@ function orthogonalize(kraus::KrausOperators; tol=1e-12) end """ - canonicalize(kraus::KrausOperators; tol=1e-12) + canonicalize(kraus::KrausOperators; tol=√eps) Transform the quantum channel into canonical form such that the kraus operators ``{A_k}`` are Hilbert–Schmidt orthorgonal: @@ -514,12 +513,12 @@ 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`. """ -canonicalize(kraus::KrausOperators; tol=1e-12) = KrausOperators(ChoiState(kraus); tol=tol) +canonicalize(kraus::KrausOperators; tol=_get_tol(kraus)) = KrausOperators(ChoiState(kraus); tol=tol) # TODO: document -function make_trace_preserving(kraus; tol=1e-12) +function make_trace_preserving(kraus; tol=_get_tol(kraus)) m = I - sum(dagger(M)*M for M in kraus.data).data - if isa(_positive_eigen(m; tol=tol), Number) + if isa(_positive_eigen(m, tol), Number) throw(ArgumentError("Channel must be trace nonincreasing")) end K = Operator(kraus.basis_l, kraus.basis_r, sqrt(Matrix(m))) @@ -539,93 +538,28 @@ _choi_state_maps_density_ops(choi::ChoiState) = (samebases(choi.basis_l[1], choi samebases(choi.basis_l[2], choi.basis_r[2])) # TODO: consider using https://github.com/jlapeyre/IsApprox.jl -_is_hermitian(M; tol=1e-12) = ishermitian(M) || isapprox(M, M', atol=tol) -_is_identity(M; tol=1e-12) = isapprox(M, I, atol=tol) +_is_hermitian(M, tol) = ishermitian(M) || isapprox(M, M', atol=tol) +_is_identity(M, tol) = isapprox(M, I, atol=tol) # TODO: document # data must be Hermitian! -# performance of dense version typically faster until underlying hilbert spaces -# have dimension on the order 60 or so? -function _positive_eigen(data; tol=1e-12) - @assert ishermitian(data) +function _positive_eigen(data, tol) # LinearAlgebra's eigen returns eigenvals sorted smallest to largest for Hermitian matrices vals, vecs = eigen(Hermitian(Matrix(data))) - println("dense eigs: ", reverse(vals)) vals[1] < -tol && return vals[1] ret = [(val, vecs[:,i]) for (i, val) in enumerate(vals) if val > tol] - #for (val, vec) in ret - # println("dense eigs: ", val) - # vec = abs.(vec) - # vec[vec .< 1e-8] .= 0 - # println("vec: ", vec) - #end return ret end -# To control the precision of eigsolve, set KrylovDefaults.tol -# this isn't controlled by tol since we genally want KrolovKit to run -# with much higher precision so the eigenvectors we get are "good" -function _positive_eigen_kryolov(data::SparseMatrixCSC; tol=1e-12) - vals, vecs, info = eigsolve(Hermitian(data), 1, :SR) - info.converged < 1 && return -Inf - vals[1] < -tol && return vals[1] - vals, vecs, info = eigsolve(Hermitian(data), 1, :LM) - println("sparse eigs: ", reverse(vals)) - vals[end] > tol && (println("bailing..."); - return _positive_eigen(Matrix(data); tol=tol)) - ret = [(vals[i], vecs[i]) for i=length(vals):-1:1 if vals[i] > tol] - #for (val, vec) in ret - # println("sparse eigs: ", val) - # vec = abs.(vec) - # vec[vec .< 1e-8] .= 0 - # println("vec: ", vec) - #end - return ret -end - -function _positive_eigen(data::SparseMatrixCSC; tol=1e-12) - @assert ishermitian(data) - dim = size(data, 1) - F, info = partialschur(data; nev=1, which=:SR); - info.converged > 0 || return _positive_eigen(Matrix(data); tol=tol) - @assert abs(imag(F.eigenvalues[1])) < tol - real(F.eigenvalues[1]) < -tol && return real(F.eigenvalues[1]) - - nev = floor(Int, sqrt(dim)) - maxdim = dim < 16 ? dim : floor(Int, dim/4) - println("dim=", dim, ", maxdim=", maxdim,", nev=", nev) - #V, H = rand(eltype(data), size(data, 1), 2*nev+1), rand(eltype(data), 2*nev+1, 2*nev) - #V, H = rand(eltype(data), dim, dim + 1), rand(eltype(data), dim + 1, dim) - #arnoldi = ArnoldiWorkspace(V, H) - arnoldi = ArnoldiWorkspace(eltype(data), dim, maxdim) - start_from = 1 - - while true - F, info = partialschur!(data, arnoldi; nev=nev, mindim=nev, maxdim=maxdim, - start_from=start_from, which=:LM) - println("num eigs: ", length(F.eigenvalues), " info: ", info) - @assert all(@. abs(imag(F.eigenvalues)) < tol) - info.nconverged == nev && (real.(F.eigenvalues))[end] < tol && break - nev *= 2 - nev > maxdim && return (println("bailing..."); _positive_eigen(Matrix(data); tol=tol)) - start_from = info.nconverged+1 - end - - println("new sparse eigs: ", real.(F.eigenvalues)) - return [(real(F.eigenvalues[i]), F.Q[:,i]) for i=1:length(F.eigenvalues) - if real(F.eigenvalues[i]) > tol] -end - - -function KrausOperators(choi::ChoiState; tol=1e-12) +function KrausOperators(choi::ChoiState; tol=_get_tol(choi)) if !_choi_state_maps_density_ops(choi) throw(DimensionMismatch("Tried to convert Choi state of something that isn't a quantum channel mapping density operators to density operators")) end - if !_is_hermitian(choi.data; tol=tol) + if !_is_hermitian(choi.data, tol) throw(ArgumentError("Tried to convert nonhermitian Choi state")) end bl, br = choi.basis_l[2], choi.basis_l[1] - eigs = _positive_eigen(choi.data; tol=tol) + eigs = _positive_eigen(choi.data, tol) if isa(eigs, Number) throw(ArgumentError("Tried to convert a non-positive semidefinite Choi state,"* "failed for smallest eigval $(eigs), consider increasing tol=$(tol)")) @@ -635,7 +569,7 @@ function KrausOperators(choi::ChoiState; tol=1e-12) return KrausOperators(bl, br, ops) end -KrausOperators(op::SuperOperator; tol=1e-12) = KrausOperators(ChoiState(op); tol=tol) +KrausOperators(op::SuperOperator; tol=_get_tol(op)) = KrausOperators(ChoiState(op); tol=tol) # TODO: document superoperator representation precident: everything of mixed type returns SuperOperator *(a::ChoiState, b::SuperOperator) = SuperOperator(a)*b @@ -645,48 +579,52 @@ KrausOperators(op::SuperOperator; tol=1e-12) = KrausOperators(ChoiState(op); tol *(a::KrausOperators, b::ChoiState) = SuperOperator(a)*SuperOperator(b) *(a::ChoiState, b::KrausOperators) = SuperOperator(a)*SuperOperator(b) +_get_tol(kraus::KrausOperators) = sqrt(eps(real(eltype(eltype(fieldtypes(typeof(kraus))[3]))))) +_get_tol(super::SuperOperator) = sqrt(eps(real(eltype(fieldtypes(typeof(super))[3])))) +_get_tol(super::ChoiState) = sqrt(eps(real(eltype(fieldtypes(typeof(super))[3])))) + # TODO: document this -is_completely_positive(choi::KrausOperators; tol=1e-12) = true +is_completely_positive(choi::KrausOperators; tol=_get_tol(choi)) = true -function is_completely_positive(choi::ChoiState; tol=1e-12) +function is_completely_positive(choi::ChoiState; tol=_get_tol(choi)) _choi_state_maps_density_ops(choi) || return false - _is_hermitian(choi.data; tol=tol) || return false - isa(_positive_eigen(Hermitian(choi.data); tol=tol), Number) && return false + _is_hermitian(choi.data, tol) || return false + isa(_positive_eigen(Hermitian(choi.data), tol), Number) && return false return true end -is_completely_positive(super::SuperOperator; tol=1e-12) = +is_completely_positive(super::SuperOperator; tol=_get_tol(super)) = is_completely_positive(ChoiState(super); tol=tol) # TODO: document this -is_trace_preserving(kraus::KrausOperators; tol=1e-12) = - _is_identity(sum(dagger(M)*M for M in kraus.data).data, tol=tol) +is_trace_preserving(kraus::KrausOperators; tol=_get_tol(kraus)) = + _is_identity(sum(dagger(M)*M for M in kraus.data).data, tol) -is_trace_preserving(choi::ChoiState; tol=1e-12) = - _is_identity(ptrace(choi_to_operator(choi), 1).data, tol=tol) +is_trace_preserving(choi::ChoiState; tol=_get_tol(choi)) = + _is_identity(ptrace(choi_to_operator(choi), 1).data, tol) -is_trace_preserving(super::SuperOperator; tol=1e-12) = +is_trace_preserving(super::SuperOperator; tol=_get_tol(super)) = is_trace_preserving(ChoiState(super); tol=tol) # TODO: document this -function is_trace_nonincreasing(kraus::KrausOperators; tol=1e-12) +function is_trace_nonincreasing(kraus::KrausOperators; tol=_get_tol(kraus)) m = I - sum(dagger(M)*M for M in kraus.data).data - _is_hermitian(m; tol=tol) || return false - return !isa(_positive_eigen(Hermitian(m); tol=tol), Number) + _is_hermitian(m, tol) || return false + return !isa(_positive_eigen(Hermitian(m), tol), Number) end -function is_trace_nonincreasing(choi::ChoiState; tol=1e-12) +function is_trace_nonincreasing(choi::ChoiState; tol=_get_tol(choi)) m = I - ptrace(choi_to_operator(choi), 1).data - _is_hermitian(m; tol=tol) || return false - return !isa(_positive_eigen(Hermitian(m); tol=tol), Number) + _is_hermitian(m, tol) || return false + return !isa(_positive_eigen(Hermitian(m), tol), Number) end -is_trace_nonincreasing(super::SuperOperator; tol=1e-12) = +is_trace_nonincreasing(super::SuperOperator; tol=_get_tol(super)) = is_trace_nonincreasing(ChoiState(super); tol=tol) # TODO: document this -is_cptp(sop; tol=1e-12) = is_completely_positive(sop; tol=tol) && is_trace_preserving(sop; tol=tol) +is_cptp(sop; tol=_get_tol(sop)) = is_completely_positive(sop; tol=tol) && is_trace_preserving(sop; tol=tol) # TODO: document this -is_cptni(sop; tol=1e-12) = is_completely_positive(sop; tol=tol) && is_trace_nonincreasing(sop; tol=tol) +is_cptni(sop; tol=_get_tol(sop)) = is_completely_positive(sop; tol=tol) && is_trace_nonincreasing(sop; tol=tol) diff --git a/test/test_superoperators.jl b/test/test_superoperators.jl index 79201e3..0843b7f 100644 --- a/test/test_superoperators.jl +++ b/test/test_superoperators.jl @@ -324,10 +324,7 @@ function test_channel(κa, κp, N, do_sparse) @test is_cptp(kraus; tol=tol) @test isapprox(SuperOperator(kraus), super; atol=tol) @test isapprox(ChoiState(kraus), ChoiState(super); atol=tol) - #@test isapprox(kraus, KrausOperators(super; tol=tol); atol=tol) - - isapprox(kraus, KrausOperators(super; tol=tol); atol=tol) || println("isapprox kraus, ", "κa=", κa, ", κp=", κp, ", N=", N, ", sparse=", sparse) - #return kraus, super + @test isapprox(kraus, KrausOperators(super); atol=tol) end for N=1:4