Local and correlated readout error characterization experiments (#611)

Added two new experiments LocalReadoutError and CorrelatedReadoutError for characterizing the readout error of devices.

docs/tutorials/readout_mitigation.rst

@@ -1,33 +1,54 @@
 Readout Mitigation
 ==================
 
-Readout mitigation is part of ``qiskit-terra``. The readout mitigators
-can be initalized based on existing backend data, or via a readout
-mitigation experiment.
-
-In a readout mitigation experiment, simple circuits are generated for
-various combinations of “0” and “1” readout values. The results give us
-a matrix describing the probability to obtain a wrong measurement. This
-matrix is used to initialize the readout mitigation object, which is
-given as the result of the expriment.
+Readout errors affect quantum computation during the measurement of the
+qubits in a quantum device. By characterizing the readout errors, it is
+possible to construct a *readout error mitigator* that is used both to
+obtain a more accurate distribution of the outputs, and more accurate
+measurements of expectation value for measurables.
+
+The readout mitigator is generated from an *assignment matrix*: a
+:math:`2^n \times 2^n` matrix :math:`A` such that :math:`A_{y,x}` is the
+probability to observe :math:`y` given the true outcome should be
+:math:`x`. The assignment matrix is used to compute the *mitigation
+matrix* used in the readout error mitigation process itself.
+
+A *Local readout mitigator* works under the assumption that readout
+errors are mostly *local*, meaning readout errors for different qubits
+are independent of each other. In this case, the assignment matrix is
+the tensor product of :math:`n` :math:`2 \times 2` matrices, one for
+each qubit, making it practical to store the assignment matrix in
+implicit form, by storing the individual :math:`2 \times 2` assignment
+matrices. The corresponding class in Qiskit is the `Local readout
+mitigator <https://qiskit.org/documentation/stubs/qiskit.result.LocalReadoutMitigator.html%3E>`__
+in ``qiskit-terra``.
+
+A *Correlated readout mitigator* uses the full :math:`2^n \times 2^n`
+assignment matrix, meaning it can only be used for small values of
+:math:`n`. The corresponding class in Qiskit is the `Correlated readout
+mitigator <https://qiskit.org/documentation/stubs/qiskit.result.CorrelatedReadoutMitigator.html>`__
+in ``qiskit-terra``.
+
+This notebook demonstrates the usage of both the local and correlated
+experiments to generate the corresponding mitigators.
 
 .. jupyter-execute::
 
     import numpy as np
     import matplotlib.pyplot as plt
     from qiskit import QuantumCircuit
     from qiskit.visualization import plot_histogram
-    from qiskit_experiments.library import ReadoutMitigationExperiment
+    from qiskit_experiments.library import LocalReadoutError, CorrelatedReadoutError
     # For simulation
     from qiskit.providers.aer import AerSimulator
     from qiskit.test.mock import FakeParis
-    
+
     from qiskit.result.mitigation.utils import (
         expval_with_stddev,
         str2diag,
         counts_probability_vector
     )
-    
+
     backend = AerSimulator.from_backend(FakeParis())
 
 .. jupyter-execute::
@@ -45,14 +66,15 @@ circuits, one for all “0” and one for all “1” results.
 
 .. jupyter-execute::
 
-    exp = ReadoutMitigationExperiment(qubits)
+    exp = LocalReadoutError(qubits)
     for c in exp.circuits():
         print(c)
 
 
 .. jupyter-execute::
 
-    result = exp.run(backend).block_for_results()
+    exp.analysis.set_options(plot=True)
+    result = exp.run(backend)
     mitigator = result.analysis_results(0).value
 
 The resulting measurement matrix can be illustrated by comparing it to
@@ -62,6 +84,7 @@ the identity.
 
     result.figure(0)
 
+
 Mitigation matrices
 -------------------
 
@@ -145,10 +168,11 @@ a few qubits.
 
     qubits = [0,3]
     num_qubits = len(qubits)
-    exp = ReadoutMitigationExperiment(qubits, method="correlated")
+    exp = CorrelatedReadoutError(qubits)
     for c in exp.circuits():
         print(c)
 
+
 .. jupyter-execute::
 
     import qiskit.tools.jupyter

qiskit_experiments/framework/json.py

@@ -34,6 +34,7 @@
 from qiskit.circuit.library import BlueprintCircuit
 from qiskit.quantum_info import DensityMatrix
 from qiskit.quantum_info.operators.channel.quantum_channel import QuantumChannel
+from qiskit.result import LocalReadoutMitigator, CorrelatedReadoutMitigator
 from qiskit_experiments.version import __version__
 
 
@@ -500,6 +501,14 @@ def default(self, obj: Any) -> Any:  # pylint: disable=arguments-differ
                 "output_dims": obj.output_dims(),
             }
             return _serialize_object(obj, settings=settings)
+        if isinstance(obj, LocalReadoutMitigator):
+            # Temporary handling until serialization is added to terra stable release
+            settings = {"assignment_matrices": obj._assignment_mats, "qubits": obj.qubits}
+            return _serialize_object(obj, settings=settings)
+        if isinstance(obj, CorrelatedReadoutMitigator):
+            # Temporary handling until serialization is added to terra stable release
+            settings = {"assignment_matrix": obj._assignment_mat, "qubits": obj.qubits}
+            return _serialize_object(obj, settings=settings)
         if isinstance(obj, DensityMatrix):
             # Temporary fix for incorrect settings in qiskit-terra
             # See https://github.com/Qiskit/qiskit-terra/pull/7194

qiskit_experiments/library/__init__.py

@@ -72,6 +72,8 @@
     ~characterization.RamseyXY
     ~characterization.FineFrequency
     ~characterization.ReadoutAngle
+    ~characterization.LocalReadoutError
+    ~characterization.CorrelatedReadoutError
     ~characterization.ResonatorSpectroscopy
 
 
@@ -105,7 +107,21 @@ class instance to manage parameters and pulse schedules.
     ~calibration.EFRoughXSXAmplitudeCal
 
 """
-
+from .calibration import (
+    RoughDragCal,
+    FineDragCal,
+    FineXDragCal,
+    FineSXDragCal,
+    RoughAmplitudeCal,
+    RoughXSXAmplitudeCal,
+    EFRoughXSXAmplitudeCal,
+    FineAmplitudeCal,
+    FineXAmplitudeCal,
+    FineSXAmplitudeCal,
+    RoughFrequencyCal,
+    FrequencyCal,
+    FineFrequencyCal,
+)
 from .characterization import (
     T1,
     T2Ramsey,
@@ -128,28 +144,12 @@ class instance to manage parameters and pulse schedules.
     FineFrequency,
     ReadoutAngle,
     ResonatorSpectroscopy,
+    LocalReadoutError,
+    CorrelatedReadoutError,
 )
-
-from .calibration import (
-    RoughDragCal,
-    FineDragCal,
-    FineXDragCal,
-    FineSXDragCal,
-    RoughAmplitudeCal,
-    RoughXSXAmplitudeCal,
-    EFRoughXSXAmplitudeCal,
-    FineAmplitudeCal,
-    FineXAmplitudeCal,
-    FineSXAmplitudeCal,
-    RoughFrequencyCal,
-    FrequencyCal,
-    FineFrequencyCal,
-)
-
 from .randomized_benchmarking import StandardRB, InterleavedRB
 from .tomography import StateTomography, ProcessTomography
 from .quantum_volume import QuantumVolume
-from .mitigation import ReadoutMitigationExperiment
 
 # Experiment Sub-modules
 from . import calibration

qiskit_experiments/library/characterization/__init__.py

@@ -43,6 +43,8 @@
     FineDrag
     FineXDrag
     FineSXDrag
+    LocalReadoutError
+    CorrelatedReadoutError
     ResonatorSpectroscopy
 
 
@@ -65,6 +67,9 @@
     RamseyXYAnalysis
     ReadoutAngleAnalysis
     ResonatorSpectroscopyAnalysis
+    LocalReadoutErrorAnalysis
+    CorrelatedReadoutErrorAnalysis
+
 """
 
 from .analysis import (
@@ -80,6 +85,8 @@
     CrossResonanceHamiltonianAnalysis,
     ReadoutAngleAnalysis,
     ResonatorSpectroscopyAnalysis,
+    LocalReadoutErrorAnalysis,
+    CorrelatedReadoutErrorAnalysis,
 )
 
 from .t1 import T1
@@ -97,4 +104,6 @@
 from .drag import RoughDrag
 from .readout_angle import ReadoutAngle
 from .fine_drag import FineDrag, FineXDrag, FineSXDrag
+from .local_readout_error import LocalReadoutError
+from .correlated_readout_error import CorrelatedReadoutError
 from .resonator_spectroscopy import ResonatorSpectroscopy

qiskit_experiments/library/characterization/analysis/__init__.py

@@ -24,4 +24,6 @@
 from .tphi_analysis import TphiAnalysis
 from .cr_hamiltonian_analysis import CrossResonanceHamiltonianAnalysis
 from .readout_angle_analysis import ReadoutAngleAnalysis
+from .local_readout_error_analysis import LocalReadoutErrorAnalysis
+from .correlated_readout_error_analysis import CorrelatedReadoutErrorAnalysis
 from .resonator_spectroscopy_analysis import ResonatorSpectroscopyAnalysis

qiskit_experiments/library/characterization/analysis/correlated_readout_error_analysis.py

@@ -0,0 +1,124 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# 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.
+"""
+Analysis class to characterize correlated readout error
+"""
+from typing import List, Tuple
+import numpy as np
+import matplotlib.pyplot as plt
+from qiskit.result import CorrelatedReadoutMitigator
+from qiskit_experiments.framework import ExperimentData
+from qiskit_experiments.framework.matplotlib import get_non_gui_ax
+from qiskit_experiments.framework import BaseAnalysis, AnalysisResultData, Options
+
+
+class CorrelatedReadoutErrorAnalysis(BaseAnalysis):
+    r"""
+    Correlated readout error characterization analysis
+
+    # section: overview
+
+        This class generates the full assignment matrix :math:`A` characterizing the
+        readout error for the given qubits from the experiment results
+        and returns the resulting :class:`~qiskit.result.CorrelatedReadoutMitigator`
+
+        :math:`A` is a :math:`2^n\times 2^n` matrix :math:`A` such that :math:`A_{y,x}`
+        is the probability to observe :math:`y` given the true outcome should be :math:`x`.
+
+        In the experiment, for each :math:`x`a circuit is constructed whose expected
+        outcome is :math:`x`. From the observed results on the circuit, the probability for
+        each :math:`y` is determined, and :math:`A_{y,x}` is set accordingly.
+
+        Analysis Results:
+           * "Local Readout Mitigator": The :class:`~qiskit.result.LocalReadoutMitigator`.
+
+        Analysis Figures:
+           * (Optional) A figure of the assignment matrix.
+
+    # section: reference
+        .. ref_arxiv:: 1 2006.14044
+    """
+
+    @classmethod
+    def _default_options(cls) -> Options:
+        """Return default analysis options.
+
+        Analysis Options:
+            plot (bool): Set ``True`` to create figure for fit result.
+            ax (AxesSubplot): Optional. A matplotlib axis object to draw.
+        """
+        options = super()._default_options()
+        options.plot = True
+        options.ax = None
+        return options
+
+    def _run_analysis(
+        self, experiment_data: ExperimentData, **options
+    ) -> Tuple[List[AnalysisResultData], List["matplotlib.figure.Figure"]]:
+        data = experiment_data.data()
+        qubits = experiment_data.metadata["physical_qubits"]
+        labels = [datum["metadata"]["state_label"] for datum in data]
+        matrix = self._generate_matrix(data, labels)
+        result_mitigator = CorrelatedReadoutMitigator(matrix, qubits=qubits)
+        analysis_results = [AnalysisResultData("Correlated Readout Mitigator", result_mitigator)]
+        if self.options.plot:
+            ax = options.get("ax", None)
+            figures = [self._assignment_matrix_visualization(matrix, labels, ax)]
+        else:
+            figures = []
+        return analysis_results, figures
+
+    def _generate_matrix(self, data, labels) -> np.array:
+        list_size = len(labels)
+        matrix = np.zeros([list_size, list_size], dtype=float)
+        # matrix[i][j] is the probability of counting i for expected j
+        for datum in data:
+            expected_outcome = datum["metadata"]["state_label"]
+            j = labels.index(expected_outcome)
+            total_counts = sum(datum["counts"].values())
+            for measured_outcome, count in datum["counts"].items():
+                i = labels.index(measured_outcome)
+                matrix[i][j] = count / total_counts
+        return matrix
+
+    def _assignment_matrix_visualization(
+        self, matrix, labels, ax=None
+    ) -> "matplotlib.figure.Figure":
+        """
+        Plot the assignment matrix (2D color grid plot).
+
+        Args:
+            matrix: assignment matrix to plot
+            ax (matplotlib.axes): settings for the graph
+
+        Returns:
+            The generated plot of the assignment matrix
+
+        Raises:
+            QiskitError: if _cal_matrices was not set.
+
+            ImportError: if matplotlib was not installed.
+
+        """
+
+        if ax is None:
+            ax = get_non_gui_ax()
+        figure = ax.get_figure()
+        ax.matshow(matrix, cmap=plt.cm.binary, clim=[0, 1])
+        ax.set_xlabel("Prepared State")
+        ax.xaxis.set_label_position("top")
+        ax.set_ylabel("Measured State")
+        ax.set_xticks(np.arange(len(labels)))
+        ax.set_yticks(np.arange(len(labels)))
+        ax.set_xticklabels(labels)
+        ax.set_yticklabels(labels)
+        return figure