Readout Error Mitigation using Operator Estimation (#19)

* added ability to return multi-qubit expectation values to all devices

* cleaned up by adding separate probability method

* typo

* added multimode hermitian expval test to pyqvm

* added suggested changes

* undo conftest change

* Loading as qvm if load_qc is False

* Switching inheritance from QVMDevice to ForestDevice more directly

* Allowing QPUDevice class to be instantiated

* Calculating PauliZ expectation using measure_observables

* Specifying symmetrization/calibration at instantiation of QPU device

* Removing readout error added to program for debugging purposes

* Double quotes for strings, where possible, as in master

* Missed a single quote

* Changing pre_expval to pre_measure

* Changing expectations to observables

* Changing to a currently available lattice

* Missed a lattice

* Changing expectation to observable, and getting rid of pre-rotations

* Directly creating expt settings for all observables, except for 'Hermitian'

* Using operator_estimation for estimating all observables except for 'Hermitian'

* Using generic qpu to prevent future errors

* Initializing a device once within the class

* Adding optional readout error for testing

* Getting rid of __init__ constructor to let pytest catch test suite

* Cleaning up code for QPU class

* Tests with readout errors but no mitigation

* Averaging over several experiments to get closer agreement

* Adding tests for readout error mitigation

* Example notebook demo-ing readout error mitigation

* Filled out notebook

* Removing unused import

* Making QPUDevice a subclass of QVMDevice

* Removing unnecessary error raising

* Re-factoring out some code to avoid duplication between QVM and QPU

* Adding authorship

README.rst

@@ -71,7 +71,7 @@ We also encourage bug reports, suggestions for new features and enhancements, an
 Authors
 =======
 
-`Josh Izaac <https://github.com/josh146>`_, `Keri A. McKiernan <https://github.com/kmckiern>`_
+`Josh Izaac <https://github.com/josh146>`_, `Keri A. McKiernan <https://github.com/kmckiern>`_, `M. Sohaib Alam <https://github.com/msohaibalam>`_
 
 
 Support

pennylane_forest/qpu.py

@@ -26,6 +26,16 @@
 from .qvm import QVMDevice
 from ._version import __version__
 
+import numpy as np
+
+from pyquil import get_qc
+from pyquil.api._quantum_computer import _get_qvm_with_topology
+from pyquil.gates import MEASURE, RESET
+from pyquil.quil import Pragma, Program
+from pyquil.paulis import sI, sX, sY, sZ
+from pyquil.operator_estimation import ExperimentSetting, TensorProductState, TomographyExperiment, measure_observables, group_experiments
+from pyquil.quilbase import Gate
+
 
 class QPUDevice(QVMDevice):
     r"""Forest QPU device for PennyLane.
@@ -40,6 +50,13 @@ class QPUDevice(QVMDevice):
             evaluation when the number of shots is larger than ~1000.
         load_qc (bool): set to False to avoid getting the quantum computing
             device on initialization. This is convenient if not currently connected to the QPU.
+        readout_error (list): specifies the conditional probabilities [p(0|0), p(1|1)], where
+            p(i|j) denotes the prob of reading out i having sampled j; can be set to `None` if no
+            readout errors need to be simulated; can only be set for the QPU-as-a-QVM
+        symmetrize_readout (str): method to perform readout symmetrization, using exhaustive
+            symmetrization by default
+        calibrate_readout (str): method to perform calibration for readout error mitigation, normalizing
+            by the expectation value in the +1-eigenstate of the observable by default
 
     Keyword args:
         forest_url (str): the Forest URL server. Can also be set by
@@ -54,8 +71,16 @@ class QPUDevice(QVMDevice):
     """
     name = "Forest QPU Device"
     short_name = "forest.qpu"
+    observables = {"PauliX", "PauliY", "PauliZ", "Identity", "Hadamard", "Hermitian"}
+
+    def __init__(self, device, *, shots=1024, active_reset=True, load_qc=True, readout_error=None,
+                 symmetrize_readout="exhaustive", calibrate_readout="plus-eig", **kwargs):
+
+        if readout_error is not None and load_qc:
+            raise ValueError("Readout error cannot be set on the physical QPU")
+
+        self.readout_error = readout_error
 
-    def __init__(self, device, *, shots=1024, active_reset=False, load_qc=True, **kwargs):
         self._eigs = {}
 
         if "wires" in kwargs:
@@ -67,11 +92,69 @@ def __init__(self, device, *, shots=1024, active_reset=False, load_qc=True, **kw
         aspen_match = re.match(r"Aspen-\d+-([\d]+)Q", device)
         num_wires = int(aspen_match.groups()[0])
 
-        super(QVMDevice, self).__init__(
-            num_wires, shots, **kwargs
-        )  # pylint: disable=bad-super-call
+        super(QVMDevice, self).__init__(num_wires, shots, **kwargs)
 
         if load_qc:
             self.qc = get_qc(device, as_qvm=False, connection=self.connection)
+        else:
+            self.qc = get_qc(device, as_qvm=True, connection=self.connection)
 
         self.active_reset = active_reset
+        self.symmetrize_readout = symmetrize_readout
+        self.calibrate_readout = calibrate_readout
+
+    def pre_rotations(self, observable, wires):
+        """
+        This is used in `QVM.pre_measure`. Since the pre-rotations are handled
+        by `measure_observables` (see `expval`), we simply don't do anything here
+        """
+        pass
+
+    def expval(self, observable, wires, par):
+        # identify Experiment Settings for each of the possible observables
+        d_expt_settings = {
+            "Identity": [ExperimentSetting(TensorProductState(), sI(0))],
+            "PauliX": [ExperimentSetting(TensorProductState(), sX(0))],
+            "PauliY": [ExperimentSetting(TensorProductState(), sY(0))],
+            "PauliZ": [ExperimentSetting(TensorProductState(), sZ(0))],
+            "Hadamard": [ExperimentSetting(TensorProductState(), float(np.sqrt(1/2)) * sX(0)),
+                         ExperimentSetting(TensorProductState(), float(np.sqrt(1/2)) * sZ(0))]
+        }
+        # expectation values for single-qubit observables
+        if len(wires) == 1:
+
+            if observable in ["PauliX", "PauliY", "PauliZ", "Identity", "Hadamard"]:
+                qubit = self.qc.qubits()[0]
+                prep_prog = Program([instr for instr in self.program if isinstance(instr, Gate)])
+                if self.readout_error is not None:
+                    prep_prog.define_noisy_readout(qubit, p00=self.readout_error[0],
+                                                          p11=self.readout_error[1])
+                tomo_expt = TomographyExperiment(settings=d_expt_settings[observable], program=prep_prog)
+                grouped_tomo_expt = group_experiments(tomo_expt)
+                meas_obs = list(measure_observables(self.qc, grouped_tomo_expt,
+                                                    active_reset=self.active_reset,
+                                                    symmetrize_readout=self.symmetrize_readout,
+                                                    calibrate_readout=self.calibrate_readout))
+                return np.sum([expt_result.expectation for expt_result in meas_obs])
+
+            elif observable == 'Hermitian':
+                # <H> = \sum_i w_i p_i
+                Hkey = tuple(par[0].flatten().tolist())
+                w = self._eigs[Hkey]['eigval']
+                return w[0]*p0 + w[1]*p1
+
+        # Multi-qubit observable
+        # ----------------------
+        # Currently, we only support qml.expval.Hermitian(A, wires),
+        # where A is a 2^N x 2^N matrix acting on N wires.
+        #
+        # Eventually, we will also support tensor products of Pauli
+        # matrices in the PennyLane UI.
+
+        probs = self.probabilities(wires)
+
+        if observable == 'Hermitian':
+            Hkey = tuple(par[0].flatten().tolist())
+            w = self._eigs[Hkey]['eigval']
+            # <A> = \sum_i w_i p_i
+            return w @ probs

pennylane_forest/qvm.py

@@ -112,25 +112,31 @@ def __init__(self, device, *, shots=1024, noisy=False, **kwargs):
 
         self.active_reset = False
 
+    def pre_rotations(self, observable, wires):
+        """Apply pre-rotations in the case of observales other than 'Hermitian'"""
+        if observable == "PauliX":
+            # X = H.Z.H
+            self.apply("Hadamard", wires, [])
+
+        elif observable == "PauliY":
+            # Y = (HS^)^.Z.(HS^) and S^=SZ
+            self.apply("PauliZ", wires, [])
+            self.apply("S", wires, [])
+            self.apply("Hadamard", wires, [])
+
+        elif observable == "Hadamard":
+            # H = Ry(-pi/4)^.Z.Ry(-pi/4)
+            self.apply("RY", wires, [-np.pi / 4])
+
+
     def pre_measure(self):
         """Run the QVM"""
         # pylint: disable=attribute-defined-outside-init
         for e in self.obs_queue:
             wires = e.wires
 
-            if e.name == "PauliX":
-                # X = H.Z.H
-                self.apply("Hadamard", wires, [])
-
-            elif e.name == "PauliY":
-                # Y = (HS^)^.Z.(HS^) and S^=SZ
-                self.apply("PauliZ", wires, [])
-                self.apply("S", wires, [])
-                self.apply("Hadamard", wires, [])
-
-            elif e.name == "Hadamard":
-                # H = Ry(-pi/4)^.Z.Ry(-pi/4)
-                self.apply("RY", wires, [-np.pi / 4])
+            if e.name in ["PauliX", "PauliY", "PauliZ", "Identity", "Hadamard"]:
+                self.pre_rotations(e.name, wires)
 
             elif e.name == "Hermitian":
                 # For arbitrary Hermitian matrix H, let U be the unitary matrix

tests/test_qpu.py

@@ -4,14 +4,16 @@
 import logging
 
 import pytest
-
+import pyquil
 import pennylane as qml
-
+from pennylane import numpy as np
 from conftest import BaseTest
 
 
 log = logging.getLogger(__name__)
 
+VALID_QPU_LATTICES = [qc for qc in pyquil.list_quantum_computers() if "qvm" not in qc]
+
 
 class TestQPUIntegration(BaseTest):
     """Test the wavefunction simulator works correctly from the PennyLane frontend."""
@@ -20,22 +22,142 @@ class TestQPUIntegration(BaseTest):
 
     def test_load_qpu_device(self):
         """Test that the QPU device loads correctly"""
-        dev = qml.device("forest.qpu", device="Aspen-1-2Q-B", load_qc=False)
+        device = [qpu for qpu in VALID_QPU_LATTICES if '2Q' in qpu][0]
+        dev = qml.device("forest.qpu", device=device, load_qc=False)
         self.assertEqual(dev.num_wires, 2)
         self.assertEqual(dev.shots, 1024)
         self.assertEqual(dev.short_name, "forest.qpu")
 
     def test_load_virtual_qpu_device(self):
         """Test that the QPU simulators load correctly"""
-        qml.device("forest.qpu", device="Aspen-1-2Q-B", load_qc=False)
+        device = np.random.choice(VALID_QPU_LATTICES)
+        qml.device("forest.qpu", device=device, load_qc=False)
 
     def test_qpu_args(self):
         """Test that the QPU plugin requires correct arguments"""
+        device = np.random.choice(VALID_QPU_LATTICES)
+
         with pytest.raises(ValueError, match="QPU device does not support a wires parameter"):
-            qml.device("forest.qpu", device="Aspen-1-7Q-B", wires=2)
+            qml.device("forest.qpu", device=device, wires=2)
 
         with pytest.raises(TypeError, match="missing 1 required positional argument"):
             qml.device("forest.qpu")
 
         with pytest.raises(ValueError, match="Number of shots must be a positive integer"):
-            qml.device("forest.qpu", "Aspen-1-7Q-B", shots=0)
+            qml.device("forest.qpu", device=device, shots=0)
+
+        with pytest.raises(ValueError, match="Readout error cannot be set on the physical QPU"):
+            qml.device("forest.qpu", device=device, load_qc=True, readout_error=[0.9, 0.75])
+
+
+class TestQPUBasic(BaseTest):
+    """Unit tests for the QPU (as a QVM)."""
+
+    # pylint: disable=protected-access
+
+    def test_no_readout_correction(self):
+        # need to find a device with qubit 0
+        found_good_device = False
+        while not found_good_device:
+            device = np.random.choice(VALID_QPU_LATTICES)
+            dev_qpu = qml.device('forest.qpu', device=device, load_qc=False, readout_error=[0.9, 0.75],
+                                symmetrize_readout=None, calibrate_readout=None)
+            if 0 in dev_qpu.qc.qubits():
+                found_good_device = True
+
+        qubit = 0
+
+        @qml.qnode(dev_qpu)
+        def circuit_Xpl():
+            qml.RY(np.pi/2, wires=qubit)
+            return qml.expval(qml.PauliX(qubit))
+
+        @qml.qnode(dev_qpu)
+        def circuit_Xmi():
+            qml.RY(-np.pi/2, wires=qubit)
+            return qml.expval(qml.PauliX(qubit))
+
+        @qml.qnode(dev_qpu)
+        def circuit_Ypl():
+            qml.RX(-np.pi/2, wires=qubit)
+            return qml.expval(qml.PauliY(qubit))
+
+        @qml.qnode(dev_qpu)
+        def circuit_Ymi():
+            qml.RX(np.pi/2, wires=qubit)
+            return qml.expval(qml.PauliY(qubit))
+
+        @qml.qnode(dev_qpu)
+        def circuit_Zpl():
+            qml.RX(0.0, wires=qubit)
+            return qml.expval(qml.PauliZ(qubit))
+
+        @qml.qnode(dev_qpu)
+        def circuit_Zmi():
+            qml.RX(np.pi, wires=qubit)
+            return qml.expval(qml.PauliZ(qubit))
+
+        num_expts = 10
+        results_unavged = np.zeros((num_expts, 6))
+
+        for i in range(num_expts):
+            results_unavged[i, :] = [circuit_Xpl(), circuit_Ypl(), circuit_Zpl(),
+                                     circuit_Xmi(), circuit_Ymi(), circuit_Zmi()]
+
+        results = np.mean(results_unavged, axis=0)
+
+        assert np.allclose(results[:3], 0.8, atol=2e-2)
+        assert np.allclose(results[3:], -0.5, atol=2e-2)
+
+    def test_no_readout_correction(self):
+        # need to find a device with qubit 0
+        found_good_device = False
+        while not found_good_device:
+            device = np.random.choice(VALID_QPU_LATTICES)
+            dev_qpu = qml.device('forest.qpu', device=device, load_qc=False, readout_error=[0.9, 0.75])
+            if 0 in dev_qpu.qc.qubits():
+                found_good_device = True
+
+        qubit = 0
+
+        @qml.qnode(dev_qpu)
+        def circuit_Xpl():
+            qml.RY(np.pi/2, wires=qubit)
+            return qml.expval(qml.PauliX(qubit))
+
+        @qml.qnode(dev_qpu)
+        def circuit_Xmi():
+            qml.RY(-np.pi/2, wires=qubit)
+            return qml.expval(qml.PauliX(qubit))
+
+        @qml.qnode(dev_qpu)
+        def circuit_Ypl():
+            qml.RX(-np.pi/2, wires=qubit)
+            return qml.expval(qml.PauliY(qubit))
+
+        @qml.qnode(dev_qpu)
+        def circuit_Ymi():
+            qml.RX(np.pi/2, wires=qubit)
+            return qml.expval(qml.PauliY(qubit))
+
+        @qml.qnode(dev_qpu)
+        def circuit_Zpl():
+            qml.RX(0.0, wires=qubit)
+            return qml.expval(qml.PauliZ(qubit))
+
+        @qml.qnode(dev_qpu)
+        def circuit_Zmi():
+            qml.RX(np.pi, wires=qubit)
+            return qml.expval(qml.PauliZ(qubit))
+
+        num_expts = 10
+        results_unavged = np.zeros((num_expts, 6))
+
+        for i in range(num_expts):
+            results_unavged[i, :] = [circuit_Xpl(), circuit_Ypl(), circuit_Zpl(),
+                                     circuit_Xmi(), circuit_Ymi(), circuit_Zmi()]
+
+        results = np.mean(results_unavged, axis=0)
+
+        assert np.allclose(results[:3], 1.0, atol=2e-2)
+        assert np.allclose(results[3:], -1.0, atol=2e-2)