Adding support for multi-qubit observable estimation (#42)

* Adding logic for single observable Tensor, adding integration test; modifying docstring for tests

* Comment

* Remove print statements

* Correct test to have PauliY in there

* Modify test parameter

* Modifying flaky parameters, modifying tolerance value, docstrings

* Adding flaky runs to hermitian test case

* Toggle flaky numbers

* Modify test such that no SWAP erros may be outputted

* Remove print stmt

* Revert tensor obs case QPUDevice.expval

* Add warning to the user, add test for it

* First pass at multi-qubit observables

* Test on wires expval with operator estimation

* Comments

* Test remvoe space

* Linting

* Test docstring

* Changing tests

* Changing PyQvm test settings

* Update pennylane_forest/qpu.py

Co-Authored-By: Josh Izaac <josh146@gmail.com>

* Update pennylane_forest/qpu.py

Co-Authored-By: Josh Izaac <josh146@gmail.com>

* Moving assert condition for multi-qub osbervables into loop condition

* Minor clean-up

* Moving check for Hermitian observables at same level of check for no. of wires

* Multi qubit observables with and without readout errors

* Testing for multi-qubit readout mitigtation

* Removing unnecessary dict

* Cleaning up single-qubit case

* Same preparation program for single and multi qubit observables

* Consolidating and cleaning up logic for expval

* Removing white space

* Keeping same format as master

* No need to separately perform logic for Hermitian observable

* Comments for clarity

* Adding flaky decorator to tests

* Spreading shots budget across 20 diff expts; adding similar test without parametric compilation

Co-authored-by: antalszava <antalszava@gmail.com>
Co-authored-by: Josh Izaac <josh146@gmail.com>

pennylane_forest/qpu.py

@@ -24,17 +24,17 @@
 
 import numpy as np
 from pyquil import get_qc
-from pyquil.operator_estimation import (
-    Experiment,
-    ExperimentSetting,
-    TensorProductState,
-    group_experiments,
-    measure_observables,
-)
+from pyquil.operator_estimation import (Experiment, ExperimentSetting,
+                                        TensorProductState, group_experiments,
+                                        measure_observables)
 from pyquil.paulis import sI, sX, sY, sZ
 from pyquil.quil import Program
 from pyquil.quilbase import Gate
 
+from pennylane.operation import Tensor
+
+from typing import List
+
 from ._version import __version__
 from .qvm import QVMDevice
 
@@ -160,60 +160,71 @@ def __init__(
     def expval(self, observable):
         wires = observable.wires
 
-        if len(wires) == 1 and not self.parametric_compilation:
-            # Single-qubit observable when parametric compilation is turned off
-
-            # identify Experiment Settings for each of the possible single-qubit observables
-            wire = wires[0]
-            qubit = self.wiring[wire]
-            d_expt_settings = {
-                "Identity": [ExperimentSetting(TensorProductState(), sI(qubit))],
-                "PauliX": [ExperimentSetting(TensorProductState(), sX(qubit))],
-                "PauliY": [ExperimentSetting(TensorProductState(), sY(qubit))],
-                "PauliZ": [ExperimentSetting(TensorProductState(), sZ(qubit))],
-                "Hadamard": [
-                    ExperimentSetting(TensorProductState(), float(np.sqrt(1 / 2)) * sX(qubit)),
-                    ExperimentSetting(TensorProductState(), float(np.sqrt(1 / 2)) * sZ(qubit)),
-                ],
-            }
-
-            if observable.name in ["PauliX", "PauliY", "PauliZ", "Identity", "Hadamard"]:
-                # expectation values for single-qubit observables
-
-                prep_prog = Program()
-                for instr in self.program.instructions:
-                    if isinstance(instr, Gate):
-                        # split gate and wires -- assumes 1q and 2q gates
-                        tup_gate_wires = instr.out().split(" ")
-                        gate = tup_gate_wires[0]
-                        str_instr = str(gate)
-                        # map wires to qubits
-                        for w in tup_gate_wires[1:]:
-                            str_instr += f" {int(w)}"
-                        prep_prog += Program(str_instr)
-
-                if self.readout_error is not None:
+        # `measure_observables` called only when parametric compilation is turned off
+        if not self.parametric_compilation:
+
+            # Single-qubit observable
+            if len(wires) == 1:
+
+                # Ensure sensible observable
+                assert observable.name in ["PauliX", "PauliY", "PauliZ", "Identity", "Hadamard"], "Unknown observable"
+
+                # Create appropriate PauliZ operator
+                wire = wires[0]
+                qubit = self.wiring[wire]
+                pauli_obs = sZ(qubit)
+
+            # Multi-qubit observable
+            elif len(wires) > 1 and isinstance(observable, Tensor) and not self.parametric_compilation:
+
+                # All observables are rotated to be measured in the Z-basis, so we just need to
+                # check which wires exist in the observable, map them to physical qubits, and measure
+                # the product of PauliZ operators on those qubits
+                pauli_obs = sI()
+                for wire in observable.wires:
+                    qubit = wire[0]
+                    pauli_obs *= sZ(self.wiring[qubit])
+
+
+            # Program preparing the state in which to measure observable
+            prep_prog = Program()
+            for instr in self.program.instructions:
+                if isinstance(instr, Gate):
+                    # split gate and wires -- assumes 1q and 2q gates
+                    tup_gate_wires = instr.out().split(" ")
+                    gate = tup_gate_wires[0]
+                    str_instr = str(gate)
+                    # map wires to qubits
+                    for w in tup_gate_wires[1:]:
+                        str_instr += f" {int(w)}"
+                    prep_prog += Program(str_instr)
+
+            if self.readout_error is not None:
+                for wire in observable.wires:
+                    if isinstance(wire, int):
+                        qubit = wire
+                    elif isinstance(wire, List):
+                        qubit = wire[0]
                     prep_prog.define_noisy_readout(
-                        qubit, p00=self.readout_error[0], p11=self.readout_error[1]
+                        self.wiring[qubit], p00=self.readout_error[0], p11=self.readout_error[1]
                     )
 
-                # All observables are rotated and can be measured in the PauliZ basis
-                tomo_expt = Experiment(settings=d_expt_settings["PauliZ"], 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,
-                    )
+            # Measure out multi-qubit observable
+            tomo_expt = Experiment(settings=[ExperimentSetting(TensorProductState(), pauli_obs)],
+                                   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])
+            )
+
+            # Return the estimated expectation value
+            return np.sum([expt_result.expectation for expt_result in meas_obs])
 
-            elif observable.name == "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
+        # Calculation of expectation value without using `measure_observables`
         return super().expval(observable)

tests/test_qpu.py

@@ -291,6 +291,82 @@ def circuit_Zmi():
         assert np.allclose(results[:3], 1.0, atol=2e-2)
         assert np.allclose(results[3:], -1.0, atol=2e-2)
 
+    @flaky(max_runs=5, min_passes=3)
+    def test_multi_qub_no_readout_errors(self):
+        """Test the QPU plugin with no readout errors or correction"""
+        device = np.random.choice(VALID_QPU_LATTICES)
+        dev_qpu = qml.device(
+            "forest.qpu",
+            device=device,
+            load_qc=False,
+            symmetrize_readout=None,
+            calibrate_readout=None,
+        )
+
+        @qml.qnode(dev_qpu)
+        def circuit():
+            qml.RY(np.pi / 2, wires=0)
+            qml.RY(np.pi / 3, wires=1)
+            return qml.expval(qml.PauliX(0) @ qml.PauliZ(1))
+
+        num_expts = 50
+        result = 0.0
+        for _ in range(num_expts):
+            result += circuit()
+        result /= num_expts
+
+        assert np.isclose(result, 0.5, atol=2e-2)
+
+    @flaky(max_runs=5, min_passes=3)
+    def test_multi_qub_readout_errors(self):
+        """Test the QPU plugin with readout errors"""
+        device = np.random.choice(VALID_QPU_LATTICES)
+        dev_qpu = qml.device(
+            "forest.qpu",
+            device=device,
+            load_qc=False,
+            shots=10_000,
+            readout_error=[0.9, 0.75],
+            symmetrize_readout=None,
+            calibrate_readout=None,
+            parametric_compilation=False
+        )
+
+        @qml.qnode(dev_qpu)
+        def circuit():
+            qml.RY(np.pi / 2, wires=0)
+            qml.RY(np.pi / 3, wires=1)
+            return qml.expval(qml.PauliX(0) @ qml.PauliZ(1))
+
+        result = circuit()
+
+        assert np.isclose(result, 0.38, atol=2e-2)
+
+    @flaky(max_runs=5, min_passes=3)
+    def test_multi_qub_readout_correction(self):
+        """Test the QPU plugin with readout errors and correction"""
+        device = np.random.choice(VALID_QPU_LATTICES)
+        dev_qpu = qml.device(
+            "forest.qpu",
+            device=device,
+            load_qc=False,
+            shots=10_000,
+            readout_error=[0.9, 0.75],
+            symmetrize_readout='exhaustive',
+            calibrate_readout='plus-eig',
+            parametric_compilation=False
+        )
+
+        @qml.qnode(dev_qpu)
+        def circuit():
+            qml.RY(np.pi / 2, wires=0)
+            qml.RY(np.pi / 3, wires=1)
+            return qml.expval(qml.PauliX(0) @ qml.PauliZ(1))
+
+        result = circuit()
+
+        assert np.isclose(result, 0.5, atol=2e-2)
+
     @flaky(max_runs=5, min_passes=3)
     def test_2q_gate(self):
         """Test that the two qubit gate with the PauliZ observable works correctly.
@@ -374,7 +450,7 @@ def circuit(x):
     @flaky(max_runs=5, min_passes=3)
     @pytest.mark.parametrize("a", np.linspace(-np.pi / 2, 0, 3))
     @pytest.mark.parametrize("b", np.linspace(0, np.pi / 2, 3))
-    def test_2q_gate_pauliz_pauliz_tensor_parametric_compilation_off(self, a, b):
+    def test_2q_circuit_pauliz_pauliz_tensor(self, a, b):
         """Test that the PauliZ tensor PauliZ observable works correctly, when parametric compilation
         is turned off.
 
@@ -390,7 +466,6 @@ def test_2q_gate_pauliz_pauliz_tensor_parametric_compilation_off(self, a, b):
             symmetrize_readout="exhaustive",
             calibrate_readout="plus-eig",
             shots=QVM_SHOTS,
-            parametric_compilation=False,
         )
 
         @qml.qnode(dev_qpu)
@@ -411,6 +486,46 @@ def circuit(x, y):
         # Check that repeated calling of the QNode works correctly
         assert np.allclose(circuit(a, b), analytic_value, atol=2e-2)
 
+    @pytest.mark.parametrize("a", np.linspace(-np.pi / 2, 0, 3))
+    @pytest.mark.parametrize("b", np.linspace(0, np.pi / 2, 3))
+    def test_2q_gate_pauliz_pauliz_tensor_parametric_compilation_off(self, a, b):
+        """Test that the PauliZ tensor PauliZ observable works correctly, when parametric compilation
+        is turned off.
+
+        As the results coming from the qvm are stochastic, a constraint of 3 out of 5 runs was added.
+        """
+
+        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="exhaustive",
+            calibrate_readout="plus-eig",
+            shots=QVM_SHOTS // 20,
+            parametric_compilation=False,
+        )
+
+        @qml.qnode(dev_qpu)
+        def circuit(x, y):
+            qml.RY(x, wires=[0])
+            qml.RY(y, wires=[1])
+            qml.CNOT(wires=[0, 1])
+            return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1))
+
+        analytic_value = (
+            np.cos(a / 2) ** 2 * np.cos(b / 2) ** 2
+            + np.cos(b / 2) ** 2 * np.sin(a / 2) ** 2
+            - np.cos(a / 2) ** 2 * np.sin(b / 2) ** 2
+            - np.sin(a / 2) ** 2 * np.sin(b / 2) ** 2
+        )
+
+        expt = np.mean([circuit(a, b) for _ in range(20)])
+        theory = analytic_value
+
+        assert np.allclose(expt, theory, atol=2e-2)
+
     def test_timeout_set_correctly(self, shots):
         """Test that the timeout attrbiute for the QuantumComputer stored by the QVMDevice
         is set correctly when passing a value as keyword argument"""