Repository: PennyLaneAI/pennylane (python)

Issue:
Implement `Hermitian.decomposition`
### Important Note

⚠️ This issue is part of an internal assignment and not meant for external contributors.

### Context

Currently, the `compute_decomposition` method of `Hermitian` is inherited from the `Operator` base class, which returns a `qml.DecompositionUndefinedError`. This means if we try to decompose it, we get:

```
>>> op = qml.X(0) + qml.Y(1) + 2 * qml.X(0) @ qml.Z(3)
>>> hermitian_op = qml.Hermitian(qml.matrix(op), wires=['a', 'b', 'aux'])
>>> hermitian_op.decomposition()
DecompositionUndefinedError
````

We want to implement a decomposition method for `Hermitian`, so that the above example instead gives a sum of Pauli operators, along the lines of:

```
>>> hermitian_op.decomposition()
Y('b') + X('a') + 2.0 * (X('a') @ Z('aux'))
```

### Implementation

The task would involve two main steps:

1) Use the existing `pauli_decompose` method to implement a custom `compute_decomposition` method on `Hermitian`. Note the relationship between `decomposition` and `compute_decomposition` as defined in the abstract base class `Operator`.
2) The `pauli_decompose` method is inefficient for larger matrices. Test performance as a function of matrix size and add a check for a reasonable limit. Above the limit we should raise a warning indicating that the decomposition may be inefficient.

In addition to these code changes:

* A clear example should be included in the `Hermitian.compute_decomposition` method
* Relevant tests should be added to [test_observables.py](https://github.com/PennyLaneAI/pennylane/blob/master/tests/ops/qubit/test_observables.py)
* A Changelog entry should be added to [the changelog](https://github.com/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) under "Improvements"


### Relevant Files:

* [`pennylane/operation.py`]( https://github.com/PennyLaneAI/pennylane/blob/master/pennylane/operation.py). Contains the base `Operation` class, for context.

* [`pennylane/pauli/conversion.py`](https://github.com/PennyLaneAI/pennylane/blob/master/pauli/conversion.py). Contains the `pauli_decompose` function for decomposing a matrix into a sum of Pauli words. 

* [`pennylane/ops/qubit/observables.py`](https://github.com/PennyLaneAI/pennylane/blob/master/pennylane/ops/qubit/observables.py) contains the `Hermitian` class



Pull request:
[Closes #6047] Implements custom `compute_decomposition` method in `Hermitian`
**Context:**

Prior to this fix, the `compute_decomposition` method of `Hermitian` returned a `qml.DecompositionUndefinedError` (due to its `Operator` base class inheritance). The desired fix, as laid out by the issue, was to use the existing `pauli_decompose` method to perform the desired decomposition. 

**Description of the Change:**

Implemented a custom `compute_decomposition` static method  in `Hermitian` to have it call the existing `pauli_decompose` method. 

<ins>Examples</ins>

```python
>>> op = qml.X(0) + qml.Y(1) + 2 * qml.X(0) @ qml.Z(3)
>>> op_matrix = qml.matrix(op)
>>> qml.Hermitian.compute_decomposition(op_matrix, wires=['a', 'b', 'aux'])
[(
      1.0 * (I('a') @ Y('b') @ I('aux'))
    + 1.0 * (X('a') @ I('b') @ I('aux'))
    + 2.0 * (X('a') @ I('b') @ Z('aux'))
)]

>>> op = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
>>> qml.Hermitian.compute_decomposition(op, wires=0)
[(
      0.7071067811865475 * X(0)
    + 0.7071067811865475 * Z(0)
)]
```

 ⏰ *Benchmarking Performance* ⏰

Performance was benchmarked as a function of matrix size (see figure below). I heuristically chose to raise an inefficiency warning if the decomposition time would take longer than a second (open to suggestions 🙂). Any observable with more than seven wires will raise a user warning indicating that the decomposition may be inefficient. 

I've also plotted the theoretical time complexity of the `pauli_decompose` method and the data is in agreement.

![benchmark.pdf](https://github.com/user-attachments/files/16459301/benchmark.pdf)

**Benefits:**

Can now compute the decomposition of any Hermitian matrix.

**Possible Drawbacks:**

N/A

**Related GitHub Issues:**

Closes #6047 


Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 1ce119364..bf26c684a 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -56,6 +56,9 @@
 * Removed `semantic_version` from the list of required packages in PennyLane. 
   [(#5836)](https://github.com/PennyLaneAI/pennylane/pull/5836)
 
+* Added the `compute_decomposition` method for `qml.Hermitian`.
+  [(#6062)](https://github.com/PennyLaneAI/pennylane/pull/6062)
+
 * During experimental program capture, the qnode can now use closure variables.
   [(#6052)](https://github.com/PennyLaneAI/pennylane/pull/6052)
 
@@ -347,6 +350,7 @@ William Maxwell,
 Vincent Michaud-Rioux,
 Anurav Modak,
 Mudit Pandey,
+Andrija Paurevic,
 Erik Schultheis,
 nate stemen,
 David Wierichs,
diff --git a/pennylane/ops/qubit/observables.py b/pennylane/ops/qubit/observables.py
index 9fd4ee43e..64e6f8e52 100644
--- a/pennylane/ops/qubit/observables.py
+++ b/pennylane/ops/qubit/observables.py
@@ -16,6 +16,7 @@ This submodule contains the discrete-variable quantum observables,
 excepting the Pauli gates and Hadamard gate in ``non_parametric_ops.py``.
 """
 
+import warnings
 from collections.abc import Sequence
 from copy import copy
 
@@ -161,6 +162,52 @@ class Hermitian(Observable):
         """
         return self.eigendecomposition["eigval"]
 
+    @staticmethod
+    def compute_decomposition(A, wires):  # pylint: disable=arguments-differ
+        r"""Decomposes a hermitian matrix as a sum of Pauli operators.
+
+        Args:
+            A (array or Sequence): hermitian matrix
+            wires (Iterable[Any], Wires): wires that the operator acts on
+        Returns:
+            list[.Operator]: decomposition of the hermitian matrix
+
+        **Examples**
+
+        >>> op = qml.X(0) + qml.Y(1) + 2 * qml.X(0) @ qml.Z(3)
+        >>> op_matrix = qml.matrix(op)
+        >>> qml.Hermitian.compute_decomposition(op_matrix, wires=['a', 'b', 'aux'])
+        [(
+              1.0 * (I('a') @ Y('b') @ I('aux'))
+            + 1.0 * (X('a') @ I('b') @ I('aux'))
+            + 2.0 * (X('a') @ I('b') @ Z('aux'))
+        )]
+        >>> op = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
+        >>> qml.Hermitian.compute_decomposition(op, wires=0)
+        [(
+              0.7071067811865475 * X(0)
+            + 0.7071067811865475 * Z(0)
+        )]
+
+        """
+        A = qml.math.asarray(A)
+
+        if isinstance(wires, (int, str)):
+            wires = Wires(wires)
+
+        if len(wires) == 0:
+            raise ValueError("Hermitian: wrong number of wires. At least one wire has to be given.")
+        Hermitian._validate_input(A, expected_mx_shape=2 ** len(wires))
+
+        # determined heuristically from test_hermitian_decomposition_performance
+        if len(wires) > 7:
+            warnings.warn(
+                "Decomposition may be inefficient for this large of a matrix.",
+                UserWarning,
+            )
+
+        return [qml.pauli.conversion.pauli_decompose(A, wire_order=wires, pauli=False)]
+
     @staticmethod
     def compute_diagonalizing_gates(eigenvectors, wires):  # pylint: disable=arguments-differ
         r"""Sequence of gates that diagonalize the operator in the computational basis (static method).


Test patch:
diff --git a/tests/ops/functions/conftest.py b/tests/ops/functions/conftest.py
index 93ff0de19..92863eb7a 100644
--- a/tests/ops/functions/conftest.py
+++ b/tests/ops/functions/conftest.py
@@ -79,7 +79,7 @@ _INSTANCES_TO_FAIL = [
     ),
     (
         qml.THermitian(np.eye(3), wires=0),
-        AssertionError,  # qutrit ops fail validation
+        (AssertionError, ValueError),  # qutrit ops fail validation
     ),
     (
         qml.ops.qubit.special_unitary.TmpPauliRot(1.1, "X", [0]),
diff --git a/tests/ops/qubit/test_observables.py b/tests/ops/qubit/test_observables.py
index 448080e8e..51c8996bb 100644
--- a/tests/ops/qubit/test_observables.py
+++ b/tests/ops/qubit/test_observables.py
@@ -65,6 +65,14 @@ EIGVALS_TEST_DATA = [
 
 EIGVALS_TEST_DATA_MULTI_WIRES = [functools.reduce(np.kron, [Y, I, Z])]
 
+DECOMPOSITION_TEST_DATA_MULTI_WIRES = [
+    functools.reduce(np.kron, [X]),
+    functools.reduce(np.kron, [X, Y]),
+    functools.reduce(np.kron, [X, Y, Z]),
+    functools.reduce(np.kron, [X, I, H, Y]),
+    functools.reduce(np.kron, [X, H, Y, Z, X]),
+]
+
 # Testing Projector observable with the basis states.
 BASISSTATEPROJECTOR_EIGVALS_TEST_DATA = [
     (np.array([0, 0])),
@@ -182,7 +190,7 @@ class TestSimpleObservables:
 # Prevents multiple threads from updating Hermitian._eigs at the same time
 @pytest.mark.xdist_group(name="hermitian_cache_group")
 @pytest.mark.usefixtures("tear_down_hermitian")
-class TestHermitian:
+class TestHermitian:  # pylint: disable=too-many-public-methods
     """Test the Hermitian observable"""
 
     def test_preserves_autograd_trainability(self):
@@ -305,6 +313,65 @@ class TestHermitian:
         assert np.allclose(qml.Hermitian._eigs[key]["eigvec"], eigvecs, atol=tol, rtol=0)
         assert len(qml.Hermitian._eigs) == 1
 
+    def test_hermitian_compute_decomposition_wire_exceptions(self):
+        """Tests user errors associated with input exceptions"""
+        single_wire_observable = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
+        double_wire_observable = qml.matrix(qml.X(0) @ qml.X(1))
+
+        # test empty wire list
+        with pytest.raises(ValueError, match="At least one wire"):
+            qml.Hermitian.compute_decomposition(single_wire_observable, wires=[])
+
+        # test wrong number of wires
+        with pytest.raises(ValueError, match="Expected input matrix to have shape"):
+            qml.Hermitian.compute_decomposition(double_wire_observable, wires=[0, 1, 2])
+
+        with pytest.raises(ValueError, match="Expected input matrix to have shape"):
+            qml.Hermitian.compute_decomposition(double_wire_observable, wires="aux")
+
+        # test int as wire type
+        A_decomp = qml.Hermitian.compute_decomposition(single_wire_observable, wires=0)
+        assert np.allclose(A_decomp[0].matrix(), single_wire_observable, rtol=0)
+
+        # test str as wire type
+        A_decomp = qml.Hermitian.compute_decomposition(single_wire_observable, wires="aux")
+        assert np.allclose(A_decomp[0].matrix(), single_wire_observable, rtol=0)
+
+    def test_hermitian_compute_decomposition_inefficiency_warning(self):
+        """Tests user inefficiency warning associated with large matrix decomposition"""
+        num_wires = 8
+        observable = qml.Identity(wires=list(range(num_wires)))
+
+        with pytest.warns(
+            UserWarning, match="Decomposition may be inefficient for this large of a matrix."
+        ):
+            qml.Hermitian.compute_decomposition(
+                qml.matrix(observable), wires=list(range(num_wires))
+            )
+
+    @pytest.mark.parametrize("test_num_wires", list(range(1, 11)))
+    def test_hermitian_compute_decomposition_performance(self, test_num_wires, benchmark):
+        """Tests the performance of the compute_decomposition method of the Hermitian class.
+        Used to determine the minimum matrix dimension to raise an inefficiency warning"""
+        observable = qml.Identity(0)
+        for i in range(test_num_wires):
+            observable = observable @ qml.X(i)
+        benchmark(
+            qml.Hermitian.compute_decomposition, qml.matrix(observable), list(range(test_num_wires))
+        )
+
+    @pytest.mark.parametrize("observable", DECOMPOSITION_TEST_DATA_MULTI_WIRES)
+    def test_hermitian_decomposition(self, observable):
+        """Tests that the compute_decomposition method of the Hermitian class returns the correct result."""
+        num_wires = int(np.log2(len(observable)))
+        A_decomp_static = qml.Hermitian.compute_decomposition(
+            observable, wires=list(range(num_wires))
+        )
+        A_decomp = qml.Hermitian(observable, wires=list(range(num_wires))).decomposition()
+
+        assert np.allclose(A_decomp_static[0].matrix(), observable, rtol=0)
+        assert np.allclose(A_decomp[0].matrix(), observable, rtol=0)
+
     @pytest.mark.parametrize("observable, eigvals, eigvecs", EIGVALS_TEST_DATA)
     def test_hermitian_diagonalizing_gates(self, observable, eigvals, eigvecs, tol, mocker):
         """Tests that the diagonalizing_gates method of the Hermitian class returns the correct results."""
@@ -380,7 +447,7 @@ class TestHermitian:
         assert len(qml.Hermitian._eigs) == 2
 
     @pytest.mark.parametrize("observable, eigvals, eigvecs", EIGVALS_TEST_DATA)
-    def test_hermitian_diagonalizing_gatesi_same_observable_twice(
+    def test_hermitian_diagonalizing_gates_same_observable_twice(
         self, observable, eigvals, eigvecs, tol
     ):
         """Tests that the diagonalizing_gates method of the Hermitian class keeps the same dictionary entries upon multiple calls."""