Repository: PennyLaneAI/pennylane
Pull request #3759 (2023-02-13 23:28:46)

`qubit.simulate` has specialized measurement method for large observables
**Context:**
This PR is furthering the simulate function for the new device API.

**Description of the Change:**
* Added `measure_hamiltonian_expval` method to `devices.qubit.simulate` for measuring `qml.Hamiltonian` or `qml.SparseHamiltonian` expectation values more efficiently.


github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3759?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#3759](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3759?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (2166934) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/023e001cd437e5283a4ffa1b2f73bda38ebf4bc9?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (023e001) will **increase** coverage by `0.00%`.
> The diff coverage is `100.00%`.

```diff
@@           Coverage Diff           @@
##           master    #3759   +/-   ##
=======================================
  Coverage   99.81%   99.81%           
=======================================
  Files         331      331           
  Lines       28982    28994   +12     
=======================================
+ Hits        28927    28939   +12     
  Misses         55       55           
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3759?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/devices/qubit/simulate.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3759?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2RldmljZXMvcXViaXQvc2ltdWxhdGUucHk=) | `100.00% <100.00%> (ø)` | |

Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

trbromley: [SC-31121]

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 8a266cdcf..6db4e1ebe 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -8,6 +8,10 @@
   state vector and a number of shots.
   [(#3720)](https://github.com/PennyLaneAI/pennylane/pull/3720)
 
+* The `simulate` function added to `devices/qubit` now supports measuring expectation values of large observables such as
+  `qml.Hamiltonian`, `qml.SparseHamiltonian`, `qml.Sum`.
+  [(#3759)](https://github.com/PennyLaneAI/pennylane/pull/3759)
+
 <h3>Improvements</h3>
 
 * `AdaptiveOptimizer` is updated to use non-default user-defined qnode arguments.
@@ -27,5 +31,6 @@ This release contains contributions from (in alphabetical order):
 
 Utkarsh Azad
 Soran Jahangiri
+Mudit Pandey
 Matthew Silverman
 Jay Soni
diff --git a/pennylane/devices/qubit/measure.py b/pennylane/devices/qubit/measure.py
new file mode 100644
index 000000000..73c08ac86
--- /dev/null
+++ b/pennylane/devices/qubit/measure.py
@@ -0,0 +1,139 @@
+# Copyright 2018-2023 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Code relevant for performing measurements on a state.
+"""
+from typing import Callable
+
+from scipy.sparse import csr_matrix
+
+from pennylane import math
+from pennylane.ops import Sum
+from pennylane.measurements import StateMeasurement, MeasurementProcess, ExpectationMP
+from pennylane.typing import TensorLike
+from pennylane.wires import Wires
+
+from .apply_operation import apply_operation
+
+
+def state_diagonalizing_gates(
+    measurementprocess: StateMeasurement, state: TensorLike
+) -> TensorLike:
+    """Apply a measurement to state when the measurement process has an observable with diagonalizing gates.
+
+    Args:
+        measurementprocess (StateMeasurement): measurement to apply to the state
+        state (TensorLike): state to apply the measurement to
+
+    Returns:
+        TensorLike: the result of the measurement
+    """
+    for op in measurementprocess.diagonalizing_gates():
+        state = apply_operation(op, state)
+
+    total_indices = len(state.shape)
+    wires = Wires(range(total_indices))
+    return measurementprocess.process_state(math.flatten(state), wires)
+
+
+def state_hamiltonian_expval(measurementprocess: ExpectationMP, state: TensorLike) -> TensorLike:
+    """Measure the expectation value of the state when the measured observable is a ``Hamiltonian`` or
+    ``SparseHamiltonian``.
+
+    Args:
+        measurementprocess (ExpectationMP): measurement process to apply to the state
+        state (TensorLike): the state to measure
+
+    Returns:
+        TensorLike: the result of the measurement
+    """
+    total_wires = len(state.shape)
+    Hmat = measurementprocess.obs.sparse_matrix(wire_order=list(range(total_wires)))
+    _state = math.toarray(state).flatten()
+
+    # Find the expectation value using the <\psi|H|\psi> matrix contraction
+    bra = csr_matrix(math.conj(_state))
+    ket = csr_matrix(_state[..., None])
+    new_ket = csr_matrix.dot(Hmat, ket)
+    res = csr_matrix.dot(bra, new_ket).toarray()[0]
+
+    return math.real(math.squeeze(res))
+
+
+def state_measurement_process(
+    measurementprocess: StateMeasurement, state: TensorLike
+) -> TensorLike:
+    """Dispatcher to `StateMeasurement.process_state` for obtaining measurement outcomes
+
+    Args:
+        measurementprocess (StateMeasurement): measurement to apply to the state
+        state (TensorLike): state to apply the measurement to
+
+    Returns:
+        TensorLike: the result of the measurement
+    """
+    total_indices = len(state.shape)
+    wires = Wires(range(total_indices))
+    return measurementprocess.process_state(math.flatten(state), wires)
+
+
+def get_measurement_function(
+    measurementprocess: MeasurementProcess,
+) -> Callable[[MeasurementProcess, TensorLike], TensorLike]:
+    """Get the appropriate method for performing a measurement.
+
+    Args:
+        measurementprocess (MeasurementProcess): measurement process to apply to the state
+
+    Returns:
+        Callable: function that returns the measurement result
+    """
+    if isinstance(measurementprocess, StateMeasurement):
+        if measurementprocess.obs is None:
+            # no need to apply diagonalizing gates
+            return state_measurement_process
+
+        if isinstance(measurementprocess, ExpectationMP) and measurementprocess.obs.name in (
+            "Hamiltonian",
+            "SparseHamiltonian",
+        ):
+            return state_hamiltonian_expval
+
+        if (
+            isinstance(measurementprocess, ExpectationMP)
+            and isinstance(measurementprocess.obs, Sum)
+            and measurementprocess.obs.has_overlapping_wires
+            and len(measurementprocess.obs.wires) > 7
+        ):
+            # Use tensor contraction for `Sum` expectation values with non-commuting summands
+            # and 8 or more wires as it's faster than using eigenvalues.
+            return state_hamiltonian_expval
+
+        if measurementprocess.obs.has_diagonalizing_gates:
+            return state_diagonalizing_gates
+
+    raise NotImplementedError
+
+
+def measure(measurementprocess: MeasurementProcess, state: TensorLike) -> TensorLike:
+    """Apply a measurement process to a state.
+
+    Args:
+        measurementprocess (MeasurementProcess): measurement process to apply to the state
+        state (TensorLike): the state to measure
+
+    Returns:
+        Tensorlike: the result of the measurement
+    """
+    return get_measurement_function(measurementprocess)(measurementprocess, state)
diff --git a/pennylane/devices/qubit/simulate.py b/pennylane/devices/qubit/simulate.py
index b4737b2a6..e591e55ad 100644
--- a/pennylane/devices/qubit/simulate.py
+++ b/pennylane/devices/qubit/simulate.py
@@ -15,56 +15,10 @@
 
 # pylint: disable=protected-access
 import pennylane as qml
-from pennylane.wires import Wires
-from pennylane.measurements import StateMeasurement, MeasurementProcess
-from pennylane.typing import TensorLike
 
 from .initialize_state import create_initial_state
 from .apply_operation import apply_operation
-
-
-def measure_state_diagonalizing_gates(
-    measurementprocess: StateMeasurement, state: TensorLike
-) -> TensorLike:
-    """Apply a measurement to state when the measurement process has an observable with diagonalizing gates.
-
-    Args:
-        measurementprocess (StateMeasurement): measurement to apply to the state
-        state (TensorLike): state to apply the measurement to
-
-    Returns:
-        TensorLike: the result of the measurement
-    """
-    total_indices = len(state.shape)
-    wires = Wires(range(total_indices))
-
-    for op in measurementprocess.obs.diagonalizing_gates():
-        state = apply_operation(op, state)
-
-    return measurementprocess.process_state(qml.math.flatten(state), wires)
-
-
-def measure(measurementprocess: MeasurementProcess, state: TensorLike) -> TensorLike:
-    """Apply a measurement process to a state.
-
-    Args:
-        measurementprocess (MeasurementProcess): measurement process to apply to the state
-        state (TensorLike): the state to measure
-
-    Returns:
-        Tensorlike: the result of the measurement
-    """
-    if isinstance(measurementprocess, StateMeasurement):
-        if measurementprocess.obs is None:
-            # no need to apply diagonalizing gates
-            total_indices = len(state.shape)
-            wires = Wires(range(total_indices))
-            return measurementprocess.process_state(qml.math.flatten(state), wires)
-
-        if measurementprocess.obs.has_diagonalizing_gates:
-            return measure_state_diagonalizing_gates(measurementprocess, state)
-
-    raise NotImplementedError
+from .measure import measure
 
 
 def simulate(circuit: qml.tape.QuantumScript) -> tuple:


Test patch:
diff --git a/tests/devices/qubit/test_measure.py b/tests/devices/qubit/test_measure.py
new file mode 100644
index 000000000..415b4a374
--- /dev/null
+++ b/tests/devices/qubit/test_measure.py
@@ -0,0 +1,132 @@
+# Copyright 2018-2023 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Unit tests for measure in devices/qubit."""
+
+import pytest
+
+import numpy as np
+from scipy.sparse import csr_matrix
+
+import pennylane as qml
+from pennylane.devices.qubit import simulate
+from pennylane.devices.qubit.measure import (
+    measure,
+    state_diagonalizing_gates,
+    state_hamiltonian_expval,
+    state_measurement_process,
+    get_measurement_function,
+)
+
+
+class TestCurrentlyUnsupportedCases:
+    # pylint: disable=too-few-public-methods
+    def test_sample_based_observable(self):
+        """Test sample-only measurements raise a notimplementedError."""
+
+        state = 0.5 * np.ones((2, 2))
+        with pytest.raises(NotImplementedError):
+            _ = measure(qml.sample(wires=0), state)
+
+
+class TestMeasurements:
+    @pytest.mark.parametrize(
+        "measurement, expected",
+        [
+            (qml.state(), -0.5j * np.ones(4)),
+            (qml.density_matrix(wires=0), 0.5 * np.ones((2, 2))),
+            (qml.probs(wires=[0]), np.array([0.5, 0.5])),
+        ],
+    )
+    def test_state_measurement_no_obs(self, measurement, expected):
+        """Test that state measurements with no observable work as expected."""
+        state = -0.5j * np.ones((2, 2))
+        res = measure(measurement, state)
+
+        assert np.allclose(res, expected)
+
+    @pytest.mark.parametrize(
+        "obs, expected",
+        [
+            (
+                qml.Hamiltonian([-0.5, 2], [qml.PauliY(0), qml.PauliZ(0)]),
+                0.5 * np.sin(0.123) + 2 * np.cos(0.123),
+            ),
+            (
+                qml.SparseHamiltonian(
+                    csr_matrix(-0.5 * qml.PauliY(0).matrix() + 2 * qml.PauliZ(0).matrix()),
+                    wires=[0],
+                ),
+                0.5 * np.sin(0.123) + 2 * np.cos(0.123),
+            ),
+        ],
+    )
+    def test_hamiltonian_expval(self, obs, expected):
+        """Test that the `measure_hamiltonian_expval` function works correctly."""
+        # Create RX(0.123)|0> state
+        state = np.array([np.cos(0.123 / 2), -1j * np.sin(0.123 / 2)])
+        res = measure(qml.expval(obs), state)
+        assert np.allclose(res, expected)
+
+    def test_sum_expval_tensor_contraction(self):
+        """Test that `Sum` expectation values are correct when tensor contraction
+        is used for computation."""
+        summands = (qml.prod(qml.PauliY(i), qml.PauliZ(i + 1)) for i in range(7))
+        obs = qml.sum(*summands)
+        ops = [qml.RX(0.123, wires=i) for i in range(8)]
+        meas = [qml.expval(obs)]
+        qs = qml.tape.QuantumScript(ops, meas)
+
+        res = simulate(qs)
+        expected = 7 * (-np.sin(0.123) * np.cos(0.123))
+        assert np.allclose(res, expected)
+
+    @pytest.mark.parametrize(
+        "obs, expected",
+        [
+            (qml.sum(qml.PauliY(0), qml.PauliZ(0)), -np.sin(0.123) + np.cos(0.123)),
+            (
+                qml.sum(*(qml.PauliZ(i) for i in range(8))),
+                sum(np.sin(i * np.pi / 2 + 0.123) for i in range(8)),
+            ),
+        ],
+    )
+    def test_sum_expval_eigs(self, obs, expected):
+        """Test that `Sum` expectation values are correct when eigenvalues are used
+        for computation."""
+        ops = [qml.RX(i * np.pi / 2 + 0.123, wires=i) for i in range(8)]
+        meas = [qml.expval(obs)]
+        qs = qml.tape.QuantumScript(ops, meas)
+
+        res = simulate(qs)
+        assert np.allclose(res, expected)
+
+    def test_correct_measurement_used(self):
+        """Test that the correct internal function is used for a given measurement process."""
+        # Test a case where state_measurement_process is used
+        mp1 = qml.state()
+        assert get_measurement_function(mp1) == state_measurement_process
+
+        # Test cases where state_diagonalizing_gates is used
+        mp2 = qml.var(qml.PauliZ(0))
+        assert get_measurement_function(mp2) == state_diagonalizing_gates
+        mp3 = qml.expval(qml.sum(qml.PauliZ(0), qml.PauliX(0)))
+        assert get_measurement_function(mp3) == state_diagonalizing_gates
+        mp4 = qml.expval(qml.sum(*(qml.PauliZ(i) for i in range(8))))
+        assert get_measurement_function(mp4) == state_diagonalizing_gates
+
+        # Test cases where state_hamiltonian_expval is used
+        mp5 = qml.expval(qml.Hamiltonian([], []))
+        assert get_measurement_function(mp5) is state_hamiltonian_expval
+        mp6 = qml.expval(qml.sum(*(qml.prod(qml.PauliY(i), qml.PauliZ(i + 1)) for i in range(7))))
+        assert get_measurement_function(mp6) is state_hamiltonian_expval
diff --git a/tests/devices/qubit/test_simulate.py b/tests/devices/qubit/test_simulate.py
index b7160d2f6..717c89df0 100644
--- a/tests/devices/qubit/test_simulate.py
+++ b/tests/devices/qubit/test_simulate.py
@@ -11,7 +11,7 @@
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
-"""Unit tests for create_initial_state in devices/qubit."""
+"""Unit tests for simulate in devices/qubit."""
 
 import pytest
 
@@ -22,14 +22,7 @@ from pennylane.devices.qubit import simulate
 
 
 class TestCurrentlyUnsupportedCases:
-    def test_hamiltonian_observable(self):
-        """Test that measuring hamiltonians gives a NotImplementedError."""
-
-        H = qml.Hamiltonian([2.0], [qml.PauliX(0)])
-        qs = qml.tape.QuantumScript(measurements=[qml.expval(H)])
-        with pytest.raises(NotImplementedError):
-            simulate(qs)
-
+    # pylint: disable=too-few-public-methods
     def test_sample_based_observable(self):
         """Test sample-only measurements raise a notimplementedError."""