Implement sparse _positive_eigen for conversion to Kraus

Project.toml

@@ -8,6 +8,7 @@ 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"
@@ -24,6 +25,7 @@ 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"

src/superoperators.jl

@@ -1,6 +1,21 @@
 import Base: isapprox
 import QuantumInterface: AbstractSuperOperator
 import FastExpm: fastExpm
+import KrylovKit: eigsolve
+
+# TODO: this should belong in QuantumInterface.jl
+abstract type OperatorBasis{BL<:Basis,BR<:Basis} end
+abstract type SuperOperatorBasis{BL<:OperatorBasis,BR<:OperatorBasis} end
+
+"""
+    tensor(E::AbstractSuperOperator, F::AbstractSuperOperator, G::AbstractSuperOperator...)
+
+Tensor product ``\\mathcal{E}⊗\\mathcal{F}⊗\\mathcal{G}⊗…`` of the given super operators.
+"""
+tensor(a::AbstractSuperOperator, b::AbstractSuperOperator) = arithmetic_binary_error("Tensor product", a, b)
+tensor(op::AbstractSuperOperator) = op
+tensor(operators::AbstractSuperOperator...) = reduce(tensor, operators)
+
 
 """
     SuperOperator <: AbstractSuperOperator
@@ -355,8 +370,25 @@ dagger(a::ChoiState) = ChoiState(dagger(SuperOperator(a)))
 ==(a::ChoiState, b::ChoiState) = (SuperOperator(a) == SuperOperator(b))
 isapprox(a::ChoiState, b::ChoiState; kwargs...) = isapprox(SuperOperator(a), SuperOperator(b); kwargs...)
 
-# TOOD: decide whether to document and export this
-choi_to_operator(c::ChoiState) = Operator(c.basis_l[2]⊗c.basis_l[1], c.basis_r[2]⊗c.basis_r[1], c.data)
+# Container to hold each of the four bases for a Choi operator when converting it to
+# an operator so that if any are CompositeBases tensor doesn't lossily collapse them
+struct ChoiSubBasis{S,B<:Basis} <: Basis
+    shape::S
+    basis::B
+end
+ChoiSubBasis(b::Basis) = ChoiSubBasis(b.shape, b)
+
+# TODO: decide whether to document and export this
+choi_to_operator(c::ChoiState) = Operator(
+    ChoiSubBasis(c.basis_l[2])⊗ChoiSubBasis(c.basis_l[1]), ChoiSubBasis(c.basis_r[2])⊗ChoiSubBasis(c.basis_r[1]), c.data)
+
+function tensor(a::ChoiState, b::ChoiState)
+    op = choi_to_operator(a) ⊗ choi_to_operator(b)
+    op = permutesystems(op, [1,3,2,4])
+    ChoiState((a.basis_l[1] ⊗ b.basis_l[1], a.basis_l[2] ⊗ b.basis_l[2]),
+              (a.basis_r[1] ⊗ b.basis_r[1], a.basis_r[2] ⊗ b.basis_r[2]), op.data)
+end
+tensor(a::SuperOperator, b::SuperOperator) = SuperOperator(tensor(ChoiState(a), ChoiState(b)))
 
 # reshape swaps within systems due to colum major ordering
 # https://docs.qojulia.org/quantumobjects/operators/#tensor_order
@@ -425,16 +457,18 @@ end
 KrausOperators{BL,BR}(b1::BL,b2::BR,data::Vector{T}) where {BL,BR,T} = KrausOperators{BL,BR,T}(b1,b2,data)
 KrausOperators(b1::BL,b2::BR,data::Vector{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])
+dense(a::KrausOperators) = KrausOperators(a.basis_l, a.basis_r, [dense(op) for op in a.data])
+sparse(a::KrausOperators) = KrausOperators(a.basis_l, a.basis_r, [sparse(op) for op in a.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)
 ==(a::KrausOperators, b::KrausOperators) = (SuperOperator(a) == SuperOperator(b))
 isapprox(a::KrausOperators, b::KrausOperators; kwargs...) = isapprox(SuperOperator(a), SuperOperator(b); kwargs...)
+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])
 
 """
     orthogonalize(kraus::KrausOperators; tol=1e-12)
@@ -508,18 +542,33 @@ _choi_state_maps_density_ops(choi::ChoiState) = (samebases(choi.basis_l[1], choi
 _is_hermitian(M; tol=1e-12) = ishermitian(M) || isapprox(M, M', atol=tol)
 _is_identity(M; tol=1e-12) = isapprox(M, I, atol=tol)
 
-# TODO: document
+# TODO: document, data must be Hermitian!
 function _positive_eigen(data; tol=1e-12)
-    # TODO: figure out how to do this with sparse matrices using e.g. Arpack.jl or ArnoldiMethod.jl
-    # I will want to run twice, first asking for smallest eigenvalue to check it is above -tol
-    # Then run a second time with asking for maybe sqrt(N) largest eigenvalues?
-    # If smallest of these is not smaller than tol, bail do dense method? 
     # LinearAlgebra's eigen returns eigenvals sorted smallest to largest for Hermitian matrices
     vals, vecs = eigen(Hermitian(Matrix(data)))
     vals[1] < -tol && return vals[1]
     return [(val, vecs[:,i]) for (i, val) in enumerate(vals) if val > tol]
 end
 
+# first asks for smallest eigenvalue to check it is above -tol
+# Then run a second time with asking for sqrt(N) largest eigenvalues
+# If smallest of these is not smaller than tol, bail to dense method 
+function _positive_eigen(data::SparseMatrixCSC; tol=1e-12)
+    vals, vecs, info = eigsolve(Hermitian(data), 1, :SR; tol=tol)
+    info.converged < 1 && return -Inf
+    vals[1] < -tol && return vals[1]
+    n = 1 + floor(Int, sqrt(size(data)[1]))
+    vals, vecs, info = eigsolve(Hermitian(data), n, :LM; tol=tol)
+    println(vals)
+    if info.converged < n || vals[n] > tol
+        println("bailed to dense method!!! n=", n, "  size(data)=", size(data))
+        return _positive_eigen(Matrix(data))
+    else
+        return [(val, vec) for (val, vec) in zip(vals, vecs) if val > tol]
+    end
+end
+
+
 function KrausOperators(choi::ChoiState; tol=1e-12)
     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"))

test/test_superoperators.jl

@@ -306,31 +306,41 @@ function phase_damp_kraus(γ)
     return KrausOperators(b,b,[K0, K1])
 end
 
-function test_channel(κa, κp)
+tensor_it(x, N) = tensor(fill(x, N)...)
+
+function test_channel(κa, κp, N, do_sparse)
     tol = 1e-7
     b = SpinBasis(1//2)
     La = liouvillian(identityoperator(b), [sigmam(b)]; rates=[κa])
     Lp = liouvillian(identityoperator(b), [sigmaz(b)]; rates=[κp])
-    super = exp(La + Lp)
+    super = tensor_it(exp(La + Lp), N)
+    super = do_sparse ? sparse(super) : dense(super)
     ka = ampl_damp_kraus(1-exp(-κa))
     kp = phase_damp_kraus(1-exp(-4*κp))
-    kraus = ka*kp
+    kraus = tensor_it(ka*kp, N)
+    kraus = do_sparse ? sparse(kraus) : dense(kraus)
+    println(typeof(super))
+    #println("is_cptp(super; tol=tol)")
     @test is_cptp(super; tol=tol)
+    #println("is_cptp(kraus;")
     @test is_cptp(kraus; tol=tol)
+    #println("isapprox(SuperOperator(kraus),")
     @test isapprox(SuperOperator(kraus), super; atol=tol)
+    #println("isapprox(ChoiState(kraus),")
     @test isapprox(ChoiState(kraus), ChoiState(super); atol=tol)
+    #println("isapprox(kraus,")
     @test isapprox(kraus, KrausOperators(super; tol=tol); atol=tol)
     return kraus, super
 end
 
-test_channel(1e-1, 0)
-test_channel(1e-2, 0)
-test_channel(1e-3, 0)
-test_channel(0, 1e-1)
-test_channel(0, 1e-2)
-test_channel(0, 1e-3)
-test_channel(1e-1, 1e-1)
-test_channel(1e-2, 1e-2)
-test_channel(1e-3, 1e-3)
+for N=1:2
+    for κa=[0, 1e-1, 1e-2, 1e-3]
+        for κp=[0, 1e-1, 1e-2, 1e-3]
+            test_channel(κa, κp, N, false)
+            #test_channel(κa, κp, N, true)
+        end
+    end
+end
+
 
 end # testset