Repository: PennyLaneAI/pennylane

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



Discussion:
andrijapau: Hi @lillian542,

I've opened up this PR as a draft in case you have some early comments. 

Two things I could use some clarification on: 

(1) Any feedback on my approach to benchmarking the method's performance? To me over a second of computation time could warrant some warning from the program. Although, perhaps my experience is a bit biased since I'm used to relatively small quantum system sizes 😁.

(2) Could you confirm the desired output type of `qml.Hermitian.compute_decomposition`? From the example in #6047 it seems like `PauliSentence` would be ideal, but that goes against the expected base class return type of `list[Operator]` causing the `jax-tests` above to fail. I have a fix ready to go which passes the test suite locally but it essentially acknowledges and waives the error in the relevant `ops/functions/conftest.py` file.

Just wanted to double check, thanks!
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/6062?dropdown=coverage&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:
> Project coverage is 99.66%. Comparing base [(`1f3dbc7`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/1f3dbc77ca8e0dcc3f8d085f42c5291c5930060f?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) to head [(`e2c9d6d`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/e2c9d6d235bfd6366da6286c2df8db78b33281c7?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> Report is 1 commits behind head on master.

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


```diff
@@            Coverage Diff             @@
##           master    #6062      +/-   ##
==========================================
- Coverage   99.66%   99.66%   -0.01%     
==========================================
  Files         430      430              
  Lines       41762    41480     -282     
==========================================
- Hits        41623    41340     -283     
- Misses        139      140       +1     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/6062?dropdown=coverage&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).

andrijapau: Great! Looks like the PR is ready to be merged whenever you are.

Thank you all for your help and guidance and I'm excited to be a Pennylane contributor! 😊
Alex-Preciado: Thank you so much, @andrijapau! I appreciate the effort you put into the implementation. I have a few questions about your plot: I noticed you plotted the theoretical time complexity O(n 4^n). Could you explain why this is the expected complexity? Also, how should I interpret the horizontal bars on your plot?
andrijapau: Certainly, @Alex-Preciado!

> Could you explain why this is the expected complexity? 

The short answer is that I read it in the [documentation](https://docs.pennylane.ai/en/stable/code/api/pennylane.pauli_decompose.html). 

<details closed>
<summary>⏰ Open for Time Complexity Analysis ⏰</summary>
<br>
The long answer is that we can analyze the operations in `_generalized_pauli_decompose` that won't take constant time and determine the dominant time complexity.

```python
# Permute by XORing
indices = [qml.math.array(range(shape[0]))] 
for idx in range(shape[0] - 1): # O(2^n)
    indices.append(qml.math.bitwise_xor(indices[-1], (idx + 1) ^ (idx))) 

term_mat = qml.math.cast(
    qml.math.stack(
        [qml.math.gather(matrix[idx], indice) for idx, indice in enumerate(indices)]
    ),
    complex,
) # O((2^n)^2), stacking 2^n arrays of size 2^n.
```
Permuting indices would be overall $\mathcal{O}(4^n)$.

```python
# Perform Hadamard transformation on columns
hadamard_transform_mat = _walsh_hadamard_transform(qml.math.transpose(term_mat)) # O(n 2^n), from documentation

# Account for the phases from Y
phase_mat = qml.math.ones(shape, dtype=complex).reshape((2,) * (2 * num_qubits))
for idx in range(num_qubits):
    index = [slice(None)] * (2 * num_qubits)
    index[idx] = index[idx + num_qubits] = 1
    phase_mat[tuple(index)] *= 1j
phase_mat = qml.math.convert_like(qml.math.reshape(phase_mat, shape), matrix) # O(n)
```
WHT would be dominated by $\mathcal{O}(n 2^n)$.

```python
# c_00 + c_11 -> I; c_00 - c_11 -> Z; c_01 + c_10 -> X; 1j*(c_10 - c_01) -> Y
# https://quantumcomputing.stackexchange.com/a/31790
term_mat = qml.math.transpose(qml.math.multiply(hadamard_transform_mat, phase_mat)) # O((2^n)^2)
```
Element wise multiplication has complexity $\mathcal{O}(4^n)$.

```python
# Obtain the coefficients for each Pauli word
coeffs, obs = [], []
for pauli_rep in product("IXYZ", repeat=num_qubits):
    bit_array = qml.math.array(
        [[(rep in "YZ"), (rep in "XY")] for rep in pauli_rep], dtype=int
    ).T
    coefficient = term_mat[tuple(int("".join(map(str, x)), 2) for x in bit_array)]

    if not is_abstract(matrix) and qml.math.allclose(coefficient, 0):
        continue

    observables = (
        [(o, w) for w, o in zip(wire_order, pauli_rep) if o != I]
        if hide_identity and not all(t == I for t in pauli_rep)
        else [(o, w) for w, o in zip(wire_order, pauli_rep)]
    )
    if observables:
        coeffs.append(coefficient)
        obs.append(observables)

coeffs = qml.math.stack(coeffs)
```
I believe this is where the poor time complexity comes from. We need to (1) loop through all possible `pauli_rep` tensor products ($4^n$, four choices per qubit) and (2) create the `bit_array` by looping through the `pauli_rep` tuple of length $n$ 😱. With indexing being a constant time, the complexity of this section would be $\mathcal{O}(n 4^n)$. 

The total time complexity of this algorithm would be the sum of each of these operations $\mathcal{O}(4^n + n2^n + 4^n + n4^n)$ but would be dominated by $\mathcal{O}(n4^n)$ overall.

</details>

> Also, how should I interpret the horizontal bars on your plot?

The horizontal bars are auto-generated box-and-whisker plots for the performance of each hermitian decomposition test case.

When the benchmarking test is called, the following table is automatically generated,

![benchmark-table.pdf](https://github.com/user-attachments/files/16582371/benchmark-table.pdf)

The reason there is uncertainty in the histogram bars is due to the `Rounds` metric. The benchmarking protocol will execute the test case multiple times to account for any variations in performance. As an illustration, I've outlined a test in white which took on average $\Delta t = 257\ \mu s$ to complete (statistics generated from the 826 rounds of execution).


Alex-Preciado: Awesome! Thank you so much for the detailed explanation, @andrijapau!! I’ll go ahead and update this branch and merge It 😎

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."""