Repository: PennyLaneAI/pennylane
Pull request #4866 (2023-11-21 19:27:24)

Use full dot product for hermitian expvals on `default.qubit`
Solves #4864 
Solves #2269 
Solves #1344 
[sc-50565]

On `default.qubit`, the expectation values of hermitians are now calculated by apply the observable to the state and performing and inner product.  This process is backprop compatible.

This also has substantial performance improvements.

For:

```
H6atasets = qml.data.load("qchem", molname="H6", basis="STO-3G", bondlength=1.16)
h = H6datasets[0].hamiltonian
mat = qml.matrix(h)

n_layers = 1
n_wires = len(h.wires)

shape = qml.StronglyEntanglingLayers.shape(n_wires=n_wires, n_layers=n_layers)
rng = qml.numpy.random.default_rng(seed=42)
params = rng.random(shape)

@qml.qnode(qml.device("default.qubit"))
def circuit(params, h):
    qml.StronglyEntanglingLayers(params, wires=range(n_wires))
    return qml.expval(qml.Hermitian(h, wires=range(n_wires)))
```
![Screenshot 2023-11-21 at 2 44 27 PM](https://github.com/PennyLaneAI/pennylane/assets/6364575/c8d7c992-bdc9-4c03-96c7-895d39001470)



codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/4866?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Comparison is base [(`34aba12`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/34aba1216116dca37fed5afde7cb41d820a042b9?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.65% compared to head [(`19ad58e`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/4866?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.64%.


<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #4866      +/-   ##
==========================================
- Coverage   99.65%   99.64%   -0.01%     
==========================================
  Files         383      383              
  Lines       34518    34267     -251     
==========================================
- Hits        34398    34146     -252     
- Misses        120      121       +1     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/4866?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).


Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index dade66c79..270918021 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -85,6 +85,9 @@
 * Simplified the logic for re-arranging states before returning.
   [(#4817)](https://github.com/PennyLaneAI/pennylane/pull/4817)
 
+* `default.qubit` now calculates the expectation value of Hermitians in a differentiable manner.
+  [(#4866)](https://github.com/PennyLaneAI/pennylane/pull/4866)
+
 * `GlobalPhase` now decomposes to nothing, in case devices do not support global phases.
   [(#4855)](https://github.com/PennyLaneAI/pennylane/pull/4855)
 
diff --git a/pennylane/devices/qubit/measure.py b/pennylane/devices/qubit/measure.py
index 9b72def18..b7b5ab477 100644
--- a/pennylane/devices/qubit/measure.py
+++ b/pennylane/devices/qubit/measure.py
@@ -108,6 +108,28 @@ def csr_dot_products(
     return math.real(math.squeeze(res))
 
 
+def full_dot_products(
+    measurementprocess: ExpectationMP, state: TensorLike, is_state_batched: bool = False
+) -> TensorLike:
+    """Measure the expectation value of an observable using the dot product between full matrix
+    representations.
+
+    Args:
+        measurementprocess (ExpectationMP): measurement process to apply to the state
+        state (TensorLike): the state to measure
+        is_state_batched (bool): whether the state is batched or not
+
+    Returns:
+        TensorLike: the result of the measurement
+    """
+    ket = apply_operation(measurementprocess.obs, state, is_state_batched=is_state_batched)
+    total_indices = len(state.shape) - is_state_batched
+    flattened_state = flatten_state(state, total_indices)
+    flattened_ket = flatten_state(ket, total_indices)
+    dot_product = math.sum(math.conj(flattened_state) * flattened_ket, axis=int(is_state_batched))
+    return math.real(dot_product)
+
+
 def sum_of_terms_method(
     measurementprocess: ExpectationMP, state: TensorLike, is_state_batched: bool = False
 ) -> TensorLike:
@@ -136,6 +158,7 @@ def sum_of_terms_method(
     )
 
 
+# pylint: disable=too-many-return-statements
 def get_measurement_function(
     measurementprocess: MeasurementProcess, state: TensorLike
 ) -> Callable[[MeasurementProcess, TensorLike], TensorLike]:
@@ -157,6 +180,9 @@ def get_measurement_function(
             if measurementprocess.obs.name == "SparseHamiltonian":
                 return csr_dot_products
 
+            if measurementprocess.obs.name == "Hermitian":
+                return full_dot_products
+
             backprop_mode = math.get_interface(state, *measurementprocess.obs.data) != "numpy"
             if isinstance(measurementprocess.obs, Hamiltonian):
                 # need to work out thresholds for when its faster to use "backprop mode" measurements


Test patch:
diff --git a/tests/devices/default_qubit/test_default_qubit.py b/tests/devices/default_qubit/test_default_qubit.py
index a605e5cd5..3d56d299e 100644
--- a/tests/devices/default_qubit/test_default_qubit.py
+++ b/tests/devices/default_qubit/test_default_qubit.py
@@ -785,32 +785,38 @@ class TestExecutingBatches:
         assert qml.math.allclose(g1, g3)
 
 
+@pytest.mark.slow
 class TestSumOfTermsDifferentiability:
     """Basically a copy of the `qubit.simulate` test but using the device instead."""
 
     @staticmethod
-    def f(dev, scale, n_wires=10, offset=0.1, convert_to_hamiltonian=False):
+    def f(dev, scale, n_wires=10, offset=0.1, style="sum"):
         """Execute a quantum script with a large Hamiltonian."""
         ops = [qml.RX(offset + scale * i, wires=i) for i in range(n_wires)]
 
         t1 = 2.5 * qml.prod(*(qml.PauliZ(i) for i in range(n_wires)))
         t2 = 6.2 * qml.prod(*(qml.PauliY(i) for i in range(n_wires)))
         H = t1 + t2
-        if convert_to_hamiltonian:
+        if style == "hamiltonian":
             H = H._pauli_rep.hamiltonian()  # pylint: disable=protected-access
+        elif style == "hermitian":
+            H = qml.Hermitian(H.matrix(), wires=H.wires)
         qs = qml.tape.QuantumScript(ops, [qml.expval(H)])
         return dev.execute(qs)
 
     @staticmethod
-    def f_hashable(scale, n_wires=10, offset=0.1, convert_to_hamiltonian=False):
+    def f_hashable(scale, n_wires=10, offset=0.1, style="sum"):
         """Execute a quantum script with a large Hamiltonian."""
         ops = [qml.RX(offset + scale * i, wires=i) for i in range(n_wires)]
 
         t1 = 2.5 * qml.prod(*(qml.PauliZ(i) for i in range(n_wires)))
         t2 = 6.2 * qml.prod(*(qml.PauliY(i) for i in range(n_wires)))
         H = t1 + t2
-        if convert_to_hamiltonian:
+        if style == "hamiltonian":
             H = H._pauli_rep.hamiltonian()  # pylint: disable=protected-access
+        elif style == "hermitian":
+            H = qml.Hermitian(H.matrix(), wires=H.wires)
+        qs = qml.tape.QuantumScript(ops, [qml.expval(H)])
         qs = qml.tape.QuantumScript(ops, [qml.expval(H)])
         return DefaultQubit().execute(qs)
 
@@ -823,23 +829,23 @@ class TestSumOfTermsDifferentiability:
         return 2.5 * qml.math.prod(cosines) + 6.2 * qml.math.prod(sines)
 
     @pytest.mark.autograd
-    @pytest.mark.parametrize("convert_to_hamiltonian", (True, False))
-    def test_autograd_backprop(self, convert_to_hamiltonian):
+    @pytest.mark.parametrize("style", ("sum", "hamiltonian", "hermitian"))
+    def test_autograd_backprop(self, style):
         """Test that backpropagation derivatives work in autograd with hamiltonians and large sums."""
         dev = DefaultQubit()
         x = qml.numpy.array(0.52)
-        out = self.f(dev, x, convert_to_hamiltonian=convert_to_hamiltonian)
+        out = self.f(dev, x, style=style)
         expected_out = self.expected(x)
         assert qml.math.allclose(out, expected_out)
 
-        g = qml.grad(self.f)(dev, x, convert_to_hamiltonian=convert_to_hamiltonian)
+        g = qml.grad(self.f)(dev, x, style=style)
         expected_g = qml.grad(self.expected)(x)
         assert qml.math.allclose(g, expected_g)
 
     @pytest.mark.jax
     @pytest.mark.parametrize("use_jit", (True, False))
-    @pytest.mark.parametrize("convert_to_hamiltonian", (True, False))
-    def test_jax_backprop(self, convert_to_hamiltonian, use_jit):
+    @pytest.mark.parametrize("style", ("sum", "hamiltonian", "hermitian"))
+    def test_jax_backprop(self, style, use_jit):
         """Test that backpropagation derivatives work with jax with hamiltonians and large sums."""
         import jax
         from jax.config import config
@@ -848,17 +854,17 @@ class TestSumOfTermsDifferentiability:
         x = jax.numpy.array(0.52, dtype=jax.numpy.float64)
         f = jax.jit(self.f_hashable, static_argnums=(1, 2, 3)) if use_jit else self.f_hashable
 
-        out = f(x, convert_to_hamiltonian=convert_to_hamiltonian)
+        out = f(x, style=style)
         expected_out = self.expected(x)
         assert qml.math.allclose(out, expected_out, atol=1e-6)
 
-        g = jax.grad(f)(x, convert_to_hamiltonian=convert_to_hamiltonian)
+        g = jax.grad(f)(x, style=style)
         expected_g = jax.grad(self.expected)(x)
         assert qml.math.allclose(g, expected_g)
 
     @pytest.mark.torch
-    @pytest.mark.parametrize("convert_to_hamiltonian", (True, False))
-    def test_torch_backprop(self, convert_to_hamiltonian):
+    @pytest.mark.parametrize("style", ("sum", "hamiltonian", "hermitian"))
+    def test_torch_backprop(self, style):
         """Test that backpropagation derivatives work with torch with hamiltonians and large sums."""
         import torch
 
@@ -866,7 +872,7 @@ class TestSumOfTermsDifferentiability:
 
         x = torch.tensor(-0.289, requires_grad=True)
         x2 = torch.tensor(-0.289, requires_grad=True)
-        out = self.f(dev, x, convert_to_hamiltonian=convert_to_hamiltonian)
+        out = self.f(dev, x, style=style)
         expected_out = self.expected(x2, like="torch")
         assert qml.math.allclose(out, expected_out)
 
@@ -875,8 +881,8 @@ class TestSumOfTermsDifferentiability:
         assert qml.math.allclose(x.grad, x2.grad)
 
     @pytest.mark.tf
-    @pytest.mark.parametrize("convert_to_hamiltonian", (True, False))
-    def test_tf_backprop(self, convert_to_hamiltonian):
+    @pytest.mark.parametrize("style", ("sum", "hamiltonian", "hermitian"))
+    def test_tf_backprop(self, style):
         """Test that backpropagation derivatives work with tensorflow with hamiltonians and large sums."""
         import tensorflow as tf
 
@@ -885,7 +891,7 @@ class TestSumOfTermsDifferentiability:
         x = tf.Variable(0.5)
 
         with tf.GradientTape() as tape1:
-            out = self.f(dev, x, convert_to_hamiltonian=convert_to_hamiltonian)
+            out = self.f(dev, x, style=style)
 
         with tf.GradientTape() as tape2:
             expected_out = self.expected(x)
diff --git a/tests/devices/qubit/test_measure.py b/tests/devices/qubit/test_measure.py
index 3f9573d18..1348c8250 100644
--- a/tests/devices/qubit/test_measure.py
+++ b/tests/devices/qubit/test_measure.py
@@ -24,6 +24,7 @@ from pennylane.devices.qubit.measure import (
     measure,
     state_diagonalizing_gates,
     csr_dot_products,
+    full_dot_products,
     get_measurement_function,
     sum_of_terms_method,
 )
@@ -62,6 +63,11 @@ class TestMeasurementDispatch:
         gates and allows the measurement to be efficiently computed with them."""
         assert get_measurement_function(m, state=1) is state_diagonalizing_gates
 
+    def test_hermitian_full_dot_product(self):
+        """Test that the expectation value of a hermitian uses the full dot products method."""
+        mp = qml.expval(qml.Hermitian(np.eye(2), wires=0))
+        assert get_measurement_function(mp, state=1) is full_dot_products
+
     def test_hamiltonian_sparse_method(self):
         """Check that the sparse expectation value method is used if the state is numpy."""
         H = qml.Hamiltonian([2], [qml.PauliX(0)])
@@ -122,10 +128,14 @@ class TestMeasurements:
                 ),
                 0.5 * np.sin(0.123) + 2 * np.cos(0.123),
             ),
+            (
+                qml.Hermitian(-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."""
+        """Test that measurements of hamiltonian/ sparse hamiltonian/ hermitians work 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)