Error mitigation#

NISQ devices are noisy due to the imperfect gate operations, unwanted interactions with environment, and measurement errors. Error mitigation is a promising way of reducing the effect of noise.

In this tutorial, we will explore three error mitigation techniques: Reqdout error mitigation, Zero-noise Extrapolation (ZNE), and Clifford Data Regression (CDR).

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. The circuit consists of the identity part and the non-trivial part. Applying the circuit to \(|000\rangle\) results in \(\frac{1}{\sqrt{2}}\left(|000\rangle + |111\rangle\right)\).

[1]:

from quri_parts.circuit import QuantumCircuit

qubit_count = 3

identity_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_circuit
circuit.add_H_gate(0)
circuit.add_CNOT_gate(0, 1)
circuit.add_CNOT_gate(0, 2)

Readout error mitigation#

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.

Create a noise model#

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

[2]:

import quri_parts.circuit.gate_names as gate_names
from quri_parts.circuit.noise import BitFlipNoise, MeasurementNoise, NoiseModel

noises = [MeasurementNoise([BitFlipNoise(0.01)])]
noise_model = NoiseModel(noises)

This noise model introduces bit flip error during the measurement.

Sampling simulation with Qulacs#

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

[3]:

from quri_parts.qulacs.sampler import create_qulacs_density_matrix_sampler

sampler = create_qulacs_density_matrix_sampler(noise_model)
counts = sampler(circuit, shots=10000)
counts

[3]:

{0: 4804, 1: 60, 2: 53, 3: 45, 4: 61, 5: 54, 6: 54, 7: 4869}

Create filter matrix#

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}}\]

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}}.\]

We can use create_filter_matrix to create filter matrix.

[4]:

from quri_parts.algo.mitigation.readout_mitigation import create_filter_matrix
from quri_parts.qulacs.sampler import create_qulacs_density_matrix_concurrent_sampler

concurernt_sampler = create_qulacs_density_matrix_concurrent_sampler(noise_model)

filter_matrix = create_filter_matrix(qubit_count, concurernt_sampler, shots=1000000)
filter_matrix

[4]:

array([[ 1.03081586e+00, -1.03608280e-02, -1.04848053e-02,
         9.09442655e-05, -1.05649343e-02,  9.64611117e-05,
         9.16189081e-05, -2.91328907e-06],
       [-1.04430603e-02,  1.03079763e+00,  1.10411980e-04,
        -1.04499958e-02,  1.20831915e-04, -1.04183144e-02,
        -2.61233451e-07,  1.11304749e-04],
       [-1.04488106e-02,  1.12301168e-04,  1.03087900e+00,
        -1.05218772e-02,  8.23001295e-05, -6.61044448e-07,
        -1.03370308e-02,  1.05931183e-04],
       [ 9.99118580e-05, -1.03935916e-02, -1.03721711e-02,
         1.03099364e+00,  1.43055172e-07,  9.84088448e-05,
         1.22719759e-04, -1.04128304e-02],
       [-1.02422550e-02,  9.88933597e-05,  9.77132129e-05,
         1.37852935e-07,  1.03102844e+00, -1.04927069e-02,
        -1.02896969e-02,  1.12983508e-04],
       [ 1.09499868e-04, -1.03545737e-02,  8.77658401e-08,
         1.11908063e-04, -1.03037655e-02,  1.03090374e+00,
         1.00718846e-04, -1.04533595e-02],
       [ 1.08584290e-04, -1.01949433e-06, -1.03179092e-02,
         1.04857721e-04, -1.04552250e-02,  1.00217507e-04,
         1.03079468e+00, -1.06677857e-02],
       [ 2.67961597e-07,  1.01186145e-04,  8.76776332e-05,
        -1.03296171e-02,  9.22081860e-05, -1.02871410e-02,
        -1.04827521e-02,  1.03120667e+00]])

Execute readout error mitigation#

Now we can get error-mitigated counts by calling readout_mitigation.

[5]:

from quri_parts.algo.mitigation.readout_mitigation import readout_mitigation

mitigated_counts = readout_mitigation([counts], filter_matrix)
next(mitigated_counts)

[5]:

{0: 4950.217659379003,
 1: 11.201708757265243,
 2: 3.9363179303488525,
 4: 13.127927514830965,
 5: 4.058106695396841,
 6: 3.068554432823876,
 7: 5019.37649704677}

Create readout error mitigation sampler#

We can also create a ConcurrentSampler that samples from noisy circuit and performs readout error mitigation behind the scenes.

[6]:

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 ConcurrentSampler
rem_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)

[6]:

{0: 4948.185784886734,
 1: 2.7399812032761752,
 5: 6.073894168192816,
 7: 5054.491959143681}

Zero-noise extrapolation#

ZNE is an error mitigation method that extrapolates the noiseless value from the multiple noise level values. To scale the noise level, whole or part of original circuit is extended. The simplest case is noise level 3.0, where we generate the scaled circuit

\[U_{\text{original}} \rightarrow U_{\text{original}}U^{-1}_{\text{original}}U_{\text{original}}.\]

We can expect that this corresponds to noise level 3.0 since this circuit has three times as many gates as original gate has.

Create a noise model#

First, we create a noise model with some NoiseInstructions. Here we consider BitFlipNoise and DepolarizingNoise.

[7]:

from quri_parts.circuit.noise import (
    BitFlipNoise,
    DepolarizingNoise,
    NoiseModel,
)

noises = [
    BitFlipNoise(
        error_prob=0.01,
        qubit_indices=[],
        target_gates=[],
    ),
    DepolarizingNoise(
        error_prob=0.01,
        qubit_indices=[],
        target_gates=[],
    ),
]
noise_model = NoiseModel(noises)

Estimate the expectation value of an operator with Qulacs#

Next, we create noisy and noiseless estimator and compare the results.

[8]:

from quri_parts.circuit.noise import NoiseModel
from quri_parts.core.operator import Operator, pauli_label
from quri_parts.core.state import GeneralCircuitQuantumState
from quri_parts.qulacs.estimator import (
    create_qulacs_density_matrix_concurrent_estimator,
)

operator = Operator({pauli_label("X0 X1 X2"): 1.0})
state = GeneralCircuitQuantumState(qubit_count, circuit)

noiseless_estimator = create_qulacs_density_matrix_concurrent_estimator(NoiseModel())
estimator = create_qulacs_density_matrix_concurrent_estimator(noise_model)

noiseless_estimate = noiseless_estimator([operator], [state])
print(f"without noise: {noiseless_estimate[0].value}")

estimate = estimator([operator], [state])
print(f"with noise: {estimate[0].value}")

without noise: (1.0000000000000016+0j)
with noise: (0.7725230020220358+0j)

Create scaled circuits#

Next, create scaled circuits. To scale the circuits scaling_circuit_folding can be used. QURI Parts has multiple options for circuit folding. Here we use random folding.

[9]:

from quri_parts.algo.mitigation.zne import (
    create_folding_random,
    scaling_circuit_folding,
)

scale_factors = [1.0, 2.0, 3.0]

random_folding = create_folding_random()
scaled_circuits = [
    scaling_circuit_folding(circuit, scale_factor, random_folding)
    for scale_factor in scale_factors
]

estimates = [
    estimator([operator], [GeneralCircuitQuantumState(qubit_count, scaled_circuit)])[
        0
    ].value
    for scaled_circuit in scaled_circuits
]

for i, scale_factor in enumerate(scale_factors):
    print(
        f"scale factor {scale_factor}: {len(scaled_circuits[i].gates)=}, estimate: {estimates[i]}"
    )

scale factor 1.0: len(scaled_circuits[i].gates)=19, estimate: (0.7725230020220358+0j)
scale factor 2.0: len(scaled_circuits[i].gates)=37, estimate: (0.6202190040024502+0j)
scale factor 3.0: len(scaled_circuits[i].gates)=57, estimate: (0.464645554575118+0j)

Extrapolate zero-noise value#

The last step is extrapolation. QURI Parts has multiple options for extrapolation. Here we use second-order polynomial extrapolation.

[10]:

from quri_parts.algo.mitigation.zne import create_polynomial_extrapolate

poly_extrapolation = create_polynomial_extrapolate(order=2)
exp_vals = poly_extrapolation(scale_factors, estimates)
print(f"mitigated estimate: {exp_vals}")

mitigated estimate: (0.9215575486338767+0j)

Create ZNE estimator#

You can also create a QuantumEstimator that returns the error-mitigated estimate of a given operator and state by performing ZNE behind the scenes.

[11]:

from quri_parts.algo.mitigation.zne import create_zne_estimator

zne_estimator = create_zne_estimator(
    estimator, scale_factors, poly_extrapolation, random_folding
)

mitigated_estimate = zne_estimator(operator, state)
mitigated_estimate.value

[11]:

0.9215575486338767

Clifford data regression#

CDR is an error mitigation method which predicts the noiseless value using training data, which can be generated by exact estimator such as a simulator and a noisy estimator such as a real device. This technique is available even for large systems in the sense that Clifford+T circuits with moderate number of T gates can be simulated efficiently on a classical computer.

Create noise model#

First, we create a noise model with some NoiseInstructions. Here we consider BitFlipNoise and DepolarizingNoise.

[12]:

from quri_parts.circuit.noise import (
    BitFlipNoise,
    DepolarizingNoise,
    NoiseModel,
)

noises = [
    BitFlipNoise(
        error_prob=0.01,
        qubit_indices=[],
        target_gates=[],
    ),
    DepolarizingNoise(
        error_prob=0.01,
        qubit_indices=[],
        target_gates=[],
    ),
]
noise_model = NoiseModel(noises)

Estimate the expectation value of an operator with Qulacs#

Let’s see how noisy result differs from exact one.

[13]:

from quri_parts.circuit.noise import NoiseModel
from quri_parts.core.operator import Operator, pauli_label
from quri_parts.core.state import GeneralCircuitQuantumState
from quri_parts.qulacs.estimator import (
    create_qulacs_density_matrix_concurrent_estimator,
    create_qulacs_vector_concurrent_estimator
)

operator = Operator({pauli_label("X0 X1 X2"): 1.0})
state = GeneralCircuitQuantumState(qubit_count, circuit)

exact_estimator = create_qulacs_vector_concurrent_estimator()
noisy_estimator = create_qulacs_density_matrix_concurrent_estimator(noise_model)

exact_estimate = exact_estimator([operator], [state])
print(f"without noise: {exact_estimate[0].value}")

noisy_estimate = noisy_estimator([operator], [state])
print(f"with noise: {noisy_estimate[0].value}")

without noise: (0.9999999999999996+0j)
with noise: (0.7725230020220358+0j)

Create training data#

Next, we create training circuits and data. Training circuits are generated by replacing randomly chosen non-Clifford gates in the circuit with the closest Clifford gates. The number of non-Clifford gates determines the computational cost for Clifford+T simulator. Here, we generate 8 training circuits and each circuit has 6 non-Clifford gates.

After generating training circuit, we calculate the estimates for exact and noisy estimators.

[14]:

from quri_parts.algo.mitigation.cdr import make_training_circuits

training_circuits = make_training_circuits(circuit, 6, 8)

exact_estimates = [
    exact_estimator(
        [operator], [GeneralCircuitQuantumState(qubit_count, training_circuit)]
    )[0].value
    for training_circuit in training_circuits
]

print("exact estimates:")
print(exact_estimates)

noisy_estimates = [
    noisy_estimator(
        [operator], [GeneralCircuitQuantumState(qubit_count, training_circuit)]
    )[0].value
    for training_circuit in training_circuits
]
print("noisy estimates:")
print(noisy_estimates)

exact estimates:
[(0.056164919507356165+0j), (0.6813385269763016+0j), (0.9574153468466569+0j), (0.7639014342253044+0j), (0.6813385269763014+0j), (0.6813385269763014+0j), (0.7639014342253043+0j), (0.7457052121767197+0j)]
noisy estimates:
[(0.04879702742708523+0j), (0.5069152232916604+0j), (0.7184099115226523+0j), (0.5662337509054978+0j), (0.508024909921587+0j), (0.5069152232916605+0j), (0.5662337509054978+0j), (0.5873135594216767+0j)]

Define the regression function and get the mitigated value#

Finally, we predict noiseless value by regression. QURI Parts has multiple options for regression. Here we perform second-order polynomial regression.

[15]:

from quri_parts.algo.mitigation.cdr import (
    create_polynomial_regression,
)

poly_regression = create_polynomial_regression(order=2)

mitigated_val = poly_regression(
    noisy_estimate[0].value, noisy_estimates, exact_estimates
).real

print(f"mitigated value: {mitigated_val}")

mitigated value: 1.0189823237477753

Create the CDR estimator#

You can also create a CDR estimator that performs CDR behind the scenes. Once it’s created, you can use it as a QuantumEstimator, which accepts Estimatable and ~QuantumState and returns Estimate.

[16]:

from quri_parts.algo.mitigation.cdr import create_cdr_estimator

cdr_estimator = create_cdr_estimator(noisy_estimator, exact_estimator, poly_regression, 10, 0.5)

estimate = cdr_estimator(operator, state)
estimate

[16]:

_Estimate(value=1.0328940345469877, error=nan)