Implementation of a pauli observable structure

src/classical_shadows.jl

@@ -210,52 +210,29 @@ end
     observables_bit_strings(observables)   
 
 Return a set containing only the bitstrings of {0,1}ⁿ that contribute to the estimation of the expected value of
-at least one observable in the set `observables`. A bitstring contribute to the estimation of an observable's expected value
-if there is no 1 in the bitstring in the same position as a I in the observable.
+at least one observable in the set `observables`. For a bit string to contribute to the estimation of one observable, 
+it must contains 0 in the same position as the identity terms of this observable and 1 in the other positions.
 
 # Examples
 
 ```jldoctest
-julia> QuantumMeasurements.observables_bit_strings(Set(["XIII", "IIZY"]))
+julia> QuantumMeasurements.observables_bit_strings(Set([QuantumMeasurements.PauliObs(o) for o in ["XIII", "IIZY"]]))
 Set{Array{Int64}} with 5 elements:
-  [0, 0, 0, 1]
-  [0, 0, 0, 0]
   [0, 0, 1, 1]
-  [0, 0, 1, 0]
   [1, 0, 0, 0]
 ```
 """
-function observables_bit_strings(observables::Set{String})::Set{Array{Int}}
-    s = Set{Array{Int}}()
-    n = length(first(observables))
-    for o in observables #For each observable in observables we add to the set s the bitstrings that contribute to the estimation of its expected value.
-        k = weight(o)
-        bit_string = [zeros(Int, n) for _ = 1:(2^k)]
-
-        non_I = Array{Int}(undef, k) #Index of qubits where the observable o acts non trivally.
-        i = 1
-        for j = 1:n
-            if o[j] != 'I'
-                non_I[i] = j
-                i += 1
-            end
-        end
-        # At an index where the observable's factor is not I, the bitstring
-        # contribute to the inverse channel whether it is 0 or 1, so, after
-        # setting the bits in the same poistion as I to 0, we generate all the
-        # possibilities for the other positions.
-        bit = bitarray(collect(0:(2^k - 1)), k)
-        for z = 1:(2^k)
-            for j = 1:k
-                if bit[j, z]
-                    bit_string[z][non_I[j]] = 1
-                end
-            end
-        end
+function observables_bit_strings(observables::Set{PauliObs})::Set{Array{Int}}
+    return Set([observable_bit_string(o) for o in observables])
+end
 
-        push!(s, bit_string...)
-    end
-    return s
+function observable_bit_string(observable::PauliObs)::Array{Int}
+    n = observable.n_qubits
+    bit_string = zeros(Int, n)
+        for j in keys(observable.pauli)
+            bit_string[j] = 1
+        end
+    return bit_string
 end
 
 """

src/estimate.jl

@@ -108,8 +108,9 @@ function estimate_expect_shadows(
     robust::Bool=false,
     rng::AbstractRNG=GLOBAL_RNG,
 )::Dict{String,Float64}
-    M = length(observables)
-    k = max(weight.(observables)...)
+    pauli_observables = Set([PauliObs(o) for o in observables])
+    M = length(pauli_observables)
+    k = max(weight.(pauli_observables)...)
     n = nqubits(state)
 
     Ne, Ke, Nc, Kc = number_samples(k, precision, probability, M, n; robust)
@@ -125,14 +126,14 @@ function estimate_expect_shadows(
         noise_shadows_list = classical_shadows(zero_state(n), Nc * Kc; noise_model, rng)
         println("Noise shadows generated")
         println("Computing calibration coefficients")
-        bit_strings = observables_bit_strings(observables)
+        bit_strings = observables_bit_strings(pauli_observables)
         f = Calibration(noise_shadows_list, bit_strings, Nc, Kc)
         println("Calibration done")
 
         println("Estimating expectation values")
-        return robust_estimate_set_obs(shadows_list, f, observables, Ne, Ke)
+        return robust_estimate_set_obs(shadows_list, f, pauli_observables, Ne, Ke)
     else
         println("Estimating expectation values")
-        return estimate_set_obs(shadows_list, observables, Ne, Ke)
+        return estimate_set_obs(shadows_list, pauli_observables, Ne, Ke)
     end
 end

src/estimate_from_shadows.jl

@@ -26,9 +26,9 @@ function estimate_from_shadows(
     nb_shadows = length(shadows)
     Re > nb_shadows ? Re = nb_shadows : Re
     Ne = Re ÷ Ke
-
+    pauli_observables = Set([PauliObs(o) for o in observables])
     println("Using $(Ne*Ke) classical shadows split into $Ke equally-sized parts")
-    estimate_set_obs(shadows, observables, Ne, Ke)
+    estimate_set_obs(shadows, pauli_observables, Ne, Ke)
 end
 
 """
@@ -98,12 +98,13 @@ used to generate the array `shadows`. Uses a median of mean estimation using `Ke
 function estimate_from_shadows(
     shadows::Vector{Shadow}, calibration::Calibration, observables::Set{String}, Re::Int, Ke::Int=1
 )::Dict{String,Float64}
+    pauli_observables = Set([PauliObs(o) for o in observables])
     nb_shadows = length(shadows)
     Re > nb_shadows ? Re = nb_shadows : Re
     Ne = Re ÷ Ke
 
     println("Using $(Ne*Ke) classical shadows split into $Ke equally-sized parts")
-    robust_estimate_set_obs(shadows, calibration, observables, Ne, Ke)
+    robust_estimate_set_obs(shadows, calibration, pauli_observables, Ne, Ke)
 end
 
 """
@@ -307,7 +308,8 @@ the estimation of observables' expected values are computed.
 function calibrate(
     observables::Set{String}, noise_shadows::Array{Shadow}, Nc::Int, Kc::Int
 )::Calibration
-    bit_strings = observables_bit_strings(observables)
+    pauli_observables = Set([PauliObs(o) for o in observables])
+    bit_strings = observables_bit_strings(pauli_observables)
     calibrate(bit_strings, noise_shadows, Nc, Kc)
 end
 

src/observable_estimation.jl

@@ -49,12 +49,12 @@ Based on one `shadows`, give an estimation of the observable `obs`'s expected va
 which equals the expected value in average.
 Uses the kronecker factorization.
 """
-function single_estimate_obs(shadow::Shadow, obs::String)::Float64
+function single_estimate_obs(shadow::Shadow, obs::PauliObs)::Float64
     U = shadow.gate
     b = shadow.measure
     prod([
-        real(tr(Matrix(mat(pauli_from_char[obs[j]]))' * chain_inverse(Matrix(mat(U[j])), b[j]))) for
-        j = 1:(shadow.n_qubits)
+        real(tr(Matrix(mat(pauli_from_char[obs.pauli[j]]))' * chain_inverse(Matrix(mat(U[j])), b[j]))) for
+        j in keys(obs.pauli)
     ])
 end
 
@@ -64,7 +64,7 @@ end
 Based on `shadows` of a state, return an estimation of the observable `obs`'s expected value.
 Uses the median of mean protocol to have a good estimation.
 """
-function estimate_obs(shadows::Array{Shadow}, obs::String, Ne::Int, Ke::Int)::Float64
+function estimate_obs(shadows::Array{Shadow}, obs::PauliObs, Ne::Int, Ke::Int)::Float64
     R = Ne * Ke
     median_of_mean([single_estimate_obs(shadows[r], obs) for r = 1:R], Ne, Ke)
 end
@@ -76,9 +76,9 @@ Return a dictionary which assciociates to each observable of the set `observable
 an estimation of its expected value.
 """
 function estimate_set_obs(
-    shadows::Array{Shadow}, observables::Set{String}, Ne::Int, Ke::Int
+    shadows::Array{Shadow}, observables::Set{PauliObs}, Ne::Int, Ke::Int
 )::Dict{String,Float64}
-    Dict(o => estimate_obs(shadows, o, Ne, Ke) for o in observables)
+    Dict(PauliObs_to_string(o) => estimate_obs(shadows, o, Ne, Ke) for o in observables)
 end
 
 """
@@ -90,9 +90,9 @@ an estimation of its expected value.
 NOTE(rmartin): Calculating each observable's weight (instead of taking the max) and compute only the usefull bit strings may be faster in some cases.
 """
 function robust_estimate_set_obs(
-    shadows::Array{Shadow}, calibration::Calibration, observables::Set{String}, Ne::Int, Ke::Int
+    shadows::Array{Shadow}, calibration::Calibration, observables::Set{PauliObs}, Ne::Int, Ke::Int
 )::Dict{String,Float64}
-    Dict(o => robust_estimate_obs(calibration, shadows, o, Ne, Ke) for o in observables)
+    Dict(PauliObs_to_string(o) => robust_estimate_obs(calibration, shadows, o, Ne, Ke) for o in observables)
 end
 
 """
@@ -102,7 +102,7 @@ Return an estimation of the `observable`'s expected value taking noise into acco
 It is computed as the median of `K` means, each one being a mean of `N` single estimation.
 """
 function robust_estimate_obs(
-    calibration::Calibration, shadows::Array{Shadow}, observable::String, Ne::Int, Ke::Int
+    calibration::Calibration, shadows::Array{Shadow}, observable::PauliObs, Ne::Int, Ke::Int
 )::Float64
     R = Ne * Ke
     median_of_mean(
@@ -115,42 +115,24 @@ end
 
 Based on a single `shadow` and noise correlation factors `f`, compute a single estimation
 which equals the `observable`'s expected value in average.
-
-The inverse of the quantum channel representing the measurement protocol is a sum of terms.
-Each term is associated to a bitstring of {0, 1}ⁿ and must be multiplied by the corresponding correlation factor of `f`.
-The only terms that contribute to the estimation are those for which the corresponding bitstring
-does not have a "1" in the same position as an "I" in the `observable`.
-This justifies why only bitstring with less than k "1" have been considered (where k is the maximum weight of observables)
-
-NOTE(rmartin): Computing robust_single_estimate for each shadow is the bottleneck of the
-implementation because its complexity is proportional to nᵏ.
 """
 function robust_single_estimate(
-    calibration::Calibration, shadow::Shadow, observable::String
+    calibration::Calibration, shadow::Shadow, observable::PauliObs
 )::Float64
     n = shadow.n_qubits
-    estimations = 0
     f = calibration.coefficients
-    for bit_string in keys(f)
-        product = 1 / f[bit_string]
-        for j = 1:n
-            if observable[j] == 'I'
-                if bit_string[j] == 1 #The corresponding term does not contribute.
-                    product = 0
-                    break
-                end #No "else" beacause in the other case we only have to multiply product by one (I expected value)
-            else
-                if bit_string[j] == 0
-                    product *= tr(Matrix(mat(pauli_from_char[observable[j]])))
-                else
-                    U = Matrix(mat(shadow.gate[j]))
-                    b = shadow.measure[j] == 0 ? [1, 0] : [0, 1]
-                    Oj = Matrix(mat(pauli_from_char[observable[j]]))
-                    product *= b' * U * Oj * U' * b - tr(Oj) / 2
-                end
-            end
-        end
-        estimations += product
+    bit_string = zeros(Int, n)
+    terms = observable.pauli
+    product = 1
+    for j in keys(terms)
+        bit_string[j] = 1
+        U = Matrix(mat(shadow.gate[j]))
+        b = shadow.measure[j] == 0 ? [1, 0] : [0, 1]
+        Oj = Matrix(mat(pauli_from_char[terms[j]]))
+        product *= b' * U * Oj * U' * b 
     end
-    return real(estimations)
+
+    product /= f[bit_string]
+
+    return real(product)
 end

src/observables.jl

@@ -18,6 +18,46 @@ Enforces that it consists of Pauli operators acting on different qubits.
 pauli_string(n::Int, vs::Vararg{Pair{Int,<:PauliGate}}) = Yao.kron(n, vs...)
 pauli_string(vs::Vararg{<:PauliGate}) = Yao.kron(vs...)
 
+"""
+    PauliObs
+
+Represent a pauli observable storing only non identity terms.
+`n_qubits` is the number of qubits on which the observable acts.
+`pauli` is a dictionary which keys are the position of the non identity terms
+and values are the character refering to the pauli matrix (namely 'X', 'Y' or 'Z').
+"""
+struct PauliObs
+    n_qubits::Int
+    pauli::Dict{Int64, Char}
+
+    function PauliObs(obs::String)::PauliObs
+        n_qubits = length(obs)
+        pauli = Dict{Int64, Char}()
+        for i in 1:n_qubits
+            if obs[i] != 'I'
+                pauli[i] = obs[i]
+            end
+        end
+        new(n_qubits, pauli)
+    end
+
+    function Pauli(n_qubits::Int, pauli::Dict{Int64, Char})
+        new(n_qubits, pauli)
+    end
+end
+
+function PauliObs_to_string(obs::PauliObs)::String
+    string = collect("I"^obs.n_qubits)
+    for p in obs.pauli
+        string[p[1]] = p[2]
+    end
+    join(string)
+end
+
+function weight(pauli_observable::PauliObs)::Int
+    length(pauli_observable.pauli)
+end
+
 struct PauliObservable
     n_qubits::Int
     terms::Vector{AbstractBlock}