From 697c296e2ae3cd82df393e89de8c96c8ebb9de9f Mon Sep 17 00:00:00 2001 From: Akira Kyle Date: Tue, 29 Oct 2024 15:57:49 -0600 Subject: [PATCH] No more sparse superoperator eigen decomposition --- Project.toml | 4 -- src/superoperators.jl | 138 ++++++++++-------------------------- test/test_superoperators.jl | 5 +- 3 files changed, 39 insertions(+), 108 deletions(-) diff --git a/Project.toml b/Project.toml index a5db2c4..6dc04b6 100644 --- a/Project.toml +++ b/Project.toml @@ -4,12 +4,10 @@ version = "0.5.3" [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 9f7cb0e..0404c34 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 @@ -472,7 +470,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 @@ -482,9 +480,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) @@ -502,7 +501,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: @@ -515,12 +514,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))) @@ -540,93 +539,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)")) @@ -636,7 +570,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 @@ -646,48 +580,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 f0bc370..67c3cbd 100644 --- a/test/test_superoperators.jl +++ b/test/test_superoperators.jl @@ -323,10 +323,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