[non-Clifford] In-place Pauli measurements (projectrand!) for GeneralizedStabilizer with expanded test suite

docs/src/references.bib

@@ -191,6 +191,16 @@ @article{nahum2017quantum
   year = {2017}
 }
 
+% non-Clifford formalism
+
+@article{yoder2012generalization,
+  title={A generalization of the stabilizer formalism for simulating arbitrary quantum circuits},
+  author={Yoder, Theodore J},
+  journal={See http://www. scottaaronson. com/showcase2/report/ted-yoder. pdf},
+  year={2012},
+  publisher={Citeseer}
+}
+
 % codes
 
 @article{mackay2004sparse,

src/nonclifford.jl

@@ -177,6 +177,18 @@ function _allthreesumtozero(a,b,c)
     true
 end
 
+"""Creates a binary vector `k` of the same length as the input vector `b`, with exactly one nonzero element.
+The single nonzero element in `k` is positioned at the index of the first nonzero element in `b`. Example:
+If `k = BitVector([0, 1, 0, 0, 1])`, then `_create_k(k)` returns `BitVector([0, 1, 0, 0, 0])`."""
+function _create_k(b::BitVector)
+    k = falses(length(b))
+    pos = findfirst(b)
+    if pos !== nothing
+        k[pos] = true
+    end
+    return k
+end
+
 """$(TYPEDSIGNATURES)
 
 Performs a randomized projection of the state represented by the [`GeneralizedStabilizer`](@ref) `sm`,
@@ -224,16 +236,66 @@ julia> prob₁ = (real(χ′)+1)/2
 
 See also: [`expect`](@ref)
 """
-function projectrand!(sm::GeneralizedStabilizer, p::PauliOperator)
-    χ′ = expect(p, sm)
-    # Compute the probability of measuring in the +1 eigenstate
-    prob₁ = (real(χ′)+1)/2
-    # Randomly choose projection based on this probability
-    return _proj(sm, rand() < prob₁ ? p : -p)
+function _projectrand_notnorm(sm::GeneralizedStabilizer, p::PauliOperator)
+    # Returns the updated `GeneralizedStabilizer` state sm′ = (χ′, B(S′, D′)),
+    # where (S′, D′) is derived from (S, D) through the traditional stabilizer update,
+    # and χ′ is the updated density matrix after measurement. Note: Λ(χ′) ≤ Λ(χ).
+    dict = sm.destabweights
+    dtype = valtype(dict)
+    tzero = zero(dtype)
+    tone = one(dtype)
+    stab = sm.stab
+    newdict = typeof(dict)(tzero)
+    phase, b, c = rowdecompose(p, stab)
+
+    # Implementation of the in-place Pauli measurement quantum operation (Algorithm 2)
+    # on a generalized stabilizer by Ted Yoder (Page 8) from [yoder2012generalization](@cite).
+    if all(x -> x == 0, b)
+        # (Eq. 14-17)
+        for ((dᵢ, dⱼ), χ) in dict
+            if (im^phase * (-tone)^(dot(dᵢ, c)) == 1) && (im^phase * (-tone)^(dot(dⱼ, c)) == 1) # (Eq. 16)
+                newdict[(dᵢ,dⱼ)] += χ
+            end
+        end
+        sm.destabweights = newdict # in-place
+        state, res = projectrand!(stab, p)
+        sm.stab = state # in-place
+        return sm, res # the stabilizer basis (S, D) is not updated (Eq. 17)
+    else
+        # (Eq. 18-26)
+        k = _create_k(b)
+        for ((dᵢ, dⱼ), χ) in dict
+            x, y = dᵢ, dⱼ
+            q = 1
+            if dot(dᵢ, k) == 1
+                q *= im^phase * (-tone)^dot(dᵢ, c)
+                x = dᵢ .⊻ b
+            end
+            if dot(dⱼ, k) == 1
+                q *= conj(im^(phase)) * (-tone)^dot(dⱼ, c) # α* is conj(im^(phase))
+                y = dⱼ .⊻ b
+            end
+            χ′ = 1/2 * χ * q
+            newdict[(x,y)] += χ′
+        end
+        sm.destabweights = newdict # in-place
+        state, res = projectrand!(stab, p) # traditional stabilizer update
+        sm.stab = state # in-place
+        return sm, res
+    end
 end
 
-function _proj(sm::GeneralizedStabilizer, p::PauliOperator)
-    error("This functionality is not implemented yet")
+function projectrand!(sm::GeneralizedStabilizer, p::PauliOperator)
+    sm, res = _projectrand_notnorm(sm, p)
+    dict = sm.destabweights
+    if sm |> QuantumClifford.invsparsity == 1
+        for ((dᵢ, dⱼ), χ) in dict
+            χ′ = χ/LinearAlgebra.tr(χ) # normalization
+            sm.destabweights[(dᵢ, dⱼ)] = χ′
+        end
+        return sm, res
+    end
+    return sm, res
 end
 
 function project!(s::GeneralizedStabilizer, p::PauliOperator)

test/test_jet.jl

@@ -40,5 +40,5 @@ end
     )
     @show rep
     @test_broken length(JET.get_reports(rep)) == 0
-    @test length(JET.get_reports(rep)) <= 24
+    @test length(JET.get_reports(rep)) <= 30
 end

test/test_nonclifford.jl

@@ -54,6 +54,67 @@ end
     end
 end
 
+##
+
+@testset "PauliChannel: inverse sparsity of k-qubit channels" begin
+    # In general, the increase in Λ(χ) due to k-qubit channels is limited to a maximum factor of 16ᵏ.
+    k = 10
+    for n in 1:10
+        for repetition in 1:1000
+            stab = random_stabilizer(n)
+            pauli = random_pauli(n)
+            genstab = GeneralizedStabilizer(stab)
+            Λ_χ = genstab |> invsparsity # Λ(χ)
+            i = rand(1:n)
+            nc = embed(n, i, pcT)
+            for i in 1:k
+                apply!(genstab, nc) # in-place
+            end
+            Λ_χ′ = genstab |> invsparsity # Λ(χ′)
+            @test (Λ_χ′/Λ_χ) <= 16^k
+        end
+    end
+end
+
+##
+
+@testset "PauliChannel: Λ(χ) ≤ Λ(χ′) ≤ Λ(E)Λ(χ)" begin
+    # In general, the complexity of χ increases in the case of evolution through a channel.
+    for n in 1:10
+        for repetition in 1:1000
+            stab = random_stabilizer(n)
+            pauli = random_pauli(n)
+            genstab = GeneralizedStabilizer(stab)
+            Λ_χ = genstab |> invsparsity # Λ(χ)
+            i = rand(1:n)
+            nc = embed(n, i, pcT)
+            Λ_E = nc |> invsparsity # Λ(E)
+            apply!(genstab, nc) # in-place
+            Λ_χ′ = genstab |> invsparsity # Λ(χ′)
+            @test Λ_χ <= Λ_χ′ <= Λ_E*Λ_χ # Corollary 14, page 9 of [yoder2012generalization](@cite).
+        end
+    end
+end
+
+##
+
+@testset "PauliMeasurement: Λ(χ′) ≤ Λ(χ)" begin
+    # In general, the complexity of χ may decrease in the case of evolution through a measurement.
+    for n in 1:10
+        for repetition in 1:1000
+            stab = random_stabilizer(n)
+            pauli = random_pauli(n)
+            genstab = GeneralizedStabilizer(stab)
+            Λ_χ = genstab |> invsparsity # Λ(χ)
+            projectrand!(genstab, pauli)[1] # in-place
+            Λ_χ′ = genstab |> invsparsity # Λ(χ′)
+            @test Λ_χ′ <= Λ_χ # Corollary 14, page 9 of [yoder2012generalization](@cite).
+        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,5 +1,5 @@
 using QuantumClifford
-using QuantumClifford: GeneralizedStabilizer, rowdecompose, PauliChannel, mul_left!, mul_right!
+using QuantumClifford: GeneralizedStabilizer, rowdecompose, PauliChannel, mul_left!, mul_right!, invsparsity, _projectrand_notnorm
 using QuantumClifford: @S_str, random_stabilizer
 using QuantumOpticsBase
 using LinearAlgebra
@@ -86,3 +86,107 @@ end
         @test copy(nc) == nc
     end
 end
+
+function _projrand(τ,p)
+    qo_state = Operator(τ)
+    projectrand!(τ, p)[1] # in-place
+    qo_state_after_proj = Operator(τ)
+    qo_pauli = Operator(p)
+    qo_proj1 = (identityoperator(qo_pauli) - qo_pauli)/2
+    qo_proj2 = (identityoperator(qo_pauli) + qo_pauli)/2
+    result1 = qo_proj1*qo_state*qo_proj1'
+    result2 = qo_proj2*qo_state*qo_proj2'
+    # Normalize to ensure consistent comparison of the projected state
+    norm_qo_state_after_proj = iszero(qo_state_after_proj) ? qo_state_after_proj : qo_state_after_proj/tr(qo_state_after_proj)
+    norm_result1 = iszero(result1) ? result1 : result1/tr(result1)
+    norm_result2 = iszero(result2) ? result2 : result2/tr(result2)
+    return norm_qo_state_after_proj, norm_result1, norm_result2
+end
+
+@testset "Single-qubit projections using for stabilizer states" begin
+    for s in [S"X", S"Y", S"Z", S"-X", S"-Y", S"-Z"]
+        for p in [P"X", P"Y", P"Z", P"-X", P"-Y", P"-Z"]
+            genstab = GeneralizedStabilizer(s)
+            apply!(genstab, pcT) # in-place
+            norm_qo_state_after_proj, norm_result1, norm_result2 = _projrand(genstab,p) # in-place
+            !(iszero(norm_qo_state_after_proj)) && @test real(tr(norm_qo_state_after_proj)) ≈ 1
+            @test projectrand!(genstab, p)[1] |> invsparsity <= genstab |> invsparsity # Λ(χ′) ≤ Λ(χ)
+            @test norm_qo_state_after_proj ≈ norm_result2 || norm_qo_state_after_proj ≈ norm_result1
+       end
+    end
+
+    for _ in 1:100
+        for n in 1:1
+            stab = random_stabilizer(n)
+            genstab = GeneralizedStabilizer(stab)
+            pauli = random_pauli(n)
+            apply!(genstab, pcT) # in-place
+            norm_qo_state_after_proj, norm_result1, norm_result2 = _projrand(genstab,pauli) # in-place
+            !(iszero(norm_qo_state_after_proj)) && @test real(tr(norm_qo_state_after_proj)) ≈ 1
+            @test projectrand!(genstab, pauli)[1] |> invsparsity <= genstab |> invsparsity # Λ(χ′) ≤ Λ(χ)
+            @test norm_qo_state_after_proj ≈ norm_result2 || norm_qo_state_after_proj ≈ norm_result1
+        end
+    end
+end
+
+@testset "Multi-qubit projections using GeneralizedStabilizer for stabilizer states" begin
+    for n in 1:5
+        for repetition in 1:3
+            for state in [Stabilizer, MixedDestabilizer, GeneralizedStabilizer]
+                s = random_stabilizer(n)
+                p = random_pauli(n)
+                τ = state(s)
+                apply!(τ, random_clifford(n)) # in-place
+                norm_qo_state_after_proj, norm_result1, norm_result2 = _projrand(τ,p) # in-place
+                !(iszero(norm_qo_state_after_proj)) && @test real(tr(norm_qo_state_after_proj)) ≈ 1
+                @test norm_qo_state_after_proj ≈ norm_result2 || norm_qo_state_after_proj ≈ norm_result1
+                isa(τ, GeneralizedStabilizer) && @test projectrand!(τ, p)[1] |> invsparsity <= τ |> invsparsity # Λ(χ′) ≤ Λ(χ)
+            end
+        end
+    end
+end
+
+@testset "The trace Tr[χ′] is the probability of measuring an outcome" begin
+    # The trace Tr[χ′] represents the probability of obtaining an outcome of 0.
+    # TODO Document - Since ((real(expect(P"Z", apply!(genstab(S"-Z"), pcT))))+1)/2 is 0,
+    # it triggers Eq. 16, where (I+(-1)^(i·c)*M)ρₛ(I+(-1)^(j·c)*M) evaluates to 0. Thus,
+    # genstab after projectrand!(apply!(genstab(S"-Z"), pcT), P"Z")[1] has no meaning.
+    for s in [S"X", S"Y", S"Z", S"-X", S"-Y", S"Z"]
+        for p in [P"X", P"Y", P"Z", P"-X", P"-Y"]
+            genstab = GeneralizedStabilizer(s)
+            apply!(genstab, pcT) # in-place
+            prob1 = (real(expect(p, genstab))+1)/2
+            _projectrand_notnorm(genstab, p)[1] # in-place
+            dict = genstab.destabweights
+            trace_χ′ = real(tr(collect(values(dict))[1])) # Tr[χ′]
+            @test isapprox(prob1, trace_χ′; atol=1e-5)
+        end
+    end
+end
+
+@testset "Multi-qubit projections using GeneralizedStabilizer with multiple non-Clifford gates" begin
+    # TODO Analyze some multi-qubit genstab states that are unsimulable due to very complex
+    # destabweights, which also exhibit an inverse sparsity relation (Λ(χ′) = Λ(χ) = 4).
+    count = 0
+    num_trials = 20
+    num_qubits = [2,3,4,5] # exclusively multi-qubit
+    for n in num_qubits # exponential cost in this term
+        for repetition in 1:num_trials
+            stab = random_stabilizer(n)
+            pauli = random_pauli(n)
+            genstab = GeneralizedStabilizer(stab)
+            # some (repeated) non-clifford ops
+            i = rand(1:n)
+            nc = embed(n, i, pcT)
+            apply!(genstab, nc) # in-place
+            apply!(genstab, nc) # in-place
+            apply!(genstab, nc) # in-place
+            norm_qo_state_after_proj, norm_result1, norm_result2 = _projrand(genstab,pauli) # in-place
+            !(iszero(norm_qo_state_after_proj)) && @test real(tr(norm_qo_state_after_proj)) ≈ 1
+            @test projectrand!(genstab, pauli)[1] |> invsparsity <= genstab |> invsparsity # Λ(χ′) ≤ Λ(χ)
+            count += (norm_qo_state_after_proj ≈ norm_result2 || norm_qo_state_after_proj ≈ norm_result1) ? 1 : 0
+        end
+    end
+    prob = count/(num_trials*(length(num_qubits)))
+    @test prob > 0.7
+end