tensor product of GeneralizedStabilizer (#345)

Co-authored-by: Stefan Krastanov <[email protected]>
Co-authored-by: Stefan Krastanov <[email protected]>

Author: Fe-r-oz

docs/src/references.bib

@@ -199,11 +199,13 @@ @article{nahum2017quantum
 
 % non-Clifford formalism
 
-@inproceedings{Yoder2012AGO,
-  title={A generalization of the stabilizer formalism for simulating arbitrary quantum circuits},
-  author={Theodore J. Yoder},
-  year={2012},
-  url={https://api.semanticscholar.org/CorpusID:18406560}
+@misc{Yoder2012AGO,
+  title = {A Generalization of the Stabilizer Formalism for Simulating Arbitrary Quantum Circuits},
+  author = {Theodore J. Yoder},
+  year = {2012},
+  howpublished = {Semantic Scholar Corpus ID: 18406560},
+  url = {https://api.semanticscholar.org/CorpusID:18406560},
+  note = {Also available at \url{http://www.scottaaronson.com/showcase2/report/ted-yoder.pdf}}
 }
 
 % codes

src/nonclifford.jl

@@ -443,6 +443,121 @@ function apply!(state::GeneralizedStabilizer, gate::AbstractPauliChannel; prune_
     state
 end
 
+"""
+$(TYPEDSIGNATURES)
+
+Tensor product of [`GeneralizedStabilizer`](@ref) states.
+
+# Stabilizer state
+
+```jldoctest
+julia> sm = GeneralizedStabilizer(S"-X")
+A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
+𝒟ℯ𝓈𝓉𝒶𝒷
++ Z
+𝒮𝓉𝒶𝒷
+- X
+with ϕᵢⱼ | Pᵢ | Pⱼ:
+ 1.0+0.0im | + _ | + _
+
+julia> sm ⊗ sm
+A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
+𝒟ℯ𝓈𝓉𝒶𝒷
++ Z_
++ _Z
+𝒮𝓉𝒶𝒷
+- X_
+- _X
+with ϕᵢⱼ | Pᵢ | Pⱼ:
+ 1.0+0.0im | + __ | + __
+```
+
+# Arbitrary state
+
+```jldoctest
+julia> using LinearAlgebra; # hide
+
+julia> sm = GeneralizedStabilizer(ghz(2))
+A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
+𝒟ℯ𝓈𝓉𝒶𝒷
++ Z_
++ _X
+𝒮𝓉𝒶𝒷
++ XX
++ ZZ
+with ϕᵢⱼ | Pᵢ | Pⱼ:
+ 1.0+0.0im | + __ | + __
+
+julia> apply!(sm, embed(2, 2, pcT))
+A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
+𝒟ℯ𝓈𝓉𝒶𝒷
++ Z_
++ _X
+𝒮𝓉𝒶𝒷
++ XX
++ ZZ
+with ϕᵢⱼ | Pᵢ | Pⱼ:
+ 0.853553+0.0im | + __ | + __
+ 0.0+0.353553im | + __ | + Z_
+ 0.0-0.353553im | + Z_ | + __
+ 0.146447+0.0im | + Z_ | + Z_
+
+julia> newsm = sm ⊗ sm
+A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
+𝒟ℯ𝓈𝓉𝒶𝒷
++ Z___
++ _X__
++ __Z_
++ ___X
+𝒮𝓉𝒶𝒷━━
++ XX__
++ ZZ__
++ __XX
++ __ZZ
+with ϕᵢⱼ | Pᵢ | Pⱼ:
+ 0.0-0.301777im | + Z___ | + ____
+ -0.125+0.0im | + Z_Z_ | + ____
+ 0.125+0.0im | + Z___ | + Z___
+ 0.728553+0.0im | + ____ | + ____
+ 0.0-0.0517767im | + Z_Z_ | + Z___
+ 0.0-0.301777im | + __Z_ | + ____
+ 0.0+0.301777im | + ____ | + Z___
+ 0.125+0.0im | + __Z_ | + Z___
+ 0.125+0.0im | + Z___ | + __Z_
+ 0.0-0.0517767im | + Z_Z_ | + __Z_
+ 0.0+0.0517767im | + Z___ | + Z_Z_
+ 0.0+0.301777im | + ____ | + __Z_
+ 0.0214466+0.0im | + Z_Z_ | + Z_Z_
+ 0.125+0.0im | + __Z_ | + __Z_
+ -0.125+0.0im | + ____ | + Z_Z_
+ 0.0+0.0517767im | + __Z_ | + Z_Z_
+
+julia> real(tr(newsm))
+1.0
+```
+"""
+function (⊗)(state₁::GeneralizedStabilizer, state₂::GeneralizedStabilizer)
+    dict₁ = state₁.destabweights
+    dict₂ = state₂.destabweights
+    dtype = valtype(dict₁)
+    tzero = zero(dtype)
+    newdict = DefaultDict{Tuple{BitVector,BitVector},dtype}(tzero)
+    newstab = state₁.stab ⊗ state₂.stab
+    for ((d1_i, d1_j), χ) in dict₁ # χ = ϕᵢⱼ for state₁
+        for ((d2_i, d2_j), χ′) in dict₂ # χ′ = ϕₖₗ for state₂
+            # Combine the Pauli operators via tensor product: Pᵢ ⊗ Pₖ
+            # and Pⱼ ⊗ Pₗ. vcat implements P₁ ⊗ P₂ as bitwise concatenation.
+            new_key_i = vcat(d1_i, d2_i)
+            new_key_j = vcat(d1_j, d2_j)
+            # The new coefficient is ϕᵢⱼ * ϕₖₗ because:
+            # ∑ϕᵢⱼ*ϕₖₗ(Pᵢρ₁Pⱼ†) ⊗ (Pₖρ₂Pₗ†) = ∑ϕᵢⱼ*ϕₖₗ(Pᵢ ⊗ Pₖ)(ρ₁ ⊗ ρ₂)(Pⱼ ⊗ Pₗ)† and
+            # thus the combined weight is the product of the individual weights.
+            newdict[(new_key_i, new_key_j)] += χ * χ′
+        end
+    end
+    return GeneralizedStabilizer(newstab, newdict)
+end
+
 """Decompose a Pauli ``P`` in terms of stabilizer and destabilizer rows from a given tableaux.
 
 For given tableaux of rows destabilizer rows ``\\{d_i\\}`` and stabilizer rows ``\\{s_i\\}``,

test/test_nonclifford.jl

@@ -1,6 +1,6 @@
 @testitem "NonClifford" begin
     using QuantumClifford
-    using QuantumClifford: GeneralizedStabilizer, rowdecompose, PauliChannel, invsparsity, mul_left!, mul_right!
+    using QuantumClifford: GeneralizedStabilizer, rowdecompose, PauliChannel, invsparsity, mul_left!, mul_right!, mixed_destab_looks_good, tr
     using Test
     using InteractiveUtils
     using Random
@@ -102,6 +102,27 @@
         end
     end
 
+    @testset "Tensor products over generalized stabilizers" begin
+        num_trials = 10
+        num_qubits = [2,3,4,5] # exclusively multi-qubit
+        for n in num_qubits
+            for repetition in 1:num_trials
+                stab = random_stabilizer(n)
+                pauli = random_pauli(n)
+                genstab = GeneralizedStabilizer(stab)
+                # Apply some (repeated) non-Clifford operations
+                i = rand(1:n)
+                nc = embed(n, i, pcT)
+                apply!(genstab, nc) # in-place
+                apply!(genstab, nc) # in-place
+                apply!(genstab, nc) # in-place
+                newsm = genstab ⊗ genstab
+                @test mixed_destab_looks_good(newsm.stab)
+                @test real(tr(newsm)) ≈ 1
+            end
+        end
+    end
+
     @test_throws ArgumentError GeneralizedStabilizer(S"XX")
     @test_throws ArgumentError PauliChannel(((P"X", P"Z"), (P"X", P"ZZ")), (1,2))
     @test_throws ArgumentError PauliChannel(((P"X", P"Z"), (P"X", P"Z")), (1,))

test/test_nonclifford_quantumoptics.jl

@@ -1,6 +1,6 @@
 @testitem "Non-Clifford Quantum Optics" begin
     using QuantumClifford
-    using QuantumClifford: GeneralizedStabilizer, rowdecompose, PauliChannel, mul_left!, mul_right!, invsparsity, _projectrand_notnorm
+    using QuantumClifford: GeneralizedStabilizer, rowdecompose, PauliChannel, mul_left!, mul_right!, invsparsity, _projectrand_notnorm, mixed_destab_looks_good, tr
     using QuantumClifford: @S_str, random_stabilizer
     using QuantumOpticsBase
     using LinearAlgebra
@@ -220,6 +220,35 @@
         @test any(projectrand!(GeneralizedStabilizer(S"X"), P"Z")[1].stab != s1 for _ in 1:10)
     end
 
+    @testset "Tensor products of generalized stabilizers" begin
+        num_trials = 3
+        num_qubits = [2,3] # exclusively multi-qubit
+        for n in num_qubits
+            for repetition in 1:num_trials
+                stab1 = random_stabilizer(n)
+                genstab1 = GeneralizedStabilizer(stab1)
+                stab2 = random_stabilizer(n)
+                genstab2 = GeneralizedStabilizer(stab2)
+                # apply some (repeated) non-Clifford operations to genstab1
+                for _ in 1:rand(1:5)
+                    i = rand(1:n)
+                    nc = embed(n, i, pcT)
+                    apply!(genstab1, nc)
+                end
+                # apply some (repeated) non-Clifford operations to genstab2
+                for _ in 1:rand(1:5)
+                    i = rand(1:n)
+                    nc = embed(n, i, pcT)
+                    apply!(genstab2, nc)
+                end
+                @test Operator(genstab1 ⊗ genstab2) ≈ Operator(genstab1) ⊗ Operator(genstab2)
+                @test Operator(genstab2 ⊗ genstab1) ≈ Operator(genstab2) ⊗ Operator(genstab1)
+                @test Operator(genstab1 ⊗ genstab1) ≈ Operator(genstab1) ⊗ Operator(genstab1)
+                @test Operator(genstab2 ⊗ genstab2) ≈ Operator(genstab2) ⊗ Operator(genstab2)
+            end
+        end
+    end
+
     @testset "smaller test redundant to the ones above" begin
         for n in 1:5
             for rep in 1:2

test/test_paulis.jl

@@ -1,6 +1,6 @@
 @testitem "Pauli Operators" begin
-  using QuantumClifford: apply_single_x!, apply_single_y!, apply_single_z!
-  test_sizes = [1,2,10,63,64,65,127,128,129] # Including sizes that would test off-by-one errors in the bit encoding.
+    using QuantumClifford: apply_single_x!, apply_single_y!, apply_single_z!
+    test_sizes = [1,2,10,63,64,65,127,128,129] # Including sizes that would test off-by-one errors in the bit encoding.
 
     @testset "Parsing, constructors, and properties" begin
         @test P"-iXYZ" == PauliOperator(0x3, 3, vcat(BitArray([1,1,0]).chunks, BitArray([0,1,1]).chunks))

test/test_stabs.jl

@@ -1,5 +1,6 @@
 @testitem "Stabilizers" begin
     using QuantumClifford
+    using LinearAlgebra
     using QuantumClifford: stab_looks_good, destab_looks_good, mixed_stab_looks_good, mixed_destab_looks_good
     test_sizes = [1,2,10,63,64,65,127,128,129] # Including sizes that would test off-by-one errors in the bit encoding.
     @testset "Pure and Mixed state initialization" begin