Repository: PennyLaneAI/pennylane
Pull request #3384 (2022-11-17 21:51:12)

Update `decompose` utility function and move into Pauli module (PauliUtility #1)
**Context:**
`decompose` is a fairly un-used utility function which takes a general hermitian matrix, decomposes it into a linear combination of pauli words and produces an instance of `Hamiltonian` representing it. This function is moved to the pauli module and updated to support the new `PauliSentence` representation.

**Description of the Change:**
- Move `decompose` from `pennylane/utils.py` --> `pennylane/pauli/conversion.py`
- Move the associated tests 
- Add a keyword argument `pauli=False` which when true produces an instance of `PauliSentence`
- Fixed (independent) bug where calling for wires of an empty  PauliWord caused errors in the to_mat method for PauliSentence 

**Benefits:**
- Reviving old functionality!

codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#3384](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (bc62dcb) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/883dad06c6047b1e2c795e84265c4329588623eb?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (883dad0) will **decrease** coverage by `0.00%`.
> The diff coverage is `100.00%`.

```diff
@@            Coverage Diff             @@
##           master    #3384      +/-   ##
==========================================
- Coverage   99.70%   99.69%   -0.01%     
==========================================
  Files         299      299              
  Lines       26214    26159      -55     
==========================================
- Hits        26136    26080      -56     
- Misses         78       79       +1     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/queuing.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3F1ZXVpbmcucHk=) | `100.00% <ø> (ø)` | |
| [pennylane/tape/tape.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RhcGUvdGFwZS5weQ==) | `97.50% <ø> (-0.54%)` | :arrow_down: |
| [pennylane/transforms/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvX19pbml0X18ucHk=) | `100.00% <ø> (ø)` | |
| [pennylane/transforms/qcut.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvcWN1dC5weQ==) | `100.00% <ø> (ø)` | |
| [pennylane/utils.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3V0aWxzLnB5) | `100.00% <ø> (ø)` | |
| [pennylane/ops/functions/generator.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9mdW5jdGlvbnMvZ2VuZXJhdG9yLnB5) | `100.00% <100.00%> (ø)` | |
| [pennylane/ops/op\_math/exp.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9vcF9tYXRoL2V4cC5weQ==) | `100.00% <100.00%> (ø)` | |
| [pennylane/pauli/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3BhdWxpL19faW5pdF9fLnB5) | `100.00% <100.00%> (ø)` | |
| [pennylane/pauli/conversion.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3BhdWxpL2NvbnZlcnNpb24ucHk=) | `100.00% <100.00%> (ø)` | |
| [pennylane/pauli/pauli\_arithmetic.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3BhdWxpL3BhdWxpX2FyaXRobWV0aWMucHk=) | `100.00% <100.00%> (ø)` | |
| ... and [3 more](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3384/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | |

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)

josh146: Minor comment, but `qml.utils.decompose_hamiltonian` is actually quite a relatively used function from chatting to users. The main problem is that most users could not find it due to it not being top level nor listed in https://docs.pennylane.ai/en/latest/introduction/operations.html, so we get a _lot_ of questions and end up directing users here.

I would also say that 'Pauli decomposition' is the main term users would know it by/try to search for it, so in this particular case I have the following suggestions:

1. Name/location: either
   - We import it top level and name it `qml.pauli_decomposition`/`pauli_decompose`, or
   - We use namespacing, and document it as `qml.pauli.decomposition`/`qml.pauli.decompose`
   (not sure what tense is best?)

2. We should document it in one of the quickstarts, likely operator. It doesn't really fit into existing categories, since this is a rare case of a `matrix` -> `operator` function. However we may have more of these in the future.

AlbertMitjans: Why do we want to allow empty `PauliWord`s?
Jaybsoni: @AlbertMitjans, 

> Why do we want to allow empty `PauliWord`s?

In order to be more efficient, we don't keep track of wires on which only the Identity is applied, So a PauliWord which acts on 2 out of 3 wires is represented as `pw1 = PauliWord({0: X, 2: Z})` and the 1 wire is ignored. Now when we multiply PauliWords, it could be the case that we end up with a Pauli word that has only Identities acting on it (eg. `pw1 * pw1 == PauliWord({0: I, 2:I})  # same as PauliWord({})`. This is also important when we are combining identity terms in a pauli sentence, if we allowed for multiple PauliWords all representing the identity (each a separate key), we would need additional logic to group them. Here we can just map all instances of the identity to the same key: PauliWord({}). 
Jaybsoni: @josh146, I decided to go with option 1. `pauli_decompose(matrix, pauli=False, wire_order=None)` is now top level imported 👍🏼  
Jaybsoni: Not sure why exactly the docs are failing, my suspicion is that it is due to the new top level import.
AlbertMitjans: > Not sure why exactly the docs are failing, my suspicion is that it is due to the new top level import.

There is one warning raised when building the docs:
```
WARNING: Tried to skip objects {'matmul'} in module pennylane.utils, but they were not present.  Ignoring.
```

I believe this is raised because we removed the `matmul` import from `pennylane.utils`, but we are still skipping it in `doc/code/qml_utils.rst`:
```
:skip: Iterable, matmul
```

Removing the `matmul` should fix it.
Jaybsoni: Any recommendations for the name `pauli_decompose`
Jaybsoni: We have decided to keep the UI as is instead of moving this method into a class method for now. Ready to be re-reviewed ! 

Patch:
diff --git a/doc/code/qml_utils.rst b/doc/code/qml_utils.rst
index c5441e595..d4c0ae924 100644
--- a/doc/code/qml_utils.rst
+++ b/doc/code/qml_utils.rst
@@ -12,4 +12,4 @@ qml.utils
     :no-heading:
     :include-all-objects:
     :no-inheritance-diagram:
-    :skip: Iterable, matmul
+    :skip: Iterable
diff --git a/doc/development/deprecations.rst b/doc/development/deprecations.rst
index 43b571e5d..134089962 100644
--- a/doc/development/deprecations.rst
+++ b/doc/development/deprecations.rst
@@ -103,6 +103,12 @@ Pending deprecations
 Completed deprecation cycles
 ----------------------------
 
+* The ``qml.utils.decompose_hamiltonian()`` method is removed. Please
+  use ``qml.pauli_decompose()``.
+
+  - Still accessible in v0.27
+  - Removed in v0.28
+
 * ``qml.tape.get_active_tape`` is deprecated. Please use ``qml.QueuingManager.active_context()`` instead.
 
   - Deprecated in v0.27
diff --git a/doc/introduction/operations.rst b/doc/introduction/operations.rst
index 3b69e7a04..2ae55d828 100644
--- a/doc/introduction/operations.rst
+++ b/doc/introduction/operations.rst
@@ -114,6 +114,24 @@ array([[ 0.92387953+0.j,  0.+0.j ,  0.-0.38268343j,  0.+0.j],
 [ 0.+0.j,  0.+0.38268343j,  0.+0.j,  -0.92387953+0.j]])
 
 
+Matrix to Operator functions
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. autosummary::
+
+    ~pennylane.pauli_decompose
+
+These functions take the matrix representation of an operator and return
+the equivalent native PennyLane operator.
+
+>>> mat = np.array([[1, 1], [1, -1]])
+>>> h = qml.pauli_decompose(mat)
+>>> type(h)
+<class 'pennylane.ops.qubit.hamiltonian.Hamiltonian'>
+>>> print(h)
+(1.0) [X0]
++ (1.0) [Z0]
+
 .. _intro_ref_ops_qubit:
 
 Qubit operators
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 59a54ebd9..fa637f2b7 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -93,6 +93,22 @@
 * New basis sets, `6-311g` and `CC-PVDZ`, are added to the qchem basis set repo.
   [#3279](https://github.com/PennyLaneAI/pennylane/pull/3279)
 
+* Added a `pauli_decompose()` which takes a hermitian matrix and decomposes it in the 
+  Pauli basis, returning it either as a `Hamiltonian` or `PauliSentence` instance.
+  [(#3384)](https://github.com/PennyLaneAI/pennylane/pull/3384)
+
+  ```pycon
+  >>> mat = np.array([[1, 1], [1, -1]])
+  >>> h = qml.pauli_decompose(mat)
+  >>> print(h)
+    (1.0) [X0]
+  + (1.0) [Z0]
+  >>> ps = qml.pauli_decompose(mat, pauli=True, wire_order=["a"])
+  >>> print(ps)
+  1.0 * X(a)
+  + 1.0 * Z(a)
+  ```
+
 * New `pauli_sentence()` function which takes native `Operator` or `Hamiltonian`
   instances representing a linear combination of Pauli words and returns
   the equivalent `PauliSentence`.
@@ -392,6 +408,9 @@
 * Python 3.7 support is no longer maintained.
   [(#3276)](https://github.com/PennyLaneAI/pennylane/pull/3276)
 
+* Removed `qml.utils.decompose_hamiltonian()`, please use `qml.pauli_decompose()` instead.
+  [(#3384)](https://github.com/PennyLaneAI/pennylane/pull/3384)
+
 * Instead of having an `OrderedDict` attribute called `_queue`, `AnnotatedQueue` now inherits from
   `OrderedDict` and encapsulates the queue. Consequentially, this also applies to the `QuantumTape`
   class which inherits from `AnnotatedQueue`.
@@ -518,6 +537,7 @@ Christina Lee
 Albert Mitjans Coma
 Romain Moyard
 Matthew Silverman
+Jay Soni
 Antal Száva
 David Wierichs
 Moritz Willmann
diff --git a/pennylane/__init__.py b/pennylane/__init__.py
index 91a2de096..fc4b4408b 100644
--- a/pennylane/__init__.py
+++ b/pennylane/__init__.py
@@ -33,7 +33,8 @@ import pennylane.math
 import pennylane.operation
 import pennylane.qnn
 import pennylane.templates
-from pennylane import pauli
+import pennylane.pauli
+from pennylane.pauli import pauli_decompose
 import pennylane.resource
 import pennylane.qchem
 from pennylane.qchem import taper, symmetry_generators, paulix_ops, taper_operation, import_operator
diff --git a/pennylane/ops/functions/generator.py b/pennylane/ops/functions/generator.py
index 431c8a4f5..089ff68c7 100644
--- a/pennylane/ops/functions/generator.py
+++ b/pennylane/ops/functions/generator.py
@@ -58,8 +58,7 @@ def _generator_hamiltonian(gen, op):
         elif isinstance(gen, qml.SparseHamiltonian):
             mat = gen.parameters[0].toarray()
 
-        coeffs, obs = qml.utils.decompose_hamiltonian(mat, wire_order=wires, hide_identity=True)
-        H = qml.Hamiltonian(coeffs, obs)
+        H = qml.pauli_decompose(mat, wire_order=wires, hide_identity=True)
 
     elif isinstance(gen, qml.operation.Observable):
         H = 1.0 * gen
diff --git a/pennylane/pauli/__init__.py b/pennylane/pauli/__init__.py
index 776747684..0a1b4715e 100644
--- a/pennylane/pauli/__init__.py
+++ b/pennylane/pauli/__init__.py
@@ -39,6 +39,7 @@ from .utils import (
 )
 
 from .conversion import (
+    pauli_decompose,
     pauli_sentence,
 )
 
diff --git a/pennylane/pauli/conversion.py b/pennylane/pauli/conversion.py
index 42d761237..8e054f9ad 100644
--- a/pennylane/pauli/conversion.py
+++ b/pennylane/pauli/conversion.py
@@ -14,13 +14,133 @@
 """
 Utility functions to convert between ``~.PauliSentence`` and other PennyLane operators.
 """
+from typing import Union
+from operator import matmul
+from itertools import product
 from functools import reduce, singledispatch
 
+import numpy as np
+
 from pennylane.operation import Tensor
 from pennylane.ops import Hamiltonian, Identity, PauliX, PauliY, PauliZ, Prod, SProd, Sum
 
-from .pauli_arithmetic import PauliWord, PauliSentence, X, Y, Z
 from .utils import is_pauli_word
+from .pauli_arithmetic import PauliWord, PauliSentence, I, X, Y, Z, mat_map, op_map
+
+
+def pauli_decompose(
+    H, hide_identity=False, wire_order=None, pauli=False
+) -> Union[Hamiltonian, PauliSentence]:
+    r"""Decomposes a Hermitian matrix into a linear combination of Pauli operators.
+
+    Args:
+        H (array[complex]): a Hermitian matrix of dimension :math:`2^n\times 2^n`.
+        hide_identity (bool): does not include the Identity observable within
+            the tensor products of the decomposition if ``True``.
+        wire_order (list[Union[int, str]]): the ordered list of wires with respect
+            to which the operator is represented as a matrix.
+        pauli (bool): return a PauliSentence instance if ``True``.
+
+    Returns:
+        Union[~.Hamiltonian, ~.PauliSentence]: the matrix decomposed as a linear combination
+            of Pauli operators, either as a :class:``~.Hamiltonian`` or :class:``~.PauliSentence`` instance.
+
+    **Example:**
+
+    We can use this function to compute the Pauli operator decomposition of an arbitrary Hermitian
+    matrix:
+
+    >>> A = np.array(
+    ... [[-2, -2+1j, -2, -2], [-2-1j,  0,  0, -1], [-2,  0, -2, -1], [-2, -1, -1,  0]])
+    >>> H = pauli_decompose(A)
+    >>> print(H)
+    (-1.0) [I0 I1]
+    + (-1.5) [X1]
+    + (-0.5) [Y1]
+    + (-1.0) [Z1]
+    + (-1.5) [X0]
+    + (-1.0) [X0 X1]
+    + (-0.5) [X0 Z1]
+    + (1.0) [Y0 Y1]
+    + (-0.5) [Z0 X1]
+    + (-0.5) [Z0 Y1]
+
+    We can return a ``~.PauliSentence`` instance by using the keyword argument ``pauli=True``:
+
+    >>> ps = pauli_decompose(A, pauli=True)
+    >>> print(ps)
+    -1.0 * I
+    + -1.5 * X(1)
+    + -0.5 * Y(1)
+    + -1.0 * Z(1)
+    + -1.5 * X(0)
+    + -1.0 * X(0) @ X(1)
+    + -0.5 * X(0) @ Z(1)
+    + 1.0 * Y(0) @ Y(1)
+    + -0.5 * Z(0) @ X(1)
+    + -0.5 * Z(0) @ Y(1)
+
+    We can also set custom wires using the ``wire_order`` argument:
+
+    >>> ps = pauli_decompose(A, pauli=True, wire_order=['a', 'b'])
+    >>> print(ps)
+    -1.0 * I
+    + -1.5 * X(b)
+    + -0.5 * Y(b)
+    + -1.0 * Z(b)
+    + -1.5 * X(a)
+    + -1.0 * X(a) @ X(b)
+    + -0.5 * X(a) @ Z(b)
+    + 1.0 * Y(a) @ Y(b)
+    + -0.5 * Z(a) @ X(b)
+    + -0.5 * Z(a) @ Y(b)
+    """
+    n = int(np.log2(len(H)))
+    N = 2**n
+
+    if wire_order is not None and len(wire_order) != n:
+        raise ValueError(
+            f"number of wires {len(wire_order)} is not compatible with number of qubits {n}"
+        )
+
+    if wire_order is None:
+        wire_order = range(n)
+
+    if H.shape != (N, N):
+        raise ValueError("The matrix should have shape (2**n, 2**n), for any qubit number n>=1")
+
+    if not np.allclose(H, H.conj().T):
+        raise ValueError("The matrix is not Hermitian")
+
+    obs_lst = []
+    coeffs = []
+
+    for term in product([I, X, Y, Z], repeat=n):
+        matrices = [mat_map[i] for i in term]
+        coeff = np.trace(reduce(np.kron, matrices) @ H) / N
+        coeff = np.real_if_close(coeff).item()
+
+        if not np.allclose(coeff, 0):
+            obs_term = (
+                [(o, w) for w, o in zip(wire_order, term) if o != I]
+                if hide_identity and not all(t == I for t in term)
+                else [(o, w) for w, o in zip(wire_order, term)]
+            )
+
+            if obs_term:
+                coeffs.append(coeff)
+                obs_lst.append(obs_term)
+
+    if pauli:
+        return PauliSentence(
+            {
+                PauliWord({w: o for o, w in obs_n_wires}): coeff
+                for coeff, obs_n_wires in zip(coeffs, obs_lst)
+            }
+        )
+
+    obs = [reduce(matmul, [op_map[o](w) for o, w in obs_term]) for obs_term in obs_lst]
+    return Hamiltonian(coeffs, obs)
 
 
 @singledispatch
diff --git a/pennylane/pauli/pauli_arithmetic.py b/pennylane/pauli/pauli_arithmetic.py
index e7fc6a512..5e73e74bd 100644
--- a/pennylane/pauli/pauli_arithmetic.py
+++ b/pennylane/pauli/pauli_arithmetic.py
@@ -13,6 +13,7 @@
 # limitations under the License.
 """The Pauli arithmetic abstract reduced representation classes"""
 from copy import copy
+from typing import Iterable
 from functools import reduce
 
 import numpy as np
@@ -295,17 +296,16 @@ class PauliSentence(dict):
             ValueError: Can't get the matrix of an empty PauliSentence.
         """
 
-        def _pw_wires(w):
-            """To account for empty pauli_words which represent identity operations."""
+        def _pw_wires(w: Iterable) -> wires.Wires:
+            """Return the native Wires instance for a list of wire labels.
+            w represents the wires of the PauliWord being processed. In case
+            the PauliWord is empty ({}), choose any arbitrary wire from the
+            PauliSentence it is composed in.
+            """
             if w:
                 return wires.Wires(w)
 
-            ps_wires = self.wires
-            if len(ps_wires) > 0:
-                return wires.Wires(
-                    list(ps_wires)[0]
-                )  # return any wire from the Pauli sentence's wires
-            return wires.Wires([])
+            return wires.Wires(list(self.wires)[0]) if len(self.wires) > 0 else wires.Wires([])
 
         if len(self) == 0:
             if wire_order is None or wire_order == wires.Wires([]):
@@ -316,8 +316,8 @@ class PauliSentence(dict):
 
         mats_and_wires_gen = (
             (
-                coeff * pw.to_mat(wire_order=_pw_wires(list(pw.wires)), format=format),
-                _pw_wires(list(pw.wires)),
+                coeff * pw.to_mat(wire_order=_pw_wires(pw.wires), format=format),
+                _pw_wires(pw.wires),
             )
             for pw, coeff in self.items()
         )
diff --git a/pennylane/qchem/tapering.py b/pennylane/qchem/tapering.py
index c2b4c1343..1331d1290 100644
--- a/pennylane/qchem/tapering.py
+++ b/pennylane/qchem/tapering.py
@@ -574,9 +574,8 @@ def _build_generator(operation, wire_order, op_gen=None):
     if op_gen is None:
         if operation.num_params < 1:  # Non-parameterized gates
             gen_mat = 1j * scipy.linalg.logm(qml.matrix(operation, wire_order=wire_order))
-            op_gen = qml.Hamiltonian(
-                *qml.utils.decompose_hamiltonian(gen_mat, wire_order=wire_order, hide_identity=True)
-            )
+            op_gen = qml.pauli_decompose(gen_mat, wire_order=wire_order, hide_identity=True)
+
             qml.simplify(op_gen)
             if op_gen.ops[0].label() == qml.Identity(wires=[wire_order[0]]).label():
                 op_gen -= qml.Hamiltonian([op_gen.coeffs[0]], [qml.Identity(wires=wire_order[0])])
diff --git a/pennylane/utils.py b/pennylane/utils.py
index dd6bab03e..5268b8796 100644
--- a/pennylane/utils.py
+++ b/pennylane/utils.py
@@ -19,9 +19,7 @@ across the PennyLane submodules.
 from collections.abc import Iterable
 import functools
 import inspect
-import itertools
 import numbers
-from operator import matmul
 
 import numpy as np
 import scipy
@@ -29,85 +27,6 @@ import scipy
 import pennylane as qml
 
 
-def decompose_hamiltonian(H, hide_identity=False, wire_order=None):
-    r"""Decomposes a Hermitian matrix into a linear combination of Pauli operators.
-
-    Args:
-        H (array[complex]): a Hermitian matrix of dimension :math:`2^n\times 2^n`
-        hide_identity (bool): does not include the :class:`~.Identity` observable within
-            the tensor products of the decomposition if ``True``
-
-    Returns:
-        tuple[list[float], list[~.Observable]]: a list of coefficients and a list
-        of corresponding tensor products of Pauli observables that decompose the Hamiltonian.
-
-    **Example:**
-
-    We can use this function to compute the Pauli operator decomposition of an arbitrary Hermitian
-    matrix:
-
-    >>> A = np.array(
-    ... [[-2, -2+1j, -2, -2], [-2-1j,  0,  0, -1], [-2,  0, -2, -1], [-2, -1, -1,  0]])
-    >>> coeffs, obs_list = decompose_hamiltonian(A)
-    >>> coeffs
-    [-1.0, -1.5, -0.5, -1.0, -1.5, -1.0, -0.5, 1.0, -0.5, -0.5]
-
-    We can use the output coefficients and tensor Pauli terms to construct a :class:`~.Hamiltonian`:
-
-    >>> H = qml.Hamiltonian(coeffs, obs_list)
-    >>> print(H)
-    (-1.0) [I0 I1]
-    + (-1.5) [X1]
-    + (-0.5) [Y1]
-    + (-1.0) [Z1]
-    + (-1.5) [X0]
-    + (-1.0) [X0 X1]
-    + (-0.5) [X0 Z1]
-    + (1.0) [Y0 Y1]
-    + (-0.5) [Z0 X1]
-    + (-0.5) [Z0 Y1]
-
-    This Hamiltonian can then be used in defining VQE problems using :class:`~.ExpvalCost`.
-    """
-    n = int(np.log2(len(H)))
-    N = 2**n
-
-    if wire_order is None:
-        wire_order = range(n)
-
-    if H.shape != (N, N):
-        raise ValueError(
-            "The Hamiltonian should have shape (2**n, 2**n), for any qubit number n>=1"
-        )
-
-    if not np.allclose(H, H.conj().T):
-        raise ValueError("The Hamiltonian is not Hermitian")
-
-    paulis = [qml.Identity, qml.PauliX, qml.PauliY, qml.PauliZ]
-    obs = []
-    coeffs = []
-
-    for term in itertools.product(paulis, repeat=n):
-        matrices = [i.compute_matrix() for i in term]
-        coeff = np.trace(functools.reduce(np.kron, matrices) @ H) / N
-        coeff = np.real_if_close(coeff).item()
-
-        if not np.allclose(coeff, 0):
-            coeffs.append(coeff)
-
-            if not all(t is qml.Identity for t in term) and hide_identity:
-                obs.append(
-                    functools.reduce(
-                        matmul,
-                        [t(i) for i, t in zip(wire_order, term) if t is not qml.Identity],
-                    )
-                )
-            else:
-                obs.append(functools.reduce(matmul, [t(i) for i, t in enumerate(term)]))
-
-    return coeffs, obs
-
-
 def sparse_hamiltonian(H, wires=None):
     r"""Computes the sparse matrix representation a Hamiltonian in the computational basis.
 


Test patch:
diff --git a/tests/ops/functions/test_generator.py b/tests/ops/functions/test_generator.py
index 9cc5dbcf3..b01f596ce 100644
--- a/tests/ops/functions/test_generator.py
+++ b/tests/ops/functions/test_generator.py
@@ -322,7 +322,7 @@ class TestHamiltonianReturn:
         gen = qml.generator(HermitianOp, format="hamiltonian")(0.5, wires=0)
         assert gen.name == "Hamiltonian"
 
-        expected = qml.Hamiltonian(*qml.utils.decompose_hamiltonian(HermitianOp.H))
+        expected = qml.pauli_decompose(HermitianOp.H, hide_identity=True)
         assert gen.compare(expected)
 
     def test_sparse_hamiltonian(self):
@@ -331,12 +331,12 @@ class TestHamiltonianReturn:
         gen = qml.generator(SparseOp, format="hamiltonian")(0.5, wires=0)
         assert gen.name == "Hamiltonian"
 
-        expected = qml.Hamiltonian(*qml.utils.decompose_hamiltonian(SparseOp.H.toarray()))
+        expected = qml.pauli_decompose(SparseOp.H.toarray(), hide_identity=True)
         assert gen.compare(expected)
 
     def test_hermitian_inverse(self):
         """Test a Hermitian inverted generator is correct"""
-        expected = qml.Hamiltonian(*qml.utils.decompose_hamiltonian(HermitianOp.H))
+        expected = qml.pauli_decompose(HermitianOp.H, hide_identity=True)
 
         gen = qml.generator(HermitianOp(0.5, wires=0).inv(), format="hamiltonian")
         assert gen.name == "Hamiltonian"
@@ -344,7 +344,7 @@ class TestHamiltonianReturn:
 
     def test_sparse_hamiltonian_inverse(self):
         """Test a SparseHamiltonian inverted generator is correct"""
-        expected = qml.Hamiltonian(*qml.utils.decompose_hamiltonian(SparseOp.H.toarray()))
+        expected = qml.pauli_decompose(SparseOp.H.toarray(), hide_identity=True)
 
         gen = qml.generator(SparseOp(0.5, wires=0).inv(), format="hamiltonian")
         assert gen.name == "Hamiltonian"
diff --git a/tests/pauli/test_conversion.py b/tests/pauli/test_conversion.py
index e46ae068a..719981507 100644
--- a/tests/pauli/test_conversion.py
+++ b/tests/pauli/test_conversion.py
@@ -17,6 +17,8 @@ import pytest
 
 import numpy as np
 import pennylane as qml
+from pennylane.operation import Tensor
+from pennylane.ops import Identity, PauliX, PauliY, PauliZ
 from pennylane.pauli import pauli_sentence, PauliWord, PauliSentence
 
 
@@ -27,6 +29,115 @@ test_hamiltonians = [
 ]
 
 
+class TestDecomposition:
+    """Tests the pauli_decompose function"""
+
+    @pytest.mark.parametrize("hamiltonian", [np.ones((3, 3)), np.ones((4, 2)), np.ones((2, 4))])
+    def test_wrong_shape(self, hamiltonian):
+        """Tests that an exception is raised if the Hamiltonian does not have
+        the correct shape"""
+        with pytest.raises(
+            ValueError,
+            match="The matrix should have shape",
+        ):
+            qml.pauli_decompose(hamiltonian)
+
+    def test_not_hermitian(self):
+        """Tests that an exception is raised if the Hamiltonian is not Hermitian, i.e.
+        equal to its own conjugate transpose"""
+        with pytest.raises(ValueError, match="The matrix is not Hermitian"):
+            qml.pauli_decompose(np.array([[1, 2], [3, 4]]))
+
+    def test_hide_identity_true(self):
+        """Tests that there are no Identity observables in the tensor products
+        when hide_identity=True"""
+        H = np.array(np.diag([0, 0, 0, 1]))
+        coeff, obs_list = qml.pauli_decompose(H, hide_identity=True).terms()
+        tensors = filter(lambda obs: isinstance(obs, Tensor), obs_list)
+
+        for tensor in tensors:
+            all_identities = all(isinstance(o, Identity) for o in tensor.obs)
+            no_identities = not any(isinstance(o, Identity) for o in tensor.obs)
+            assert all_identities or no_identities
+
+    def test_hide_identity_true_all_identities(self):
+        """Tests that the all identity operator remains even with hide_identity = True."""
+        H = np.eye(4)
+        coeff, obs_list = qml.pauli_decompose(H, hide_identity=True).terms()
+        tensors = filter(lambda obs: isinstance(obs, Tensor), obs_list)
+
+        for tensor in tensors:
+            assert all(isinstance(o, Identity) for o in tensor.obs)
+
+    @pytest.mark.parametrize("hide_identity", [True, False])
+    @pytest.mark.parametrize("hamiltonian", test_hamiltonians)
+    def test_observable_types(self, hamiltonian, hide_identity):
+        """Tests that the Hamiltonian decomposes into a linear combination of tensors,
+        the identity matrix, and Pauli matrices."""
+        allowed_obs = (Tensor, Identity, PauliX, PauliY, PauliZ)
+
+        decomposed_coeff, decomposed_obs = qml.pauli_decompose(hamiltonian, hide_identity).terms()
+        assert all([isinstance(o, allowed_obs) for o in decomposed_obs])
+
+    @pytest.mark.parametrize("hamiltonian", test_hamiltonians)
+    def test_result_length(self, hamiltonian):
+        """Tests that tensors are composed of a number of terms equal to the number
+        of qubits."""
+        decomposed_coeff, decomposed_obs = qml.pauli_decompose(hamiltonian).terms()
+        n = int(np.log2(len(hamiltonian)))
+
+        tensors = filter(lambda obs: isinstance(obs, Tensor), decomposed_obs)
+        assert all(len(tensor.obs) == n for tensor in tensors)
+
+    @pytest.mark.parametrize("hamiltonian", test_hamiltonians)
+    def test_decomposition(self, hamiltonian):
+        """Tests that pauli_decompose successfully decomposes Hamiltonians into a
+        linear combination of Pauli matrices"""
+        decomposed_coeff, decomposed_obs = qml.pauli_decompose(hamiltonian).terms()
+
+        linear_comb = sum([decomposed_coeff[i] * o.matrix() for i, o in enumerate(decomposed_obs)])
+        assert np.allclose(hamiltonian, linear_comb)
+
+    @pytest.mark.parametrize("hamiltonian", test_hamiltonians)
+    def test_to_paulisentence(self, hamiltonian):
+        """Test that a PauliSentence is returned if the kwarg paulis is set to True"""
+        ps = qml.pauli_decompose(hamiltonian, pauli=True)
+        num_qubits = int(np.log2(len(hamiltonian)))
+
+        assert isinstance(ps, qml.pauli.PauliSentence)
+        assert np.allclose(hamiltonian, ps.to_mat(range(num_qubits)))
+
+    def test_wire_order(self):
+        wire_order1 = ["a", 0]
+        wire_order2 = ["auxiliary", "working"]
+        hamiltonian = np.array(
+            [[-2, -2 + 1j, -2, -2], [-2 - 1j, 0, 0, -1], [-2, 0, -2, -1], [-2, -1, -1, 0]]
+        )
+
+        for wire_order in (wire_order1, wire_order2):
+            h = qml.pauli_decompose(hamiltonian, wire_order=wire_order)
+            ps = qml.pauli_decompose(hamiltonian, pauli=True, wire_order=wire_order)
+
+            assert ps.wires == set(wire_order)
+            assert h.wires.toset() == set(wire_order)
+
+    def test_wire_error(self):
+        wire_order = [0]
+        hamiltonian = np.array(
+            [[-2, -2 + 1j, -2, -2], [-2 - 1j, 0, 0, -1], [-2, 0, -2, -1], [-2, -1, -1, 0]]
+        )
+
+        with pytest.raises(
+            ValueError, match="number of wires 1 is not compatible with number of qubits 2"
+        ):
+            h = qml.pauli_decompose(hamiltonian, wire_order=wire_order)
+
+        with pytest.raises(
+            ValueError, match="number of wires 1 is not compatible with number of qubits 2"
+        ):
+            ps = qml.pauli_decompose(hamiltonian, pauli=True, wire_order=wire_order)
+
+
 class TestPauliSentence:
     """Test the pauli_sentence function."""
 
diff --git a/tests/test_utils.py b/tests/test_utils.py
index bfdcafddf..9567bb506 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -25,9 +25,6 @@ import pennylane as qml
 import pennylane.utils as pu
 import scipy.sparse
 
-from pennylane import Identity, PauliX, PauliY, PauliZ
-from pennylane.operation import Tensor
-
 
 flat_dummy_array = np.linspace(-1, 1, 64)
 test_shapes = [
@@ -67,74 +64,6 @@ U_toffoli = np.diag([1 for i in range(8)])
 U_toffoli[6:8, 6:8] = np.array([[0, 1], [1, 0]])
 
 
-test_hamiltonians = [
-    np.array([[2.5, -0.5], [-0.5, 2.5]]),
-    np.array(np.diag([0, 0, 0, 1])),
-    np.array([[-2, -2 + 1j, -2, -2], [-2 - 1j, 0, 0, -1], [-2, 0, -2, -1], [-2, -1, -1, 0]]),
-]
-
-
-class TestDecomposition:
-    """Tests the decompose_hamiltonian function"""
-
-    @pytest.mark.parametrize("hamiltonian", [np.ones((3, 3)), np.ones((4, 2)), np.ones((2, 4))])
-    def test_wrong_shape(self, hamiltonian):
-        """Tests that an exception is raised if the Hamiltonian does not have
-        the correct shape"""
-        with pytest.raises(
-            ValueError,
-            match="The Hamiltonian should have shape",
-        ):
-            pu.decompose_hamiltonian(hamiltonian)
-
-    def test_not_hermitian(self):
-        """Tests that an exception is raised if the Hamiltonian is not Hermitian, i.e.
-        equal to its own conjugate transpose"""
-        with pytest.raises(ValueError, match="The Hamiltonian is not Hermitian"):
-            pu.decompose_hamiltonian(np.array([[1, 2], [3, 4]]))
-
-    def test_hide_identity_true(self):
-        """Tests that there are no Identity observables in the tensor products
-        when hide_identity=True"""
-        H = np.array(np.diag([0, 0, 0, 1]))
-        coeff, obs_list = pu.decompose_hamiltonian(H, hide_identity=True)
-        tensors = filter(lambda obs: isinstance(obs, Tensor), obs_list)
-
-        for tensor in tensors:
-            all_identities = all(isinstance(o, Identity) for o in tensor.obs)
-            no_identities = not any(isinstance(o, Identity) for o in tensor.obs)
-            assert all_identities or no_identities
-
-    @pytest.mark.parametrize("hide_identity", [True, False])
-    @pytest.mark.parametrize("hamiltonian", test_hamiltonians)
-    def test_observable_types(self, hamiltonian, hide_identity):
-        """Tests that the Hamiltonian decomposes into a linear combination of tensors,
-        the identity matrix, and Pauli matrices."""
-        allowed_obs = (Tensor, Identity, PauliX, PauliY, PauliZ)
-
-        decomposed_coeff, decomposed_obs = pu.decompose_hamiltonian(hamiltonian, hide_identity)
-        assert all([isinstance(o, allowed_obs) for o in decomposed_obs])
-
-    @pytest.mark.parametrize("hamiltonian", test_hamiltonians)
-    def test_result_length(self, hamiltonian):
-        """Tests that tensors are composed of a number of terms equal to the number
-        of qubits."""
-        decomposed_coeff, decomposed_obs = pu.decompose_hamiltonian(hamiltonian)
-        n = int(np.log2(len(hamiltonian)))
-
-        tensors = filter(lambda obs: isinstance(obs, Tensor), decomposed_obs)
-        assert all(len(tensor.obs) == n for tensor in tensors)
-
-    @pytest.mark.parametrize("hamiltonian", test_hamiltonians)
-    def test_decomposition(self, hamiltonian):
-        """Tests that decompose_hamiltonian successfully decomposes Hamiltonians into a
-        linear combination of Pauli matrices"""
-        decomposed_coeff, decomposed_obs = pu.decompose_hamiltonian(hamiltonian)
-
-        linear_comb = sum([decomposed_coeff[i] * o.matrix() for i, o in enumerate(decomposed_obs)])
-        assert np.allclose(hamiltonian, linear_comb)
-
-
 class TestSparse:
     """Tests the sparse_hamiltonian function"""