Even more WIP Pauli

Author: akirakyle

src/QuantumOpticsBase.jl

@@ -39,7 +39,7 @@ export Basis, GenericBasis, CompositeBasis, basis, basis_l, basis_r, dimension,
                 AbstractTimeDependentOperator, TimeDependentSum, set_time!,
                 current_time, time_shift, time_stretch, time_restrict, static_operator,
         #superoperators
-                KetBraBasis, ChoiBasis, PauliBasis,
+                KetBraBasis, ChoiBasis, PauliBasis, ChiBasis,
                 vec, unvec, super, choi, kraus, stinespring, pauli, chi,
                 spre, spost, sprepost, liouvillian, identitysuperoperator,
                 SuperOperatorType, DenseSuperOpType, SparseSuperOpType,

src/pauli.jl

@@ -1,14 +1,14 @@
 
 # comp stands for computational basis
-function _comp_pauli_1(basis_fn)
+function _pauli_comp_1(b_out_fn, b_in)
     b = SpinBasis(1//2)
-    vec_it(fn) = vec(fn(b).data)/2 # provides "standard" normalization
+    vec_it(fn) = vec((fn(b)/sqrt(2)).data) # provides "standard" normalization
     V = hcat(map(vec_it, [identityoperator, sigmax, sigmay, sigmaz])...)
-    Operator(basis_fn(b, b), PauliBasis(1), V)
+    Operator(b_out_fn(b, b), b_in, V)
 end
 
-_comp_pauli_kb1_cached = _comp_pauli_1(KetBraBasis)
-_comp_pauli_choi1_cached = _comp_pauli_1(ChoiBasis)
+_pauli_comp_kb1_cached = _pauli_comp_1(KetBraBasis, PauliBasis(1))
+_pauli_comp_choi1_cached = _pauli_comp_1(ChoiBasis, ChiBasis(2))
 
 # TODO: should this be further cached?
 """
@@ -17,52 +17,43 @@ _comp_pauli_choi1_cached = _comp_pauli_1(ChoiBasis)
 Creates a superoperator which changes from the computational `KetBra(SpinBasis(1//2))`
 to the `PauliBasis()` over `N` qubits.
 """
-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
+pauli_comp_kb(N::Integer) = tensor_pow(_pauli_comp_kb1_cached, N)
+pauli_comp_choi(N::Integer) = tensor_pow(_pauli_comp_choi1_cached, N)
 
 # 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
+pauli(op::Operator) = dagger(pauli_comp_kb(length(basis(op)))) * vec(op)
 
 function pauli(op::SOpKetBraType)
-    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(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...
+    Nl, Nr = length(basis_l(basis_l(op))), length(basis_l(basis_r(op)))
+    Vl = pauli_comp_kb(Nl)
+    Vr = Nl == Nr ? Vl : pauli_comp_kb(Nr)
+    dagger(Vl) * op * Vr
 end
 
 function chi(op::ChoiStateType)
-    (basis_l(op) == basis_r(op)) || throw(ArgumentError("Choi state must map between square operators in order to be converted to chi represenation"))
-    bl, br = basis_l(basis_l(op)), basis_r(basis_l(op))
-    foreach(_check_is_spinbasis, (bl, br))
+    Nl, Nr = length(basis_l(basis_l(op))), length(basis_r(basis_l(op)))
+    Vl = pauli_comp_choi(Nl)
+    Vr = Nl == Nr ? Vl : pauli_comp_choi(Nr)
+    dagger(Vl) * op * Vr
+end
 
+function _pauli_chi(basis_fn, op)
+    bl, br = basis_l(op), basis_r(op)
     Nl, Nr = length(bl), length(br)
-    Vl = comp_pauli_kb(Nl)
-    Vr = Nl == Nr ? Vl : comp_pauli_kb(Nr)
-    dagger(Vl) * op * Vr # TODO figure out normalization
+    data = reshape(op.data, (2^Nl, 2^Nl, 2^Nr, 2^Nr))
+    data = PermutedDimsArray(data, (4, 2, 3, 1))
+    data = reshape(data, (2^(Nl+Nr), 2^(Nl+Nr)))
+    return Operator(basis_fn(Nl+Nr), basis_fn(Nl+Nr), data)
 end
 
+pauli(op::ChiType) = _pauli_chi(PauliBasis, op)
+chi(op::SOpPauliType) = _pauli_chi(ChiBasis, op)
+pauli(op::ChoiStateType) = pauli(chi(op))
+chi(op::SOpKetBraType) = chi(pauli(op))
+
 """
 function hwpauli(op::SuperOperatorType; tol=1e-9)
     bl, br = basis_l(op), basis_r(op)

src/superoperators.jl

@@ -1,4 +1,4 @@
-import QuantumInterface: KetBraBasis, ChoiBasis, PauliBasis
+import QuantumInterface: KetBraBasis, ChoiBasis, PauliBasis, ChiBasis
 using TensorCast
 
 const SOpBasis = Union{KetBraBasis, PauliBasis}
@@ -83,7 +83,7 @@ function _super_choi(basis_fn, op)
 end
 
 choi(op::SuperOperatorType) = _super_choi(ChoiBasis, op)
-super(op::ChoiStateType) =  _super_choi(KetBraBasis, op)
+super(op::ChoiStateType) = _super_choi(KetBraBasis, op)
 
 # I'm not sure this is actually right... see sec V.C. of https://arxiv.org/abs/1111.6950
 dagger(a::ChoiStateType) = choi(dagger(super(a)))

test/test_pauli.jl

@@ -2,48 +2,68 @@ 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)
-
+bs = SpinBasis(1//2)
+bp = PauliBasis(1)
+
+# Test Pauli basis vectors are in I, X, Y, Z order
+@test pauli(identityoperator(bs)/sqrt(2)).data ≈ [1., 0, 0, 0]
+@test pauli(sigmax(bs)/sqrt(2)).data ≈ [0, 1., 0, 0]
+@test pauli(sigmay(bs)/sqrt(2)).data ≈ [0, 0, 1., 0]
+@test pauli(sigmaz(bs)/sqrt(2)).data ≈ [0, 0, 0, 1.]
+
+# Test that single qubit unitary Pauli channels are diagonal
+Isop = sprepost(identityoperator(bs), dagger(identityoperator(bs)))
+Xsop = sprepost(sigmax(bs), dagger(sigmax(bs)))
+Ysop = sprepost(sigmay(bs), dagger(sigmay(bs)))
+Zsop = sprepost(sigmaz(bs), dagger(sigmaz(bs)))
+@test pauli(Isop).data ≈ diagm([1, 1, 1, 1])
+@test pauli(Xsop).data ≈ diagm([1, 1, -1, -1])
+@test pauli(Ysop).data ≈ diagm([1, -1, 1, -1])
+@test pauli(Zsop).data ≈ diagm([1, -1, -1, 1])
+
+# Test bit flip encoder isometry
+encoder_kraus = (tensor_pow(spinup(bs), 3) ⊗ dagger(spinup(bs)) +
+                 tensor_pow(spindown(bs), 3) ⊗ dagger(spindown(bs)))
+encoder_sup = sprepost(encoder_kraus, dagger(encoder_kraus))
+decoder_sup = sprepost(dagger(encoder_kraus), encoder_kraus)
+@test super(choi(encoder_sup)).data == encoder_sup.data
+@test decoder_sup == dagger(encoder_sup)
+@test choi(decoder_sup) == dagger(choi(encoder_sup))
+@test decoder_sup*encoder_sup ≈ dense(identitysuperoperator(bs))
+@test decoder_sup*choi(encoder_sup) ≈ dense(identitysuperoperator(bs))
+@test choi(decoder_sup)*encoder_sup ≈ dense(identitysuperoperator(bs))
+@test super(choi(decoder_sup)*choi(encoder_sup)) ≈ dense(identitysuperoperator(bs))
 
 # Test conversion of unitary matrices to superoperators.
-CZ = dm(spinup(b))⊗identityoperator(b) + dm(spindown(b))⊗sigmaz(b)
+CZ = dm(spinup(bs))⊗identityoperator(bs) + dm(spindown(bs))⊗sigmaz(bs)
 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 basis_l(CZ_sop) == basis_r(CZ_sop) == KetBraBasis(b^2, b^2)
+@test CZ_sop.data ≈ diagm([1,1,1,-1,1,1,1,-1,1,1,1,-1,-1,-1,-1,1])
+@test basis_l(CZ_sop) == basis_r(CZ_sop) == KetBraBasis(bs^2, bs^2)
 
 # Test conversion of superoperator to Pauli transfer matrix.
 CZ_ptm = pauli(CZ_sop)
 
-# Test DensePauliTransferMatrix constructor.
+# Test construction of dense Pauli transfer matrix
 @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]], 4))
-@test all(isapprox.(CZ_ptm.data[[106,151]], -4))
+@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 CZ_ptm == PauliTransferMatrix(ChiMatrix(CZ))
+@test CZ_ptm == pauli(chi(CZ_sop))
 
-# Test construction of non-symmetric unitary.
-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)
+CNOT = dm(spinup(bs))⊗identityoperator(bs) + dm(spindown(bs))⊗sigmax(bs)
+CNOT_sop = sprepost(CZ,dagger(CZ))
+CNOT_chi = chi(CNOT_sop)
+CNOT_ptm = pauli(CNOT_sop)
 
-@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 basis_l(CNOT_sop) == basis_r(CNOT_sop) == KetBraBasis(bs^2, bs^2)
+@test basis_l(CNOT_chi) == basis_r(CNOT_chi) == ChoiBasis(bp^2, bp^2)
+@test basis_l(CNOT_ptm) == basis_r(CNOT_ptm) == ChoiBasis(bp^2, bp^2)
 
 @test all(isapprox.(imag.(CNOT_sop.data), 0))
 @test all(isapprox.(imag.(CNOT_chi.data), 0))
@@ -55,25 +75,25 @@ CNOT_ptm = PauliTransferMatrix(CNOT)
 @test all(isapprox.(CNOT_ptm.data[[1,18,47,64,70,85,108,138,153,183,205,222,227,244,]], 1))
 @test all(isapprox.(CNOT_ptm.data[[123,168]], -1))
 
-# Test DenseChiMatrix constructor.
-@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 construction of chi matrix
+@test_throws DimensionMismatch Operator(ChoiBasis(bp^2, bp^2), ChoiBasis(bp^3, bp^3), CNOT_chi.data)
+@test Operator(ChoiBasis(bp^2, bp^2), CNOT_chi.data) == CNOT_chi
 
-# Test equality and conversion among all three bases.
+# Test equality and conversion of identity 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(b^2, ident)
+IDENT = DenseOperator(bs^2, ident)
 
 IDENT_sop = SuperOperator(IDENT)
 IDENT_chi = ChiMatrix(IDENT)
 IDENT_ptm = PauliTransferMatrix(IDENT)
 
-@test ChiMatrix(IDENT_sop) == IDENT_chi
-@test ChiMatrix(IDENT_ptm) == IDENT_chi
-@test SuperOperator(IDENT_chi) == IDENT_sop
-@test SuperOperator(IDENT_ptm) == IDENT_sop
-@test PauliTransferMatrix(IDENT_sop) == IDENT_ptm
-@test PauliTransferMatrix(IDENT_chi) == IDENT_ptm
+@test chi(IDENT_sop) == IDENT_chi
+@test chi(IDENT_ptm) == IDENT_chi
+@test super(IDENT_chi) == IDENT_sop
+@test super(IDENT_ptm) == IDENT_sop
+@test pauli(IDENT_sop) == IDENT_ptm
+@test pauli(IDENT_chi) == IDENT_ptm
 
 # 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)]

test/test_superoperators.jl

@@ -224,7 +224,7 @@ b_fock = FockBasis(5)
 z_l = normalize(fockstate(b_fock, 0) + fockstate(b_fock, 4))
 o_l = fockstate(b_fock, 2)
 encoder_kraus = z_l ⊗ dagger(spinup(b_logical)) + o_l ⊗ dagger(spindown(b_logical))
-    encoder_sup = sprepost(encoder_kraus, dagger(encoder_kraus))
+encoder_sup = sprepost(encoder_kraus, dagger(encoder_kraus))
 decoder_sup = sprepost(dagger(encoder_kraus), encoder_kraus)
 @test super(choi(encoder_sup)).data == encoder_sup.data
 @test decoder_sup == dagger(encoder_sup)