diff --git a/src/classical_shadows.jl b/src/classical_shadows.jl index be0b6e2..208e01e 100644 --- a/src/classical_shadows.jl +++ b/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 """ diff --git a/src/estimate.jl b/src/estimate.jl index 82b8256..5916c42 100644 --- a/src/estimate.jl +++ b/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 diff --git a/src/estimate_from_shadows.jl b/src/estimate_from_shadows.jl index c264bb3..1b87005 100644 --- a/src/estimate_from_shadows.jl +++ b/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 diff --git a/src/observable_estimation.jl b/src/observable_estimation.jl index fb2d7f8..4831866 100644 --- a/src/observable_estimation.jl +++ b/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 diff --git a/src/observables.jl b/src/observables.jl index b31365e..342b921 100644 --- a/src/observables.jl +++ b/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}