# Sampling estimation

In order to estimate expectation value of operators on a quantum computer, sampling measurements are necessary. In sampling measurements, execution of a quantum circuit and a subsequent measurement of qubits are repeated multiple times. Estimation of expectation value of operators is then performed using statistics of the repeated measurements.

Before performing sampling estimation on real quantum computers, let's see how such procedures can be simulated with a quantum circuit simulator. In this tutorial, we will introduce how to

• Perform grouping for Operators.
• Construct measurement circuit from grouping result.
• Allocate shots for measurement circuits and perform sampling.
• Reconstruct expectation value from sampling result.

## Prerequisite​

QURI Parts modules used in this tutorial: quri-parts-circuit, quri-parts-core and quri-parts-qulacs. You can install them as follows:

!pip install "quri-parts[qulacs]"

## Overview​

The major purpose of this tutorial is to introduce how to perform sampling estimation with a sampling estimator. A sampling estimator is a QuantumEstimator where the expectation value is computed via a sampling experiment. As a preparation, we create an operator and a state we want to estimate the expectation value with.

from quri_parts.core.operator import Operator, pauli_label, PAULI_IDENTITYfrom quri_parts.core.state import quantum_stateop = Operator({    pauli_label("Z0"): 0.25,    pauli_label("Z1 Z2"): 2.0,    pauli_label("X1 X2"): 0.5 + 0.25j,    pauli_label("Z1 Y3"): 1.0j,    pauli_label("Z2 Y3"): 1.5 + 0.5j,    pauli_label("X1 Y3"): 2.0j,    PAULI_IDENTITY: 3.0,})print("Operator:")print(op)state = quantum_state(4, bits=0b0101)print("")print("State:")print(state)
#output    Operator:    0.25*Z0 + 2.0*Z1 Z2 + (0.5+0.25j)*X1 X2 + 1j*Z1 Y3 + (1.5+0.5j)*Z2 Y3 + 2j*X1 Y3 + 3.0*I        State:    ComputationalBasisState(qubit_count=4, bits=0b101, phase=0π/2)

Here, we provide a sample code for creating and using a sampling estimator. We create a sampling estimator with the following specification:

• Total shots: 5000
• measurement factory: Bitwise commuting Pauli grouping
• Shot allocator: proportional shots allocator

We will explain what "measurement factory" and "Shot allocator" mean in later sections. For now, let's temporarily take them for granted and create a sampling estimator with them:

from quri_parts.core.estimator.sampling import create_sampling_estimatorfrom quri_parts.core.measurement import bitwise_commuting_pauli_measurementfrom quri_parts.core.sampling.shots_allocator import create_proportional_shots_allocatorfrom quri_parts.qulacs.sampler import create_qulacs_vector_concurrent_samplertotal_shots = 5000concurrent_sampler = create_qulacs_vector_concurrent_sampler()measurement_factory = bitwise_commuting_pauli_measurementshots_allocator = create_proportional_shots_allocator()sampling_estimator = create_sampling_estimator(    total_shots,    concurrent_sampler,    measurement_factory,    shots_allocator)

Let's look at the result of the sampling estimate:

estimate = sampling_estimator(op, state)print("Estimated expectation value:", estimate.value)print("Standard error of this sampling estimation:", estimate.error)
#output    Estimated expectation value: (0.6527269558463302+0.02582389806961259j)    Standard error of this sampling estimation: 0.070823557319059

We may compare the sampling estimate's result with the exact result by a vector estimator.

from quri_parts.qulacs.estimator import create_qulacs_vector_estimatorvector_estimator = create_qulacs_vector_estimator()vector_estimator(op, state)
#output    _Estimate(value=(0.75+0j), error=0.0)

## Grouping operators into commutable Pauli strings​

In sampling experiments, one of the ways to estimate expectation value of such an operator is to estimate the expectation value of each Pauli string and then sum them up. However, it is more efficient to measure multiple Pauli strings at once if they are commutable. The first step is thus to group the Pauli strings into several sets of commutable Pauli strings. This Pauli grouping is an important research subject in context of operator estimation.

Then, from the grouping results, we need to construct a circuit for each commutable group for measurement and sampling, which we call a measurement circuit. The explicit gates involved in a measurement circuit will depend on the content of the commutable group of Pauli strings.

Follow this structure, we use the operator and state constructed above:

print("Operator:")for pl, coeff in op.items():    print(f"   + {coeff} * {pl}")print("")print("State:", state)
#output    Operator:    + 0.25 * Z0    + 2.0 * Z1 Z2    + (0.5+0.25j) * X1 X2    + 1j * Z1 Y3    + (1.5+0.5j) * Z2 Y3    + 2j * X1 Y3    + 3.0 * I        State: ComputationalBasisState(qubit_count=4, bits=0b101, phase=0π/2)

to illustrate how to

• perform operator grouping
• construct measurement circuit
• reconstruct expectation of Pauli strings from measurement results

in QURI Parts. Finally, we will wrap this section up by introducing a convenient object: CommutablePauliSetMeasurement, that encapsulates all the steps above.

### Pauli grouping​

In QURI Parts, a grouping function is represented by PauliGrouping. They are functions that take an Operator or a collection of PauliLabel into sets of CommutablePauliSet, where CommutablePauliSet is a set of Pauli strings that commute with each other. The function signature is given by:

from typing import Callable, Set, Union, Iterablefrom typing_extensions import TypeAliasfrom quri_parts.core.operator import PauliLabelCommutablePauliSet: TypeAlias = Set[PauliLabel]#: Represents a grouping functionPauliGrouping = Callable[    [Union[Operator, Iterable[PauliLabel]]], Iterable[CommutablePauliSet]]

As an example, we use one of the simplest Pauli grouping: bitwise grouping, where the groups are determined based on bitwise commutability of the Pauli strings. We can test the grouping as follows:

from quri_parts.core.operator.grouping import bitwise_pauli_groupingpauli_sets = bitwise_pauli_grouping(op)print(pauli_sets)
#output    frozenset({frozenset({PauliLabel({(1, <SinglePauli.X: 1>), (3, <SinglePauli.Y: 2>)})}), frozenset({PauliLabel()}), frozenset({PauliLabel({(1, <SinglePauli.X: 1>), (2, <SinglePauli.X: 1>)})}), frozenset({PauliLabel({(0, <SinglePauli.Z: 3>)}), PauliLabel({(2, <SinglePauli.Z: 3>), (1, <SinglePauli.Z: 3>)})}), frozenset({PauliLabel({(2, <SinglePauli.Z: 3>), (3, <SinglePauli.Y: 2>)}), PauliLabel({(3, <SinglePauli.Y: 2>), (1, <SinglePauli.Z: 3>)})})})

The above looks a bit complicated since the grouping method returns a frozenset of frozensets of Pauli labels.

print(f"Number of groups: {len(pauli_sets)}")for i, pauli_set in enumerate(pauli_sets):    labels = ", ".join([str(pauli) for pauli in pauli_set])    print(f"Group {i} contains: {labels}")
#output    Number of groups: 5    Group 0 contains: X1 Y3    Group 1 contains: I    Group 2 contains: X1 X2    Group 3 contains: Z0, Z1 Z2    Group 4 contains: Z2 Y3, Z1 Y3

Here, we provide a list of available grouping functions in QURI Parts:

Grouping functionExplanation
individual_pauli_groupingNo grouping
bitwise_pauli_groupingCheck Pauli string commutability based on bitwise commutability
sorted_injection_groupingExplained in arXiv:1908.06942

### Measurement circuit​

To perform a measurement for a commutable Pauli set, it is necessary to construct a circuit to be applied before measurement. In the case of bitwise Pauli grouping, the corresponding circuit can be constructed with the bitwise_commuting_pauli_measurement_circuit function, which we show as follows:

from quri_parts.core.measurement import bitwise_commuting_pauli_measurement_circuitfrom quri_parts.circuit import QuantumCircuitfrom quri_parts.circuit.utils.circuit_drawer import draw_circuitpauli_set = {pauli_label("X0 Z2 Y3"), pauli_label("X0 Z1 Y3")}measurement_circuit = bitwise_commuting_pauli_measurement_circuit(pauli_set)draw_circuit(QuantumCircuit(4, gates=measurement_circuit))
#output         ___                  | H |             0 --|0  |---------        |___|                                                         1 ----------------                                                                      2 ----------------                               ___     ___          |Sdg|   | H |     3 --|1  |---|2  |-        |___|   |___| 

### Sampling on a state​

In order to measure a commutable Pauli set for a quantum state, the following steps are necessary:

• Construct a circuit that prepares the state
• Concatenate the measurement circuit for the Pauli set after the state preparation circuit
• Perform sampling on the concatenated circuit

Here we use a ComputationalBasisState as the initial state for simplicity, though you can use any CircuitQuantumState.

from quri_parts.core.state import quantum_statefrom quri_parts.qulacs.sampler import create_qulacs_vector_samplersampler = create_qulacs_vector_sampler()initial_state = quantum_state(4, bits=0b0101)# Circuit for state preparationstate_prep_circuit = initial_state.circuit# Measurement circuitpauli_set = {pauli_label("Z2 Y3"), pauli_label("Z1 Y3")}measurement_circuit = bitwise_commuting_pauli_measurement_circuit(pauli_set)# Concatenate measurement circuitsampled_circuit = state_prep_circuit + measurement_circuitprint("State preparation circuit concatenated with measurement circuit:")draw_circuit(sampled_circuit)# Samplingsampling_result = sampler(sampled_circuit, shots=1000)print("\n")print("Sampling result:")print({bin(bits): count for bits, count in sampling_result.items()})
#output    State preparation circuit concatenated with measurement circuit:         ___                  | X |             0 --|0  |---------        |___|                                                         1 ----------------                               ___                  | X |             2 --|1  |---------        |___|                  ___     ___          |Sdg|   | H |     3 --|2  |---|3  |-        |___|   |___|     Sampling result:    {'0b101': 501, '0b1101': 499}

### Reconstructing Pauli string expectation values from sampled values​

It is then necessary to read off the measured eigenvalue of a Pauli string from the measurement result. In QURI Parts, this is done with a PauliReconstructorFactory and a PauliReconstructor. A PauliReconstructor is a function that returns the eigenvalue of a Pauli string that corresponds to the measurement result, and a PauliReconstructorFactory is a function that creates a PauliReconstructor from a PauliLabel. Here is their function signatures defined in QURI Parts.

#: PauliReconstructor represents a function that reconstructs a value of a Pauli#: operator from a measurement result of its measurement circuit.PauliReconstructor: TypeAlias = Callable[[int], int]#: PauliReconstructorFactory represents a factory function that returns#: a :class:~PauliReconstructor for a given Pauli operator.PauliReconstructorFactory: TypeAlias = Callable[[PauliLabel], PauliReconstructor]

In the example of the above section, we performed a sampling measurement for $Z_2 Y_3$ and $Z_1 Y_3$, and obtained two bit patterns 0b1101 and 0b0101. In the case of bitwise Pauli grouping, a value of a Pauli operator can be reconstructed from a sampled bit pattern as follows:

from quri_parts.core.measurement import bitwise_pauli_reconstructor_factory# Obtain a reconstructor for Z2 Y3reconstructor = bitwise_pauli_reconstructor_factory(pauli_label("Z2 Y3"))# Reconstruct a value of Z2 Y3 from a sampled bit pattern 0b1101pauli_value = reconstructor(0b1101)print("Measurement result: 0b1101.", f"Corresponding eigenvalue: {pauli_value}")# Reconstruct from 0b0101pauli_value = reconstructor(0b0101)print("Measurement result: 0b0101.", f"Corresponding eigenvalue: {pauli_value}")
#output    Measurement result: 0b1101. Corresponding eigenvalue: 1    Measurement result: 0b0101. Corresponding eigenvalue: -1

The expectation value of $Z_2 Y_3$ can then be calculated as:

pauli_exp = (    reconstructor(0b1101) * sampling_result[0b1101] +    reconstructor(0b0101) * sampling_result[0b0101]) / 1000print(pauli_exp)# Equivalent to abovepauli_exp = sum(    reconstructor(bits) * count for bits, count in sampling_result.items()) / sum(sampling_result.values())print(pauli_exp)
#output    -0.002    -0.002

QURI Parts also provide the trivial_pauli_expectation_estimator function as a convenient way of obtaining expection value of PauliLabel from sampling results.

from quri_parts.core.estimator.sampling import trivial_pauli_expectation_estimatorpauli_exp = trivial_pauli_expectation_estimator(sampling_result, pauli_label("Z2 Y3"))print(pauli_exp)
#output    -0.002

Here trivial_pauli_expectation_estimator can be used because we used bitwise Pauli grouping. In a more general case, general_pauli_expectation_estimator should be used with a PauliReconstructorFactory:

from quri_parts.core.estimator.sampling import general_pauli_expectation_estimatorpauli_exp = general_pauli_expectation_estimator(    sampling_result, pauli_label("Z2 Y3"), bitwise_pauli_reconstructor_factory)print(pauli_exp)
#output    -0.002

Finally, expectation value of the original operator can be estimated by summing up contributions of each Pauli term. To calculate contribution of $Z_2 Y_3$:

# Coefficient of Z2 Y3 in opcoef = op[pauli_label("Z2 Y3")]pauli_contrib = coef * pauli_expprint(pauli_contrib)
#output    (-0.003-0.001j)

Repeating this procedure for all Pauli strings, expectation value of the original operator can be estimated.

### The CommutablePauliSetMeasurement object​

The above procedure is a bit complicated, so a shortcut method is provided. To use them we first introduce CommutablePauliSetMeasurement: it is a data structure containing the following elements:

• pauli_set: A set of commutable Pauli operators that are measured together.
• measurement_circuit: A circuit required to measure the given pauli_set.
• pauli_reconstructor_factory: A factory of functions that reconstruct Pauli operator values from sampling result

A set of CommutablePauliSetMeasurement needs to be constructed depending on a specific measurement scheme. They are usually obtained by a CommutablePauliSetMeasurementFactory, which returns a sequence of CommutablePauliSetMeasurement by taking in an Operator or a sequence of PauliLabels.

For example, for bitwise Pauli grouping, you can do as follows:

from quri_parts.core.measurement import bitwise_commuting_pauli_measurementmeasurements = bitwise_commuting_pauli_measurement(op)print(f"Number of CommutablePauliSetMeasurement: {len(measurements)}")measurement = measurements[-1]
#output    Number of CommutablePauliSetMeasurement: 5
• Commutable Pauli Set
print(", ".join(map(lambda s: str(s), measurement.pauli_set)))
#output    Z2 Y3, Z1 Y3
• State preparation circuit + Measurement circuit

We show the circuit for estimating $\langle Z_2 Y_3 \rangle$ and $\langle Z_1 Y_3 \rangle$, and then sample with 1000 shots from it.

full_circuit = state.circuit + measurement.measurement_circuitdraw_circuit(full_circuit)samling_result = sampler(full_circuit, 1000)print("\n")print(f"Measurement Result: {samling_result}")
#output         ___                  | X |             0 --|0  |---------        |___|                                                         1 ----------------                               ___                  | X |             2 --|1  |---------        |___|                  ___     ___          |Sdg|   | H |     3 --|2  |---|3  |-        |___|   |___|             Measurement Result: Counter({13: 501, 5: 499})
• Reconstructor factory

We retrieve the PauliReconstructorFactory from measurement and then compute $\langle Z_2 Y_3 \rangle$ and $\langle Z_1 Y_3 \rangle$

print(f"Reconstructor factory: {measurement.pauli_reconstructor_factory}")for pl in measurement.pauli_set:    estimate = general_pauli_expectation_estimator(        sampling_result, pl, measurement.pauli_reconstructor_factory    )    print(f" <{pl}> = {estimate}")    
#output    Reconstructor factory: <function bitwise_pauli_reconstructor_factory at 0x127af98b0>     <Z2 Y3> = -0.002     <Z1 Y3> = 0.002

We now list out all the PauliReconstructorFactorys provided in QURI Parts:

• bitwise_commuting_pauli_measurement: Grouping based on bitwise commutation.
• individual_pauli_measurement: No grouping.

## Shot Allocators​

Up until this point, we have learned how to group an Operator into several CommutablePauliSets. We also learned how to evaluate expectation values of individual PauliLabel by reading off sampling result from a measurement circuit. The final step necessary for estimation is PauliSamplingShotsAllocator: it specifies how total sampling shots should be allocated to measurement of each Pauli sets.

### Interface​

In QURI Parts, shot allocators used in sampling estimation are represented by PauliSamplingShotsAllocator [1]. They distribute speficied total shot number to CommutablePauliSet according to the Operator under consideration or the set of CommutablePauliSet. The function signature is given by

from typing import Collectionfrom quri_parts.core.sampling import PauliSamplingSettingCommutablePauliSetMeasurement: TypeAlias = Callable[    [Operator, Collection[CommutablePauliSet], int], Collection[PauliSamplingSetting]]

where the returned PauliSamplingSetting is a named tuple holding a CommutablePauliSet and the number of shots assigned to it.

[1] ^ There is another type of shot allocator WeightedSamplingShotsAllocator provided by QURI Parts. However, they are not used for sampling estimation. So, we will not introduce them in this tutorial.

### Available shot allocators​

Finally, we introduce the shots allocator available in QURI Parts:

from quri_parts.core.sampling.shots_allocator import (    create_equipartition_shots_allocator,    create_proportional_shots_allocator,    create_weighted_random_shots_allocator,)# Allocates shots equally among the Pauli setsallocator = create_equipartition_shots_allocator()# Allocates shots proportional to Pauli coefficients in the operatorallocator = create_proportional_shots_allocator()# Allocates shots using random weightsallocator = create_weighted_random_shots_allocator(seed=777)

We also summarize the allocators below:

AllocatorGenerating functionExplanationReference
equipartition shots allocatorcreate_equipartition_shots_allocatorAllocates shots equally among the Pauli sets-
proportional shots allocatorcreate_proportional_shots_allocatorAllocates shots proportional to Pauli coefficients in the operator- arXiv:2112.07416 (2021). See Appendix D.
- Phys. Rev. A 92, 042303
weighted random shots allocatorcreate_weighted_random_shots_allocatorAllocates shots using random weightsarXiv:2004.06252 (2020)

## Sampling estimation​

It has been a long way introducing all the necessary components required to perform sampling estimation. In QURI Parts, we provide a sampling_estimate function for you to customize the sampling estimation computation. The required elements for estimating the expectation value of an operator for a CircuitQuantumState are:

• Total shots: The total number of shots you use to estimate the expectation value
• Concurrent sampler: A ConcurrentSampler that performs the sampling for all (state preparation circuit + measurement circuit).
• Measurement factory: A CommutablePauliSetMeasurementFactory that performs grouping for the input operator and evaluates all the expectation values of Pauli strings involved.
• Shots allocator: A PauliSamplingShotsAllocator that distributes the total shots among CommutablePauliSet generated by the measurement factory.

We now show how to use it.

from quri_parts.qulacs.sampler import create_qulacs_vector_concurrent_samplerfrom quri_parts.core.estimator.sampling import sampling_estimateconcurrent_sampler = create_qulacs_vector_concurrent_sampler()estimate = sampling_estimate(    op,            # Operator to estimate    state,         # A circuit state    5000,          # Total sampling shots    concurrent_sampler, # ConcurrentSampler    bitwise_commuting_pauli_measurement, # Factory function for CommutablePauliSetMeasurement    allocator,     # PauliSamplingShotsAllocator)print(f"Estimated expectation value: {estimate.value}")print(f"Standard error of estimation: {estimate.error}")
#output    Estimated expectation value: (0.7812513815547097+0.002275914633594121j)    Standard error of estimation: 0.0707060360741753

Most of the time, you might want to embed an estimator or a concurrent estimator into an algorithm. QURI Parts provide create_sampling_estimator and create_sampling_concurrent_estimator to create (Concurrent)QuantumEstimator that can be used like a vector estimator introduced in the estimator tutorial.

from quri_parts.core.estimator.sampling import create_sampling_estimatorestimator = create_sampling_estimator(    5000,          # Total sampling shots    concurrent_sampler, # ConcurrentSampler    bitwise_commuting_pauli_measurement, # Factory function for CommutablePauliSetMeasurement    allocator,     # PauliSamplingShotsAllocator)estimate = estimator(op, state)print(f"Estimated expectation value: {estimate.value}")print(f"Standard error of estimation: {estimate.error}")
#output    Estimated expectation value: (0.7533110100454796-0.0441731501615594j)    Standard error of estimation: 0.07035284315335183

For concurrent sampling estimators, we create another set of operator and state:

from quri_parts.core.state import GeneralCircuitQuantumStatefrom quri_parts.circuit import H, QuantumCircuitop_1 = opop_2 = Operator({    PAULI_IDENTITY: 8,    pauli_label("X0 Y1 Z2 X3"): 9 + 8j})state_1 = statestate_2 = GeneralCircuitQuantumState(4, QuantumCircuit(4, gates=[H(i) for i in range(4)]))
from quri_parts.core.estimator.sampling import create_sampling_concurrent_estimatorconcurrent_estimator = create_sampling_concurrent_estimator(    5000,          # Total sampling shots    concurrent_sampler, # ConcurrentSampler    bitwise_commuting_pauli_measurement, # Factory function for CommutablePauliSetMeasurement    allocator,     # PauliSamplingShotsAllocator)estimates = concurrent_estimator([op_1, op_2], [state_1, state_2])print(f"Estimated expectation value 1: {estimates[0].value}")print(f"Standard error of estimation 1: {estimates[0].error}")print("")print(f"Estimated expectation value 2: {estimates[1].value}")print(f"Standard error of estimation 2: {estimates[1].error}")
#output    Estimated expectation value 1: (0.7328563810878067+0.02033531023801298j)    Standard error of estimation 1: 0.0711939414351656        Estimated expectation value 2: (7.8092-0.1696j)    Standard error of estimation 2: 0.1702555909214144