Due to the imperfection of real devices, errors occur in state preparation and measurement. Readout error mitigation reduces the effect of those errors by applying inverse of such noisy operations. The inverse of the noisy operation, here we call filter matrix, plays an important role in readout error mitigation. Noisy counts $C_{\text{noisy}}$ can be considered as a product of ideal counts $C_{\text{ideal}}$ which we could get in noiseless world and the error matrix $E$.

$C_{\text{noisy}} = EC_{\text{ideal}}$

The $(i, j)$ components of $E$ represent the probability that the true state is $|j\rangle$ when the measurement reports $|i\rangle$ as the outcome, which is represented as $P(i|j)$. That is, summing vertically over the matrix $E$ yields 1. For exsample, $E$ of 2 qubits device is

$E = \left( \begin{matrix} P(00|00) & P(00|01) & P(00|10) & P(00|11) \\ P(01|00) & P(01|01) & P(01|10) & P(01|11) \\ P(10|00) & P(10|01) & P(10|10) & P(10|11) \\ P(11|00) & P(11|01) & P(11|10) & P(11|11) \end{matrix} \right)$

$E$ is estimated by repeating state preparation and measurement in each basis and obtaining probability distributions from histograms of actual measurement results.

Filter matrix is defined as the inverse of error matrix. With filter matrix we can estimate the error-free counts

$C_{\text{ideal}} = E^{-1}C_{\text{noisy}}.$

In this tutorial, we demonstrate how to build the filter matrix and predict the noise-free sampling count.

## Prerequisite​

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

# !pip install "quri-parts[qulacs]"

## Preparation and overview​

Here, we prepare the circuit and the noise model we use throughout this tutorial. The circuit we use in this tutorial consists of an identity part and a non-trivial part. The non-trivial part is responsible for converting the state $|000\rangle$ into $\frac{1}{\sqrt{2}}\left(|000\rangle + |111\rangle\right)$, while we decompose the identity circuit into multiple gates to amplify the effect of the noises.

from quri_parts.circuit import QuantumCircuitfrom quri_parts.circuit.utils.circuit_drawer import draw_circuitqubit_count = 3identity_circuit = QuantumCircuit(3)identity_circuit.add_RX_gate(0, 1.3)identity_circuit.add_RY_gate(1, 0.2)identity_circuit.add_RZ_gate(0, -2.3)identity_circuit.add_SqrtXdag_gate(1)identity_circuit.add_T_gate(0)identity_circuit.add_RX_gate(1, 0.4)identity_circuit.add_RY_gate(0, 2.7)identity_circuit.add_Tdag_gate(1)identity_circuit.add_RY_gate(0, -2.7)identity_circuit.add_T_gate(1)identity_circuit.add_Tdag_gate(0)identity_circuit.add_RX_gate(1, -0.4)identity_circuit.add_RZ_gate(0, 2.3)identity_circuit.add_SqrtX_gate(1)identity_circuit.add_RX_gate(0, -1.3)identity_circuit.add_RY_gate(1, -0.2)circuit = QuantumCircuit(3)circuit += identity_circuitcircuit.add_H_gate(0)circuit.add_CNOT_gate(0, 1)circuit.add_CNOT_gate(0, 2)print("The circuit:")draw_circuit(circuit, line_length=200)
#output    The circuit:       ___     ___     ___     ___     ___     ___     ___     ___     ___                        |RX |   |RZ |   | T |   |RY |   |RY |   |Tdg|   |RZ |   |RX |   | H |                     --|0  |---|2  |---|4  |---|6  |---|8  |---|10 |---|12 |---|14 |---|16 |-----●-------●---      |___|   |___|   |___|   |___|   |___|   |___|   |___|   |___|   |___|     |       |          ___     ___     ___     ___     ___     ___     ___     ___             _|_      |         |RY |   |sXd|   |RX |   |Tdg|   | T |   |RX |   |sqX|   |RY |           |CX |     |       --|1  |---|3  |---|5  |---|7  |---|9  |---|11 |---|13 |---|15 |-----------|17 |-----|---      |___|   |___|   |___|   |___|   |___|   |___|   |___|   |___|           |___|     |                                                                                          _|_                                                                                        |CX |     ----------------------------------------------------------------------------------|18 |-                                                                                      |___| 

Next, we create a noise model with some NoiseInstructions. Here we only consider MeasurementNoise.

from quri_parts.circuit.noise import BitFlipNoise, MeasurementNoise, NoiseModelnoise_model = NoiseModel([    MeasurementNoise([BitFlipNoise(0.01)])])

This noise model introduces bit flip error during the measurement.

### Readout error mitigation and peformance​

Here, we explicitly show how to build an estimator that performs readout error mitigation. In this simple example, we will compare the performance of a sampling estimator that runs readout error mitigation with that of noiseless and noisy esimators. We first prepare an operator for this purpose.

from quri_parts.core.operator import Operator, pauli_label, PAULI_IDENTITYop = Operator({    pauli_label("Z0"): 0.25,    pauli_label("Z1 Z2"): 2.0,    pauli_label("X1 X2"): 0.5,    pauli_label("Z1 Y2"): 1.0,    pauli_label("X1 Y2"): 2.0,    PAULI_IDENTITY: 3.0,})

Now, we demonstrate how to construct a readout error mitigation sampling estimator. As we shown in the sampling estimator tutorial, the creation of a sampling estimator would require a concurrent sampler. In the case of a readout error mitigation sampling estimator, we will build a sampler that constructs a filter matrix and performs the mitigation scheme internally. This sampler is created with the create_readout_mitigation_concurrent_sampler function provided in quri_parts.algo, which we show below.

from quri_parts.qulacs.estimator import create_qulacs_vector_estimator, create_qulacs_density_matrix_estimatorfrom quri_parts.qulacs.sampler import create_qulacs_density_matrix_concurrent_samplerfrom quri_parts.core.estimator.sampling import create_sampling_estimatorfrom quri_parts.core.estimator import QuantumEstimatorfrom quri_parts.core.sampling.shots_allocator import create_equipartition_shots_allocatorfrom quri_parts.core.measurement import bitwise_commuting_pauli_measurementfrom quri_parts.core.state import CircuitQuantumState, quantum_statefrom quri_parts.algo.mitigation.readout_mitigation import create_readout_mitigation_concurrent_samplerdef get_readout_mit_estimator(qubit_count: int, shots: int=int(1e6)) -> QuantumEstimator[CircuitQuantumState]:    noisy_concurrent_sampler = create_qulacs_density_matrix_concurrent_sampler(noise_model)    allocator = create_equipartition_shots_allocator()    readout_mit_concurrent_sampler = create_readout_mitigation_concurrent_sampler(        qubit_count,        noisy_concurrent_sampler,  # takes in a noisy concurrent sampler        shots    )    return create_sampling_estimator(        shots,        readout_mit_concurrent_sampler,        bitwise_commuting_pauli_measurement,        allocator    )# Returns the estimatorreadout_mit_estimator = get_readout_mit_estimator(qubit_count=3)

Now, we have a sampling estimator that performs readout error mitigation, which we compare its estimation result with that of a noiseless estimator and a noisy estimator.

exact_estimator = create_qulacs_vector_estimator()noisy_estimator = create_qulacs_density_matrix_estimator(noise_model)state = quantum_state(3, circuit=circuit)print("Noiseless estimate:", exact_estimator(op, state).value)print("Noisy estimate:", noisy_estimator(op, state).value)print("Readout mitigation estimate:", readout_mit_estimator(op, state).value)
#output    Noiseless estimate: (4.999999999999998+0j)    Noisy estimate: (4.920800000000014+0j)    Readout mitigation estimate: 4.998136830346459

We can see that the noise is indeed tamed and we obtain a better estimation of the expectation value compared to noisy estimation.

## Building the readout error mitigation sampler step by step​

Now, we start to explain all the steps necessary to construct a readout error mitigation (concurrent) sampler. This envolves:

• Create a filter matrix from a noisy sampler.
• Apply the filter matrix to sampling counts produced by a noisy sampler.

### Sampling simulation with Qulacs​

First, we create a noisy Sampler and execute the sampling without error mitigation.

from quri_parts.qulacs.sampler import create_qulacs_density_matrix_samplersampler = create_qulacs_density_matrix_sampler(noise_model)counts = sampler(circuit, shots=10000)counts
#output    {0: 4912, 1: 44, 2: 49, 3: 48, 4: 51, 5: 52, 6: 45, 7: 4799}

### Create filter matrix​

We can use create_filter_matrix to create filter matrix.

from quri_parts.algo.mitigation.readout_mitigation import create_filter_matrixfrom quri_parts.qulacs.sampler import create_qulacs_density_matrix_concurrent_samplerconcurernt_sampler = create_qulacs_density_matrix_concurrent_sampler(noise_model)filter_matrix = create_filter_matrix(qubit_count, concurernt_sampler, shots=1000000)filter_matrix.round(5)
#output    array([[ 1.03065e+00, -1.01600e-02, -1.05200e-02,  1.00000e-04,            -1.03600e-02,  9.00000e-05,  9.00000e-05,  0.00000e+00],           [-1.05200e-02,  1.03066e+00,  1.20000e-04, -1.05200e-02,             1.10000e-04, -1.03900e-02, -0.00000e+00,  1.10000e-04],           [-1.02800e-02,  1.00000e-04,  1.03128e+00, -1.03600e-02,             1.10000e-04, -0.00000e+00, -1.04700e-02,  1.00000e-04],           [ 1.00000e-04, -1.05300e-02, -1.04100e-02,  1.03087e+00,            -0.00000e+00,  1.00000e-04,  1.10000e-04, -1.01200e-02],           [-1.01300e-02,  1.10000e-04,  1.00000e-04, -0.00000e+00,             1.03076e+00, -1.06300e-02, -1.04400e-02,  1.10000e-04],           [ 9.00000e-05, -1.02800e-02, -0.00000e+00,  1.10000e-04,            -1.04400e-02,  1.03110e+00,  1.20000e-04, -1.02800e-02],           [ 9.00000e-05, -0.00000e+00, -1.06800e-02,  1.20000e-04,            -1.03100e-02,  1.10000e-04,  1.03115e+00, -1.03900e-02],           [-0.00000e+00,  1.00000e-04,  1.20000e-04, -1.03300e-02,             1.20000e-04, -1.03800e-02, -1.05600e-02,  1.03047e+00]])

Now we can get error-mitigated counts by calling readout_mitigation.
from quri_parts.algo.mitigation.readout_mitigation import readout_mitigationmitigated_counts = readout_mitigation([counts], filter_matrix)next(mitigated_counts)
#output    {0: 5062.934660510149,     4: 0.7909420515811665,     5: 3.290159602996353,     7: 4946.536432906023}
We can also create a ConcurrentSampler that samples from noisy circuit and performs readout error mitigation behind the scenes.
from quri_parts.algo.mitigation.readout_mitigation import (    create_readout_mitigation_concurrent_sampler,    create_readout_mitigation_sampler,)from quri_parts.qulacs.sampler import create_qulacs_density_matrix_sampler# Create a ConcurrentSamplerrem_concurrent_sampler = create_readout_mitigation_concurrent_sampler(    qubit_count, concurernt_sampler, shots=1000000)# You can also create a Sampler# rem_sampler = create_readout_mitigation_sampler(#     qubit_count, concurernt_sampler, shots=1000000# )# mitigated_counts = rem_sampler(circuit, 10000)# print(mitigated_counts)mitigated_counts_concurrent = rem_concurrent_sampler([(circuit, 10000)])next(mitigated_counts_concurrent)
#output    {0: 4993.502510704794,     1: 1.0062540276240308,     2: 5.686739903778277,     3: 2.304943262678012,     4: 6.040353491108938,     6: 0.7900380906298454,     7: 4991.698210158229}