Note

This page was generated from tutorials/noise/3_measurement_error_mitigation.ipynb.

Measurement Error Mitigation¶

Introduction¶

The measurement calibration is used to mitigate measurement errors. The main idea is to prepare all \(2^n\) basis input states and compute the probability of measuring counts in the other basis states. From these calibrations, it is possible to correct the average results of another experiment of interest. This notebook gives examples for how to use the ignis.mitigation.measurement module.

[1]:

# Import general libraries (needed for functions)
import numpy as np
import time

# Import Qiskit classes
import qiskit
from qiskit import QuantumRegister, QuantumCircuit, ClassicalRegister, Aer
from qiskit.providers.aer import noise
from qiskit.tools.visualization import plot_histogram

# Import measurement calibration functions
from qiskit.ignis.mitigation.measurement import (complete_meas_cal, tensored_meas_cal,
                                                 CompleteMeasFitter, TensoredMeasFitter)

3 Qubit Example of the Calibration Matrices¶

Assume that we would like to generate a calibration matrix for the 3 qubits Q2, Q3 and Q4 in a 5-qubit Quantum Register [Q0,Q1,Q2,Q3,Q4].

Since we have 3 qubits, there are \(2^3=8\) possible quantum states.

Generating Measurement Calibration Circuits¶

First, we generate a list of measurement calibration circuits for the full Hilbert space. Each circuit creates a basis state. If there are \(n=3\) qubits, then we get \(2^3=8\) calibration circuits.

The following function complete_meas_cal returns a list meas_calibs of QuantumCircuit objects containing the calibration circuits, and a list state_labels of the calibration state labels.

The input to this function can be given in one of the following three forms:

qubit_list: A list of qubits to perform the measurement correction on, or:
qr (QuantumRegister): A quantum register, or:
cr (ClassicalRegister): A classical register.

In addition, one can provide a string circlabel, which is added at the beginning of the circuit names for unique identification.

For example, in our case, the input is a 5-qubit QuantumRegister containing the qubits Q2,Q3,Q4:

[2]:

# Generate the calibration circuits
qr = qiskit.QuantumRegister(5)
qubit_list = [2,3,4]
meas_calibs, state_labels = complete_meas_cal(qubit_list=qubit_list, qr=qr, circlabel='mcal')

Print the \(2^3=8\) state labels (for the 3 qubits Q2,Q3,Q4):

[3]:

state_labels

[3]:

['000', '001', '010', '011', '100', '101', '110', '111']

Computing the Calibration Matrix¶

If we do not apply any noise, then the calibration matrix is expected to be the \(8 \times 8\) identity matrix.

[4]:

# Execute the calibration circuits without noise
backend = qiskit.Aer.get_backend('qasm_simulator')
job = qiskit.execute(meas_calibs, backend=backend, shots=1000)
cal_results = job.result()

[5]:

# The calibration matrix without noise is the identity matrix
meas_fitter = CompleteMeasFitter(cal_results, state_labels, circlabel='mcal')
print(meas_fitter.cal_matrix)

[[1. 0. 0. 0. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0. 0. 0.]
 [0. 0. 0. 1. 0. 0. 0. 0.]
 [0. 0. 0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 0. 0. 0. 1.]]

Assume that we apply some noise model from Qiskit Aer to the 5 qubits, then the calibration matrix will have most of its mass on the main diagonal, with some additional ‘noise’.

Alternatively, we can execute the calibration circuits using an IBMQ provider.

[6]:

# Generate a noise model for the 5 qubits
noise_model = noise.NoiseModel()
for qi in range(5):
    read_err = noise.errors.readout_error.ReadoutError([[0.9, 0.1],[0.25,0.75]])
    noise_model.add_readout_error(read_err, [qi])

[7]:

# Execute the calibration circuits
backend = qiskit.Aer.get_backend('qasm_simulator')
job = qiskit.execute(meas_calibs, backend=backend, shots=1000, noise_model=noise_model)
cal_results = job.result()

[8]:

# Calculate the calibration matrix with the noise model
meas_fitter = CompleteMeasFitter(cal_results, state_labels, qubit_list=qubit_list, circlabel='mcal')
print(meas_fitter.cal_matrix)

[[0.731 0.223 0.226 0.067 0.198 0.049 0.052 0.017]
 [0.091 0.583 0.021 0.164 0.026 0.176 0.01  0.056]
 [0.067 0.02  0.588 0.167 0.026 0.008 0.191 0.063]
 [0.005 0.057 0.062 0.502 0.001 0.025 0.025 0.137]
 [0.085 0.035 0.022 0.004 0.625 0.18  0.174 0.046]
 [0.01  0.066 0.003 0.025 0.072 0.488 0.022 0.125]
 [0.011 0.006 0.073 0.021 0.049 0.017 0.473 0.147]
 [0.    0.01  0.005 0.05  0.003 0.057 0.053 0.409]]

[9]:

# Plot the calibration matrix
meas_fitter.plot_calibration()

../../_images/tutorials_noise_3_measurement_error_mitigation_17_0.png

Analyzing the Results¶

We would like to compute the total measurement fidelity, and the measurement fidelity for a specific qubit, for example, Q0.

Since the on-diagonal elements of the calibration matrix are the probabilities of measuring state ‘x’ given preparation of state ‘x’, then the trace of this matrix is the average assignment fidelity.

[10]:

# What is the measurement fidelity?
print("Average Measurement Fidelity: %f" % meas_fitter.readout_fidelity())

# What is the measurement fidelity of Q0?
print("Average Measurement Fidelity of Q0: %f" % meas_fitter.readout_fidelity(
    label_list = [['000','001','010','011'],['100','101','110','111']]))

Average Measurement Fidelity: 0.549875
Average Measurement Fidelity of Q0: 0.814250

Applying the Calibration¶

We now perform another experiment and correct the measured results.

Correct Measurement Noise on a 3Q GHZ State¶

As an example, we start with the 3-qubit GHZ state on the qubits Q2,Q3,Q4:

\[\mid GHZ \rangle = \frac{\mid{000} \rangle + \mid{111} \rangle}{\sqrt{2}}\]

[11]:

# Make a 3Q GHZ state
cr = ClassicalRegister(3)
ghz = QuantumCircuit(qr, cr)
ghz.h(qr[2])
ghz.cx(qr[2], qr[3])
ghz.cx(qr[3], qr[4])
ghz.measure(qr[2],cr[0])
ghz.measure(qr[3],cr[1])
ghz.measure(qr[4],cr[2])

[11]:

<qiskit.circuit.instructionset.InstructionSet at 0x7f4544dacb50>

We now execute the calibration circuits (with the noise model above):

[12]:

job = qiskit.execute([ghz], backend=backend, shots=5000, noise_model=noise_model)
results = job.result()

We now compute the results without any error mitigation and with the mitigation, namely after applying the calibration matrix to the results.

There are two fitting methods for applying the calibration (if no method is defined, then ‘least_squares’ is used). - ‘pseudo_inverse’, which is a direct inversion of the calibration matrix, - ‘least_squares’, which constrains to have physical probabilities.

The raw data to be corrected can be given in a number of forms:

Form1: A counts dictionary from results.get_counts,
Form2: A list of counts of length=len(state_labels),
Form3: A list of counts of length=M*len(state_labels) where M is an integer (e.g. for use with the tomography data),
Form4: A qiskit Result (e.g. results as above).

[13]:

# Results without mitigation
raw_counts = results.get_counts()

# Get the filter object
meas_filter = meas_fitter.filter

# Results with mitigation
mitigated_results = meas_filter.apply(results)
mitigated_counts = mitigated_results.get_counts(0)

We can now plot the results with and without error mitigation:

[14]:

from qiskit.tools.visualization import *
plot_histogram([raw_counts, mitigated_counts], legend=['raw', 'mitigated'])

[14]:

../../_images/tutorials_noise_3_measurement_error_mitigation_27_0.png

Applying to a reduced subset of qubits¶

Consider now that we want to correct a 2Q Bell state, but we have the 3Q calibration matrix. We can reduce the matrix and build a new mitigation object.

[15]:

# Make a 2Q Bell state between Q2 and Q4
cr = ClassicalRegister(2)
bell = QuantumCircuit(qr, cr)
bell.h(qr[2])
bell.cx(qr[2], qr[4])
bell.measure(qr[2],cr[0])
bell.measure(qr[4],cr[1])

[15]:

<qiskit.circuit.instructionset.InstructionSet at 0x7f4544af3bd0>

[16]:

job = qiskit.execute([bell], backend=backend, shots=5000, noise_model=noise_model)
results = job.result()

[17]:

#build a fitter from the subset
meas_fitter_sub = meas_fitter.subset_fitter(qubit_sublist=[2,4])

[18]:

#The calibration matrix is now in the space Q2/Q4
meas_fitter_sub.cal_matrix

[18]:

array([[0.806 , 0.2385, 0.2335, 0.0685],
       [0.0895, 0.653 , 0.031 , 0.197 ],
       [0.0955, 0.033 , 0.6605, 0.195 ],
       [0.009 , 0.0755, 0.075 , 0.5395]])

[19]:

# Results without mitigation
raw_counts = results.get_counts()

# Get the filter object
meas_filter_sub = meas_fitter_sub.filter

# Results with mitigation
mitigated_results = meas_filter_sub.apply(results)
mitigated_counts = mitigated_results.get_counts(0)
from qiskit.tools.visualization import *
plot_histogram([raw_counts, mitigated_counts], legend=['raw', 'mitigated'])

[19]:

../../_images/tutorials_noise_3_measurement_error_mitigation_34_0.png

Tensored mitigation¶

The calibration can be simplified if the error is known to be local. By “local error” we mean that the error can be tensored to subsets of qubits. In this case, less than \(2^n\) states are needed for the computation of the calibration matrix.

Assume that the error acts locally on qubit 2 and the pair of qubits 3 and 4. Construct the calibration circuits by using the function tensored_meas_cal. Unlike before we need to explicitly divide the qubit list up into subset regions.

[20]:

# Generate the calibration circuits
qr = qiskit.QuantumRegister(5)
mit_pattern = [[2],[3,4]]
meas_calibs, state_labels = tensored_meas_cal(mit_pattern=mit_pattern, qr=qr, circlabel='mcal')

We now retrieve the names of the generated circuits. Note that in each label (of length 3), the least significant bit corresponds to qubit 2, the middle bit corresponds to qubit 3, and the most significant bit corresponds to qubit 4.

[21]:

for circ in meas_calibs:
    print(circ.name)

mcalcal_000
mcalcal_010
mcalcal_101
mcalcal_111

Let us elaborate on the circuit names. We see that there are only four circuits, instead of eight. The total number of required circuits is \(2^m\) where \(m\) is the number of qubits in the target subset (here \(m=2\)).

Each basis state of qubits 3 and 4 appears exactly once. Only two basis states are required for qubit 2, so these are split equally across the four experiments. For example, state ‘0’ of qubit 2 appears in state labels ‘000’ and ‘010’.

We now execute the calibration circuits on an Aer simulator, using the same noise model as before. This noise is in fact local to qubits 3 and 4 separately, but assume that we don’t know it, and that we only know that it is local for qubit 2.

[22]:

# Generate a noise model for the 5 qubits
noise_model = noise.NoiseModel()
for qi in range(5):
    read_err = noise.errors.readout_error.ReadoutError([[0.9, 0.1],[0.25,0.75]])
    noise_model.add_readout_error(read_err, [qi])

[23]:

# Execute the calibration circuits
backend = qiskit.Aer.get_backend('qasm_simulator')
job = qiskit.execute(meas_calibs, backend=backend, shots=5000, noise_model=noise_model)
cal_results = job.result()

[24]:

meas_fitter = TensoredMeasFitter(cal_results, mit_pattern=mit_pattern)

The fitter provides two calibration matrices. One matrix is for qubit 2, and the other matrix is for qubits 3 and 4.

[25]:

print(meas_fitter.cal_matrices)

[array([[0.9035, 0.2424],
       [0.0965, 0.7576]]), array([[0.8108, 0.2348, 0.2236, 0.065 ],
       [0.0846, 0.6644, 0.0208, 0.187 ],
       [0.094 , 0.0258, 0.6814, 0.184 ],
       [0.0106, 0.075 , 0.0742, 0.564 ]])]

We can look at the readout fidelities of the individual tensored components or qubits within a set:

[26]:

#readout fidelity of Q2
print('Readout fidelity of Q2: %f'%meas_fitter.readout_fidelity(0))

#readout fidelity of Q3/Q4
print('Readout fidelity of Q3/4 space (e.g. mean assignment '
      '\nfidelity of 00,10,01 and 11): %f'%meas_fitter.readout_fidelity(1))

#readout fidelity of Q3
print('Readout fidelity of Q3: %f'%meas_fitter.readout_fidelity(1,[['00','10'],['01','11']]))

Readout fidelity of Q2: 0.830550
Readout fidelity of Q3/4 space (e.g. mean assignment
fidelity of 00,10,01 and 11): 0.680150
Readout fidelity of Q3: 0.825050

Plot the individual calibration matrices:

[27]:

# Plot the calibration matrix
print('Q2 Calibration Matrix')
meas_fitter.plot_calibration(0)
print('Q3/Q4 Calibration Matrix')
meas_fitter.plot_calibration(1)

Q2 Calibration Matrix

../../_images/tutorials_noise_3_measurement_error_mitigation_48_1.png

Q3/Q4 Calibration Matrix

../../_images/tutorials_noise_3_measurement_error_mitigation_48_3.png

[28]:

# Make a 3Q GHZ state
cr = ClassicalRegister(3)
ghz = QuantumCircuit(qr, cr)
ghz.h(qr[2])
ghz.cx(qr[2], qr[3])
ghz.cx(qr[3], qr[4])
ghz.measure(qr[2],cr[0])
ghz.measure(qr[3],cr[1])
ghz.measure(qr[4],cr[2])

[28]:

<qiskit.circuit.instructionset.InstructionSet at 0x7f4544daca10>

We now execute the calibration circuits (with the noise model above):

[29]:

job = qiskit.execute([ghz], backend=backend, shots=5000, noise_model=noise_model)
results = job.result()

[30]:

# Results without mitigation
raw_counts = results.get_counts()

# Get the filter object
meas_filter = meas_fitter.filter

# Results with mitigation
mitigated_results = meas_filter.apply(results)
mitigated_counts = mitigated_results.get_counts(0)

Plot the raw vs corrected state:

[31]:

meas_filter = meas_fitter.filter
mitigated_results = meas_filter.apply(results)
mitigated_counts = mitigated_results.get_counts(0)
plot_histogram([raw_counts, mitigated_counts], legend=['raw', 'mitigated'])

[31]:

../../_images/tutorials_noise_3_measurement_error_mitigation_54_0.png

As a check we should get the same answer if we build the full correction matrix from a tensor product of the subspace calibration matrices:

[32]:

meas_calibs2, state_labels2 = complete_meas_cal([2,3,4])
meas_fitter2 = CompleteMeasFitter(None, state_labels2)
meas_fitter2.cal_matrix = np.kron(meas_fitter.cal_matrices[1],meas_fitter.cal_matrices[0])
meas_filter2 = meas_fitter2.filter
mitigated_results2 = meas_filter2.apply(results)
mitigated_counts2 = mitigated_results2.get_counts(0)
plot_histogram([raw_counts, mitigated_counts2], legend=['raw', 'mitigated'])

[32]:

../../_images/tutorials_noise_3_measurement_error_mitigation_56_0.png

Running Qiskit Algorithms with Measurement Error Mitigation¶

To use measurement error mitigation when running quantum circuits as part of a Qiskit algorithm, we need to include the respective measurement error fitting instance in the QuantumInstance. This object also holds the specifications for the chosen backend.

In the following, we illustrate measurement error mitigation on the example of searching the ground state of a Hamiltonian with VQE.

First, we need to import the libraries that provide backends as well as the classes that are needed to run the algorithm.

[33]:

# Import qiskit functions and libraries
from qiskit import Aer, IBMQ
from qiskit.circuit.library import TwoLocal
from qiskit.utils import QuantumInstance
from qiskit.algorithms import VQE
from qiskit.algorithms.optimizers import COBYLA
from qiskit.opflow import X, Y, Z, I, CX, T, H, S, PrimitiveOp
from qiskit.providers.aer import noise

# Import error mitigation functions
from qiskit.ignis.mitigation.measurement import CompleteMeasFitter

Then, we initialize the instances that are required to execute the algorithm.

[34]:

# Initialize Hamiltonian
h_op =  (-1.0523732 * I^I) + \
                (0.39793742 * I^Z) + \
                (-0.3979374 * Z^I) + \
                (-0.0112801 * Z^Z) + \
                (0.18093119 * X^X)
# Initialize trial state
ansatz = TwoLocal(h_op.num_qubits, ['ry', 'rz'], 'cz', reps=3, entanglement='full')
# Initialize optimizer
optimizer = COBYLA(maxiter=350)

Here, we choose the Aer qasm_simulator as backend and also add a custom noise model. The application of an actual quantum backend provided by IBMQ is outlined in the commented code.

[35]:

# Generate a noise model
noise_model = noise.NoiseModel()
for qi in range(h_op.num_qubits):
    read_err = noise.errors.readout_error.ReadoutError([[0.8, 0.2],[0.1,0.9]])
    noise_model.add_readout_error(read_err, [qi])

# Initialize the backend configuration using measurement error mitigation with a QuantumInstance
qi_noise_model_qasm = QuantumInstance(backend=Aer.get_backend('qasm_simulator'), noise_model=noise_model, shots=1000,
                                      measurement_error_mitigation_cls=CompleteMeasFitter,
                                      measurement_error_mitigation_shots=1000)

# Intialize your TOKEN and provider with
# provider = IBMQ.get_provider(...)
# qi_noise_model_ibmq = QuantumInstance(backend=provider = provider.get_backend(backend_name)), shots=8000,
#                                       measurement_error_mitigation_cls=CompleteMeasFitter, measurement_error_mitigation_shots=8000)
# Initialize algorithm to find the ground state
vqe = VQE(ansatz, optimizer, quantum_instance=qi_noise_model_qasm)

/home/computertreker/git/qiskit/qiskit/.tox/docs/lib/python3.7/site-packages/sympy/core/expr.py:3951: SymPyDeprecationWarning:

expr_free_symbols method has been deprecated since SymPy 1.9. See
https://github.com/sympy/sympy/issues/21494 for more info.

  deprecated_since_version="1.9").warn()

Finally, we can run the algorithm and check the results.

[36]:

# Run the algorithm
result = vqe.compute_minimum_eigenvalue(h_op)
print(result)

{   'aux_operator_eigenvalues': None,
    'cost_function_evals': 120,
    'eigenstate': {   '00': 0.46311154903618346,
                      '01': 0.24524181572589915,
                      '10': 0.6152773988627401,
                      '11': 0.5889124446108024},
    'eigenvalue': (-0.7792974258788177+0j),
    'optimal_parameters': {   ParameterVectorElement(θ[8]): 3.961350502541183,
                              ParameterVectorElement(θ[9]): -1.5423144397173896,
                              ParameterVectorElement(θ[10]): 5.648550809031063,
                              ParameterVectorElement(θ[11]): -3.6202876936457313,
                              ParameterVectorElement(θ[12]): -1.7411729885029326,
                              ParameterVectorElement(θ[13]): -2.9363357564805104,
                              ParameterVectorElement(θ[14]): 2.367864822222402,
                              ParameterVectorElement(θ[15]): 1.7483777745309879,
                              ParameterVectorElement(θ[0]): -3.6228128782494653,
                              ParameterVectorElement(θ[1]): 3.001032564834221,
                              ParameterVectorElement(θ[2]): -0.6789306545344673,
                              ParameterVectorElement(θ[3]): -4.1901530099471485,
                              ParameterVectorElement(θ[4]): 0.5897232487204125,
                              ParameterVectorElement(θ[5]): 0.7672718412646299,
                              ParameterVectorElement(θ[6]): 5.7181850269935035,
                              ParameterVectorElement(θ[7]): -0.4052456346342691},
    'optimal_point': array([-3.62281288,  3.00103256, -0.67893065, -4.19015301,  0.58972325,
        0.76727184,  5.71818503, -0.40524563,  3.9613505 , -1.54231444,
        5.64855081, -3.62028769, -1.74117299, -2.93633576,  2.36786482,
        1.74837777]),
    'optimal_value': -0.7792974258788177,
    'optimizer_evals': None,
    'optimizer_time': 49.09925627708435}

[37]:

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

/home/computertreker/git/qiskit/qiskit/.tox/docs/lib/python3.7/site-packages/qiskit/aqua/__init__.py:86: DeprecationWarning: The package qiskit.aqua is deprecated. It was moved/refactored to qiskit-terra For more information see <https://github.com/Qiskit/qiskit-aqua/blob/main/README.md#migration-guide>
  warn_package('aqua', 'qiskit-terra')

Version Information

Qiskit Software	Version
`qiskit-terra`	0.18.2
`qiskit-aer`	0.8.2
`qiskit-ignis`	0.6.0
`qiskit-ibmq-provider`	0.16.0
`qiskit-aqua`	0.9.5
`qiskit`	0.29.1
`qiskit-nature`	0.2.2
`qiskit-finance`	0.3.0
`qiskit-optimization`	0.2.3
`qiskit-machine-learning`	0.2.1
System information
Python	3.7.12 (default, Nov 22 2021, 14:57:10) [GCC 11.1.0]
OS	Linux
CPUs	32
Memory (Gb)	125.71650314331055
Tue Jan 04 11:09:49 2022 EST

This code is a part of Qiskit

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.

[ ]: