Sampling simulation#

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.

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]"

Prepare a circuit#

As a preparation, we create a circuit to be sampled:

[1]:

from math import pi
from quri_parts.circuit import QuantumCircuit
# A circuit with 4 qubits
circuit = QuantumCircuit(4)
circuit.add_X_gate(0)
circuit.add_H_gate(1)
circuit.add_Y_gate(2)
circuit.add_CNOT_gate(1, 2)
circuit.add_RX_gate(3, pi/4)

Sampler#

When performing a sampling measurement for a circuit, you can use a Sampler. The Sampler itself (defined in quri_parts.core.sampling) is an abstract interface and you need a concrete instance to actually perform sampling. There are several implementations of Sampler interface, some of which use a circuit simulator while others use a real quantum computer.

Here we start with a sampler using state vector simulation with Qulacs:

[2]:

from quri_parts.qulacs.sampler import create_qulacs_vector_sampler
# Create the sampler
sampler = create_qulacs_vector_sampler()
sampling_result = sampler(circuit, shots=1000)
print(sampling_result)

Counter({5: 434, 3: 411, 13: 90, 11: 65})

A sampler receives two arguments: a circuit to be sampled and a number of repeated samplings (shots). The returned value is a mapping (dict) with the following keys and values:

  • Keys Bitstrings of measurement outcomes, encoded as int. Each measured bit of each qubit are ordered from the least significant bit to the most significant bit. For example, when qubit 0 and qubit 2 are measured to be in \(|1\rangle\) while the others are in \(|0\rangle\), the bitstring is 0b0101.

  • Values Counts of times for which each bitstring is measured. A sum of all counts is equal to the specified shots.

[3]:

for bits, count in sampling_result.items():
    print(f"A bitstring '{bin(bits)}' is measured {count} times")
print(f"Total count is {sum(sampling_result.values())}")

A bitstring '0b11' is measured 411 times
A bitstring '0b101' is measured 434 times
A bitstring '0b1101' is measured 90 times
A bitstring '0b1011' is measured 65 times
Total count is 1000

In the above example, we can see that the qubit 0 is always measured as in \(|1\rangle\) since the only gate applied to the qubit is \(X\). On the other hand qubit 1, 2 and 3 can take either 0 or 1, but since qubit 1 and 2 are correlated by a CNOT gate, only 4 patterns are in the measurement result.

Pauli grouping#

Estimation of expectation value of operators involves operators and states, which cannot be handled directly by a sampler. It is thus necessary to translate the problem to circuits, which can be sampled by a sampler.

First let’s define an operator to be estimated:

[4]:

from quri_parts.core.operator import Operator, pauli_label, PAULI_IDENTITY
op = 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(op)

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

The operator is represented as a sum of Pauli operators. One of the ways to estimate expectation value of such an operator is to estimate expectation value of each Pauli term and then sum up them.

When estimating the Pauli terms, it is possible to measure multiple Pauli terms at once if they are commutable. The first step is thus to group the Pauli terms into several sets of commutable Pauli terms. This Pauli grouping is an important research subject in context of operator estimation.

One of the simplest Pauli grouping is bitwise grouping, where the groups are determined based on bitwise commutability of the Pauli terms. We can test the grouping as follows:

[5]:

from quri_parts.core.operator.grouping import bitwise_pauli_grouping
pauli_sets = bitwise_pauli_grouping(op)
print(pauli_sets)

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.

[6]:

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}")

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

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 as follows:

[7]:

from quri_parts.core.measurement import bitwise_commuting_pauli_measurement_circuit
pauli_set = {pauli_label("Z2 Y3"), pauli_label("Z1 Y3")}
measurement_circuit = bitwise_commuting_pauli_measurement_circuit(pauli_set)
print(measurement_circuit)

(QuantumGate(name='Sdag', target_indices=(3,), control_indices=(), params=(), pauli_ids=()), QuantumGate(name='H', target_indices=(3,), control_indices=(), params=(), pauli_ids=()))

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.

[8]:

from quri_parts.core.state import ComputationalBasisState
initial_state = ComputationalBasisState(4, bits=0b0101)
# Circuit for state preparation
state_prep_circuit = initial_state.circuit
# Concatenate measurement circuit
sampled_circuit = state_prep_circuit + measurement_circuit
# Sampling
sampling_result = sampler(sampled_circuit, shots=1000)
print({bin(bits): count for bits, count in sampling_result.items()})

{'0b101': 500, '0b1101': 500}

Reconstructing Pauli values from sampled values#

It is then necessary to reconstruct values of Pauli terms from the sampling result. In the above example, 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:

[9]:

from quri_parts.core.measurement import bitwise_pauli_reconstructor_factory
# Obtain a reconstructor for Z2 Y3
reconstructor = bitwise_pauli_reconstructor_factory(pauli_label("Z2 Y3"))
# Reconstruct a value of Z2 Y3 from a sampled bit pattern 0b1101
pauli_value = reconstructor(0b1101)
print(pauli_value)
# Reconstruct from 0b0101
pauli_value = reconstructor(0b0101)
print(pauli_value)

1
-1

The expectation value of \(Z_2 Y_3\) can then be calculated as:

[10]:

pauli_exp = (
    reconstructor(0b1101) * sampling_result[0b1101] +
    reconstructor(0b0101) * sampling_result[0b0101]
) / 1000
print(pauli_exp)

# Equivalent to above
pauli_exp = sum(
    reconstructor(bits) * count for bits, count in sampling_result.items()
) / sum(sampling_result.values())
print(pauli_exp)

# More convenient way
from quri_parts.core.estimator.sampling import trivial_pauli_expectation_estimator
pauli_exp = trivial_pauli_expectation_estimator(sampling_result, pauli_label("Z2 Y3"))
print(pauli_exp)

0.0
0.0
0.0

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:

[11]:

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

0.0

Estimate original operator expectation value from estimated Pauli values#

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\):

[12]:

# Coefficient of Z2 Y3 in op
coef = op[pauli_label("Z2 Y3")]
pauli_contrib = coef * pauli_exp
print(pauli_contrib)

0j

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

Sampling estimation#

The above procedure is a bit complicated, so a shortcut methods are 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: for example, for bitwise Pauli grouping, you can do as follows:

[13]:

from quri_parts.core.measurement import bitwise_commuting_pauli_measurement
measurements = bitwise_commuting_pauli_measurement(op)
print(f"Number of CommutablePauliSetMeasurement: {len(measurements)}")
measurement = measurements[0]
print(measurement.pauli_set)
print(measurement.measurement_circuit)
print(measurement.pauli_reconstructor_factory)

Number of CommutablePauliSetMeasurement: 5
frozenset({PauliLabel({(1, <SinglePauli.X: 1>), (3, <SinglePauli.Y: 2>)})})
(QuantumGate(name='H', target_indices=(1,), control_indices=(), params=(), pauli_ids=()), QuantumGate(name='Sdag', target_indices=(3,), control_indices=(), params=(), pauli_ids=()), QuantumGate(name='H', target_indices=(3,), control_indices=(), params=(), pauli_ids=()))
<function bitwise_pauli_reconstructor_factory at 0x7f567494e9d0>

Another input necessary for estimation is PauliSamplingShotsAllocator: it specifies how total sampling shots should be allocated to measurement of each Pauli sets. There are several allocators available:

[14]:

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 sets
allocator = create_equipartition_shots_allocator()
# Allocates shots proportional to Pauli coefficients in the operator
allocator = create_proportional_shots_allocator()
# Allocates shots using random weights
allocator = create_weighted_random_shots_allocator(seed=777)

With these inputs, sampling estimation can be performed as follows:

[15]:

from quri_parts.qulacs.sampler import create_qulacs_vector_concurrent_sampler
from quri_parts.core.estimator.sampling import sampling_estimate
concurrent_sampler = create_qulacs_vector_concurrent_sampler()
estimate = sampling_estimate(
    op,            # Operator to estimate
    initial_state, # Initial (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}")

Estimated expectation value: (0.6911126120594374+0.021229238383697557j)
Standard error of estimation: 0.07071572965514332

You can also create a QuantumEstimator that performs sampling estimation:

[16]:

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

Estimated expectation value: (0.7443525413563896-0.06931160253571916j)
Standard error of estimation: 0.07035696518563977