More WIP Pauli

Author: akirakyle

src/QuantumOpticsBase.jl

@@ -8,7 +8,7 @@ import QuantumInterface: Basis, GenericBasis, CompositeBasis, basis, basis_l, ba
     IncompatibleBases, @compatiblebases, samebases, check_samebases,
     addible, check_addible, multiplicable, check_multiplicable, reduced, ptrace, permutesystems,
     dagger, directsum, ⊕, dm, embed, expect, identityoperator, identitysuperoperator,
-    permutesystems, projector, ptrace, reduced, tensor, ⊗, variance, apply!,
+    permutesystems, projector, ptrace, reduced, tensor, tensor_pow, ⊗, variance, apply!,
     vec, unvec, super, choi, kraus, stinespring, pauli, chi, spre, spost, sprepost, liouvillian
 
 # metrics
@@ -18,7 +18,7 @@ import QuantumInterface: entropy_vn, fidelity, logarithmic_negativity
 import QuantumInterface: complement, remove, shiftremove, reducedindices!, check_indices, check_sortedindices, check_embed_indices
 
 export Basis, GenericBasis, CompositeBasis, basis, basis_l, basis_r, dimension, shape,
-        tensor, ⊗, permutesystems, @compatiblebases,
+        tensor, tensor_pow, ⊗, permutesystems, @compatiblebases,
         #states
                 StateVector, Bra, Ket, basisstate, sparsebasisstate, norm,
                 dagger, normalize, normalize!,
@@ -98,7 +98,7 @@ include("particle.jl")
 include("nlevel.jl")
 include("manybody.jl")
 include("transformations.jl")
-#include("pauli.jl")
+include("pauli.jl")
 include("metrics.jl")
 include("spinors.jl")
 include("phasespace.jl")

src/pauli.jl

@@ -1,68 +1,55 @@
-using TransmuteDims
 
-#function _pauli_to_ketbra1()
-#    b = SpinBasis(1//2)
-#    vec_it(fn) = vec(fn(b).data)
-#    p2kb = hcat(map(vec_it, [identityoperator, sigmax, sigmay, sigmaz])...)
-#end
-#
-#function _ketbra_to_pauli(N)
-#    T = kron((_pauli_to_ketbra1 for _=1:N)...)
-#    T = reshape(T, Tuple(4 for _=1:2N))
-#    T = PermutedDimsArray(T, ((2i-1 for i=1:N)..., (2i for i=1:N)...))
-#    reshape(T, (4^N, 4^N))
-#end
-#
-#function _op_to_pauli(Nl, Nr, data)
-#    Ul = ket_bra_to_pauli(Nl)
-#    Ur = Nl == Nr ? Ul : _ketbra_to_pauli(Nr)
-#    data = dagger(Ul) * data * Ur # TODO figure out normalization
-#    @assert isapprox(imag.(data), zero(data), atol=tol)
-#    Operator(PauliBasis()^Nl, PauliBasis()^Nr, real.(data))
-#end
-
-function _pauli_to_comp1()
+# comp stands for computational basis
+function _comp_pauli_1(basis_fn)
     b = SpinBasis(1//2)
-    vec_it(fn) = vec(fn(b).data)
+    vec_it(fn) = vec(fn(b).data)/2 # provides "standard" normalization
     V = hcat(map(vec_it, [identityoperator, sigmax, sigmay, sigmaz])...)
-    Operator(PauliBasis(), KetBraBasis(b, b), V)
+    Operator(basis_fn(b, b), PauliBasis(1), V)
 end
 
-_pauli_to_comp1_cached = _pauli_to_comp1()
+_comp_pauli_kb1_cached = _comp_pauli_1(KetBraBasis)
+_comp_pauli_choi1_cached = _comp_pauli_1(ChoiBasis)
 
-# TODO should this be further cached?
+# TODO: should this be further cached?
 """
-    comp_to_pauli(N)
+    pauli_comp_kb(N)
 
-Creates a superoperator which changes from the computational `SpinBasis(1//2)`
+Creates a superoperator which changes from the computational `KetBra(SpinBasis(1//2))`
 to the `PauliBasis()` over `N` qubits.
 """
-function pauli_to_comp(N::Integer)
-    N > 0 || throw(ArgumentError())
-    N == 1 && return _pauli_to_comp1_cached
-    V = pauli_to_comp(N÷2)
-    (N%2 == 0) && return V⊗V
-    return V⊗V⊗pauli_to_comp(1)
-end
+comp_pauli_kb(N::Integer) = tensor_pow(_comp_pauli_kb1_cached, N)
+comp_pauli_choi(N::Integer) = tensor_pow(_comp_pauli_choi1_cached, N)
 
 function _check_is_spinbasis(b)
     for i=1:length(b)
         (b[i] isa SpinBasis && dimension(b[i]) == 2) || throw(ArgumentError("Superoperator must be over systems composed of SpinBasis(1//2) to be converted to pauli representation"))
     end
 end
 
+# It's possible to get better asympotic speedups using, e.g. methods from
+# https://iopscience.iop.org/article/10.1088/1402-4896/ad6499
+# https://arxiv.org/abs/2411.00526
+# https://quantum-journal.org/papers/q-2024-09-05-1461/ (see appendices)
+function pauli(op::Operator)
+    (basis_l(op) == basis_r(op)) || throw(ArgumentError("Operator must be square in order to be vectorized to the pauli represenation"))
+    _check_is_spinbasis(basis_l(op))
+    dagger(comp_pauli_kb(length(basis_l(op)))) * vec(op)
+end
+
 function pauli(op::SOpKetBraType)
-    ((basis_l(basis_l(op)) == basis_r(basis_l(op))) && (basis_l(br) == basis_r(basis_r(op)))) || throw(ArgumentError("Superoperator must map between square operators in order to be converted to pauli represenation"))
-    bl, br = basis_l(basis_l(op)), basis_l(basis_r(op))
+    bl, br = basis_l(op), basis_r(op)
+    ((basis_l(bl) == basis_r(bl)) && (basis_l(br) == basis_r(br))) || throw(ArgumentError("Superoperator must map between square operators in order to be converted to pauli represenation"))
+    bl, br = basis_l(bl), basis_l(br)
     foreach(_check_is_spinbasis, (bl, br))
 
-    Nl, Nr = length(basis_l(bl)), length(basis_l(br))
-    Vl = pauli_to_comp(Nl)
-    Vr = Nl == Nr ? Vl : pauli_to_comp(Nr)
+    Nl, Nr = length(bl), length(br)
+    Vl = comp_pauli_kb(Nl)
+    Vr = Nl == Nr ? Vl : comp_pauli_kb(Nr)
     #data = dagger(Ul)*op.data*Ur # TODO figure out normalization
     #@assert isapprox(imag.(data), zero(data), atol=tol)
     #Operator(PauliBasis()^Nl, PauliBasis()^Nr, real.(data))
     dagger(Vl) * op * Vr # TODO figure out normalization
+    # TODO make sure dagger here is okay...
 end
 
 function chi(op::ChoiStateType)
@@ -71,11 +58,9 @@ function chi(op::ChoiStateType)
     foreach(_check_is_spinbasis, (bl, br))
 
     Nl, Nr = length(bl), length(br)
-    Ul = ket_bra_to_pauli(Nl)
-    Ur = Nl == Nr ? Ul : ket_bra_to_pauli(Nr)
-    data = dagger(Ul)*op.data*Ur # TODO figure out normalization
-    @assert isapprox(imag.(data), zero(data), atol=tol)
-    Operator(PauliBasis()^Nl, PauliBasis()^Nr, real.(data))
+    Vl = comp_pauli_kb(Nl)
+    Vr = Nl == Nr ? Vl : comp_pauli_kb(Nr)
+    dagger(Vl) * op * Vr # TODO figure out normalization
 end
 
 """

src/superoperators.jl

@@ -17,15 +17,15 @@ const SparseSuperOpAdjType{BL,BR} = SuperOperatorType{BL,BR,<:Adjoint{<:Number,<
 const SparseSuperOpType{BL,BR} = Union{SparseOpPureType{BL,BR},SparseOpAdjType{BL,BR}}
 
 vec(op::Operator) = Ket(KetBraBasis(basis_l(op), basis_r(op)), vec(op.data))
-#function unvec(k::Ket{<:KetBraBasis})
-#    bl, br = basis_l(basis(k)), basis_r(basis(k))
-#    Operator(bl, br, reshape(k.data, dimension(bl), dimension(br)))
-#end
 function unvec(k::Ket{<:KetBraBasis})
     bl, br = basis_l(basis(k)), basis_r(basis(k))
-    @cast A[n,m] |= k.data[(n,m)] (n ∈ 1:dimension(bl), m ∈ 1:dimension(br))
-    return Operator(bl, br, A)
+    Operator(bl, br, reshape(k.data, dimension(bl), dimension(br)))
 end
+#function unvec(k::Ket{<:KetBraBasis})
+#    bl, br = basis_l(basis(k)), basis_r(basis(k))
+#    @cast A[n,m] |= k.data[(n,m)] (n ∈ 1:dimension(bl), m ∈ 1:dimension(br))
+#    return Operator(bl, br, A)
+#end
 
 function spre(op::Operator)
     multiplicable(op, op) || throw(ArgumentError("It's not clear what spre of a non-square operator should be. See issue #113"))
@@ -37,18 +37,22 @@ function spost(op::Operator)
     sprepost(identityoperator(op), op)
 end
 
-#sprepost(A::Operator, B::Operator) = Operator(KetBraBasis(basis_l(A), basis_r(B)), KetBraBasis(basis_r(A), basis_l(B)), kron(permutedims(B.data), A.data))
-function sprepost(A::Operator, B::Operator)
-    @cast C[(ν,μ), (n,m)] |= A.data[ν,n] * B.data[m,μ]
-    Operator(KetBraBasis(basis_l(A), basis_r(B)), KetBraBasis(basis_r(A), basis_l(B)), C)
-end
+sprepost(A::Operator, B::Operator) = Operator(KetBraBasis(basis_l(A), basis_r(B)), KetBraBasis(basis_r(A), basis_l(B)), kron(permutedims(B.data), A.data))
+#function sprepost(A::Operator, B::Operator)
+#    @cast C[(ν,μ), (n,m)] |= A.data[ν,n] * B.data[m,μ]
+#    Operator(KetBraBasis(basis_l(A), basis_r(B)), KetBraBasis(basis_r(A), basis_l(B)), C)
+#end
 
 function tensor(A::SuperOperatorType, B::SuperOperatorType)
     a1, a2 = basis_l(A), basis_r(A)
     b1, b2 = basis_l(B), basis_r(B)
-    da1, da2, db1, db2 = map(dimension, (a1, a2, b1, b2))
-    @cast A[(ν,μ), (n,m)] |= A[m,μ] * B[n,ν] (m ∈ 1:da1, μ ∈ 1:da2, n ∈ 1:db1, ν ∈ 1:db2)
-    return Operator(a1⊗b2, a2⊗b2, A)
+    #da1, da2, db1, db2 = map(dimension, (a1, a2, b1, b2))
+    #@cast C[(ν,μ), (n,m)] |= A.data[m,μ] * B.data[n,ν] (m ∈ 1:da1, μ ∈ 1:da2, n ∈ 1:db1, ν ∈ 1:db2)
+    data = kron(B.data, A.data)
+    data = reshape(data, map(dimension, (a1, a2, b1, b2)))
+    data = PermutedDimsArray(data, (4, 2, 3, 1))
+    data = reshape(data, map(dimension, (a2⊗b2, a1⊗b1)))
+    return Operator(a1⊗b1, a2⊗b2, data)
 end
 
 #function _sch(data, l1, l2, r1, r2)
@@ -79,13 +83,13 @@ end
 function _super_choi(basis_fn, op)
     l1, l2 = basis_l(basis_l(op)), basis_r(basis_l(op))
     r1, r2 = basis_l(basis_r(op)), basis_r(basis_r(op))
-#    data = _sch(op.data, l1, l2, r1, r2)
-#    data = reshape(op.data, map(dimension, (l1, l2, r1, r2)))
-#    data = permutedims(data, (4, 2, 3, 1))
-#    data = reshape(data, map(dimension, (l2⊗r2, l1⊗r1)))
-    dl1, dl2, dr1, dr2 = map(dimension, (l1, l2, r1, r2))
-    @cast A[(ν,μ), (n,m)] |= op.data[(m,μ), (n,ν)] (m ∈ 1:dl1, μ ∈ 1:dl2, n ∈ 1:dr1, ν ∈ 1:dr2)
-    return Operator(basis_fn(r2, l2), basis_fn(r1, l1), A)
+    #data = _sch(op.data, l1, l2, r1, r2)
+    #dl1, dl2, dr1, dr2 = map(dimension, (l1, l2, r1, r2))
+    #@cast A[(ν,μ), (n,m)] |= op.data[(m,μ), (n,ν)] (m ∈ 1:dl1, μ ∈ 1:dl2, n ∈ 1:dr1, ν ∈ 1:dr2)
+    data = reshape(op.data, map(dimension, (l1, l2, r1, r2)))
+    data = PermutedDimsArray(data, (4, 2, 3, 1))
+    data = reshape(data, map(dimension, (l2⊗r2, l1⊗r1)))
+    return Operator(basis_fn(r2, l2), basis_fn(r1, l1), data)
 end
 
 choi(op::SuperOperatorType) = _super_choi(ChoiBasis, op)

test/test_pauli.jl

@@ -2,41 +2,48 @@ using LinearAlgebra
 using Test
 
 using QuantumOpticsBase
+import QuantumOpticsBase: comp_pauli_kb
 
 @testset "pauli" begin
 
 b = SpinBasis(1//2)
+Isop = sprepost(identityoperator(b), dagger(identityoperator(b)))
+Xsop = sprepost(sigmax(b), dagger(sigmax(b)))
+Ysop = sprepost(sigmay(b), dagger(sigmay(b)))
+Zsop = sprepost(sigmaz(b), dagger(sigmaz(b)))
+Xsk = vec(sigmax(b))
+V = comp_pauli_kb(1)
+
+
 # Test conversion of unitary matrices to superoperators.
-q2 = b^2
-q3 = b^3
-CZ = DenseOperator(q2, q2, diagm(0 => [1,1,1,-1]))
-CZ_sop = SuperOperator(CZ)
+CZ = dm(spinup(b))⊗identityoperator(b) + dm(spindown(b))⊗sigmaz(b)
+CZ_sop = sprepost(CZ,dagger(CZ))
 
 # Test conversion of unitary matrices to superoperators.
 @test diag(CZ_sop.data) ==  ComplexF64[1,1,1,-1,1,1,1,-1,1,1,1,-1,-1,-1,-1,1]
-@test CZ_sop.basis_l == CZ_sop.basis_r == (q2, q2)
+@test basis_l(CZ_sop) == basis_r(CZ_sop) == KetBraBasis(b^2, b^2)
 
 # Test conversion of superoperator to Pauli transfer matrix.
-CZ_ptm = PauliTransferMatrix(CZ_sop)
+CZ_ptm = pauli(CZ_sop)
 
 # Test DensePauliTransferMatrix constructor.
-@test_throws DimensionMismatch DensePauliTransferMatrix((q2, q2), (q3, q3), CZ_ptm.data)
-@test DensePauliTransferMatrix((q2, q2), (q2, q2), CZ_ptm.data) == CZ_ptm
+@test_throws DimensionMismatch Operator(PauliBasis(2), PauliBasis(3), CZ_ptm.data)
+@test Operator(PauliBasis(2), PauliBasis(2), CZ_ptm.data) == CZ_ptm
 
-@test all(isapprox.(CZ_ptm.data[[1,30,47,52,72,91,117,140,166,185,205,210,227,256]], 1))
-@test all(isapprox.(CZ_ptm.data[[106,151]], -1))
+@test all(isapprox.(CZ_ptm.data[[1,30,47,52,72,91,117,140,166,185,205,210,227,256]], 4))
+@test all(isapprox.(CZ_ptm.data[[106,151]], -4))
 
 @test CZ_ptm == PauliTransferMatrix(ChiMatrix(CZ))
 
 # Test construction of non-symmetric unitary.
-CNOT = DenseOperator(q2, q2, diagm(0 => [1,1,0,0], 1 => [0,0,1], -1 => [0,0,1]))
+CNOT = DenseOperator(b^2, b^2, diagm(0 => [1,1,0,0], 1 => [0,0,1], -1 => [0,0,1]))
 CNOT_sop = SuperOperator(CNOT)
 CNOT_chi = ChiMatrix(CNOT)
 CNOT_ptm = PauliTransferMatrix(CNOT)
 
-@test CNOT_sop.basis_l == CNOT_sop.basis_r == (q2, q2)
-@test CNOT_chi.basis_l == CNOT_chi.basis_r == (q2, q2)
-@test CNOT_ptm.basis_l == CNOT_ptm.basis_r == (q2, q2)
+@test CNOT_sop.basis_l == CNOT_sop.basis_r == (b^2, b^2)
+@test CNOT_chi.basis_l == CNOT_chi.basis_r == (b^2, b^2)
+@test CNOT_ptm.basis_l == CNOT_ptm.basis_r == (b^2, b^2)
 
 @test all(isapprox.(imag.(CNOT_sop.data), 0))
 @test all(isapprox.(imag.(CNOT_chi.data), 0))
@@ -49,13 +56,13 @@ CNOT_ptm = PauliTransferMatrix(CNOT)
 @test all(isapprox.(CNOT_ptm.data[[123,168]], -1))
 
 # Test DenseChiMatrix constructor.
-@test_throws DimensionMismatch DenseChiMatrix((q2, q2), (q3, q3), CNOT_chi.data)
-@test DenseChiMatrix((q2, q2), (q2, q2), CNOT_chi.data) == CNOT_chi
+@test_throws DimensionMismatch Operator(ChoiBasis(PauliBasis(2), PauliBasis(2)), ChoiBasis(PauliBasis(3), PauliBasis(3)), CNOT_chi.data)
+@test Operator(ChoiBasis(PauliBasis(2), PauliBasis(2)), CNOT_chi.data) == CNOT_chi
 
 # Test equality and conversion among all three bases.
 ident = Complex{Float64}[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
 
-IDENT = DenseOperator(q2, ident)
+IDENT = DenseOperator(b^2, ident)
 
 IDENT_sop = SuperOperator(IDENT)
 IDENT_chi = ChiMatrix(IDENT)
@@ -71,7 +78,7 @@ IDENT_ptm = PauliTransferMatrix(IDENT)
 # Test approximate equality and conversion among all three bases.
 cphase = Complex{Float64}[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 exp(1im*.6)]
 
-CPHASE = DenseOperator(q2, cphase)
+CPHASE = DenseOperator(b^2, cphase)
 
 CPHASE_sop = SuperOperator(CPHASE)
 CPHASE_chi = ChiMatrix(CPHASE)