Skip to content

Commit

Permalink
No more sparse superoperator eigen decomposition
Browse files Browse the repository at this point in the history
  • Loading branch information
@akirakyle
akirakyle committed Nov 19, 2024
1 parent f928ecb commit 3868939
Show file tree
Hide file tree
Showing 3 changed files with 39 additions and 108 deletions.
4 changes: 0 additions & 4 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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"
Expand Down
138 changes: 38 additions & 100 deletions src/superoperators.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)

Expand All @@ -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:
Expand All @@ -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)))
Expand All @@ -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)"))
Expand All @@ -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
Expand All @@ -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)

5 changes: 1 addition & 4 deletions test/test_superoperators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit 3868939

Please sign in to comment.