Repository: PennyLaneAI/pennylane
Pull request #3585 (2022-12-27 13:50:22)

Move code from qml.utils.sparse_hamiltonian to Hamiltonian.sparse_matrix #3236
…trix, added a deprecation warning to qml.utils.sparse_hamiltonian

### Before submitting

Please complete the following checklist when submitting a PR:

- [x] All new features must include a unit test.
      If you've fixed a bug or added code that should be tested, add a test to the
      test directory!

- [x] All new functions and code must be clearly commented and documented.
      If you do make documentation changes, make sure that the docs build and
      render correctly by running `make docs`.

- [x] Ensure that the test suite passes, by running `make test`.

- [x] Add a new entry to the `doc/releases/changelog-dev.md` file, summarizing the
      change, and including a link back to the PR.

- [x] The PennyLane source code conforms to
      [PEP8 standards](https://www.python.org/dev/peps/pep-0008/).
      We check all of our code against [Pylint](https://www.pylint.org/).
      To lint modified files, simply `pip install pylint`, and then
      run `pylint pennylane/path/to/file.py`.

When all the above are checked, delete everything above the dashed
line and fill in the pull request template.

------------------------------------------------------------------------------------------------------------

**Context:**

**Description of the Change:**

**Benefits:**

**Possible Drawbacks:**

**Related GitHub Issues:**


codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3585?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#3585](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3585?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (61ff330) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/880517228462ba01af35d94204e55fc9c47a7cac?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (8805172) will **decrease** coverage by `0.10%`.
> The diff coverage is `18.91%`.

```diff
@@            Coverage Diff             @@
##           master    #3585      +/-   ##
==========================================
- Coverage   99.87%   99.76%   -0.11%     
==========================================
  Files         314      314              
  Lines       28019    28056      +37     
==========================================
+ Hits        27983    27990       +7     
- Misses         36       66      +30     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3585?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/ops/qubit/hamiltonian.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3585/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9oYW1pbHRvbmlhbi5weQ==) | `88.80% <11.76%> (-11.20%)` | :arrow_down: |
| [pennylane/utils.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3585/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% <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)

12mB7693: @albi3ro Thank you for reviewing what I did and for your comments, that helped me a lot.
I think I get what you mean and I will push an update soon. 
Shall I actually also update the tests, where the old function qml.utils.sparse_hamiltonian(..) is called? 
albi3ro: > @albi3ro Thank you for reviewing what I did and for your comments, that helped me a lot.
> I think I get what you mean and I will push an update soon.
> Shall I actually also update the tests, where the old function qml.utils.sparse_hamiltonian(..) is called?

Updating the tests where `qml.utils.spase_hamiltonian` is called (outside of unit tests) would be great! But I won't necessarily block approval based on that being accomplished.
12mB7693: 

> > @albi3ro Thank you for reviewing what I did and for your comments, that helped me a lot.
> > I think I get what you mean and I will push an update soon.
> > Shall I actually also update the tests, where the old function qml.utils.sparse_hamiltonian(..) is called?
> 
> Updating the tests where `qml.utils.spase_hamiltonian` is called (outside of unit tests) would be great! But I won't necessarily block approval based on that being accomplished.

I'm not sure if I understood you correctly regarding the updates of the tests. Do you mean with "outside of unit tests" that I should not update the unit tests? 
12mB7693: Not sure why the documentation check fails, I ran "make doc" locally and I get the warning
`releases/changelog-dev.md:391: WARNING: Definition list ends without a blank line; unexpected unindent.
`
However, I do not see how to fix it.
12mB7693: No problem, I'm glad that I could help. This task was definitely a good starting point as a new contributor. 
Thank you for your support @albi3ro!

Patch:
diff --git a/doc/development/deprecations.rst b/doc/development/deprecations.rst
index c7b2b0f1f..1d6dbb83e 100644
--- a/doc/development/deprecations.rst
+++ b/doc/development/deprecations.rst
@@ -71,6 +71,11 @@ Pending deprecations
 
   - Deprecated in v0.29.
   - Will be removed in v0.30.
+  
+* The ``qml.utils.sparse_hamiltonian`` function is deprecated. ``~.Hamiltonian.sparse_matrix`` should be used instead.
+
+  - Deprecated in v0.29
+  - Removed in v0.31
 
 Completed deprecation cycles
 ----------------------------
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index ecc5a7940..9e8f3e39f 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -488,6 +488,9 @@
   `qml.gradients.SUPPORTED_GRADIENT_KWARGS`.
   [(#3526)](https://github.com/PennyLaneAI/pennylane/pull/3526)
 
+* Moved `qml.utils.sparse_hamiltonian` function to `~.Hamiltonian.sparse_matrix` method.
+  [(#3585)](https://github.com/PennyLaneAI/pennylane/pull/3585)
+
 * The `PauliSentence.operation()` method has been improved to avoid instantiating an `SProd` operator when
   the coefficient is equal to 1.
   [(#3595)](https://github.com/PennyLaneAI/pennylane/pull/3595)
@@ -636,6 +639,10 @@
 
 <h3>Deprecations</h3>
 
+* `qml.utils.sparse_hamiltonian` function has been deprecated, and usage will now raise a warning.
+  Instead, one should use the `~.Hamiltonian.sparse_matrix` method.
+  [(#3585)](https://github.com/PennyLaneAI/pennylane/pull/3585)
+
 * Deprecate the `collections` module.
   [(#3686)](https://github.com/PennyLaneAI/pennylane/pull/3686)
 
@@ -719,6 +726,7 @@ Guillermo Alonso-Linaje
 Juan Miguel Arrazola
 Ikko Ashimine
 Utkarsh Azad
+Miriam Beddig
 Cristian Boghiu
 Astral Cai
 Isaac De Vlugt
diff --git a/pennylane/devices/default_qubit.py b/pennylane/devices/default_qubit.py
index b58c2cee7..3760f8a34 100644
--- a/pennylane/devices/default_qubit.py
+++ b/pennylane/devices/default_qubit.py
@@ -600,10 +600,8 @@ class DefaultQubit(QubitDevice):
             else:
                 # Coefficients and the state are not trainable, we can be more
                 # efficient in how we compute the Hamiltonian sparse matrix.
-                if observable.name == "Hamiltonian":
-                    Hmat = qml.utils.sparse_hamiltonian(observable, wires=self.wires)
-                elif observable.name == "SparseHamiltonian":
-                    Hmat = observable.sparse_matrix()
+                if observable.name in {"Hamiltonian", "SparseHamiltonian"}:
+                    Hmat = observable.sparse_matrix(wire_order=self.wires)
 
                 state = qml.math.toarray(self.state)
                 if self._ndim(state) == 2:
diff --git a/pennylane/ops/functions/matrix.py b/pennylane/ops/functions/matrix.py
index 29f3b6ca9..5bc5627a2 100644
--- a/pennylane/ops/functions/matrix.py
+++ b/pennylane/ops/functions/matrix.py
@@ -122,7 +122,7 @@ def matrix(op, *, wire_order=None):
         op = 1.0 * op  # convert to a Hamiltonian
 
     if isinstance(op, qml.Hamiltonian):
-        return qml.utils.sparse_hamiltonian(op, wires=wire_order).toarray()
+        return op.sparse_matrix(wire_order=wire_order).toarray()
 
     return op.matrix(wire_order=wire_order)
 
diff --git a/pennylane/ops/qubit/hamiltonian.py b/pennylane/ops/qubit/hamiltonian.py
index c2081a7be..06a1d6de0 100644
--- a/pennylane/ops/qubit/hamiltonian.py
+++ b/pennylane/ops/qubit/hamiltonian.py
@@ -20,7 +20,10 @@ import itertools
 import numbers
 from collections.abc import Iterable
 from copy import copy
+import functools
 from typing import List
+import scipy
+
 
 import pennylane as qml
 from pennylane import numpy as np
@@ -340,6 +343,92 @@ class Hamiltonian(Observable):
                 self.ops, grouping_type=grouping_type, method=method
             )
 
+    def sparse_matrix(self, wire_order=None):
+        r"""Computes the sparse matrix representation of a Hamiltonian in the computational basis.
+
+        Args:
+            wire_order (Iterable): global wire order, must contain all wire labels from the operator's wires.
+                If not provided, the default order of the wires (self.wires) of the Hamiltonian is used.
+
+        Returns:
+            csr_matrix: a sparse matrix in scipy Compressed Sparse Row (CSR) format with dimension
+            :math:`(2^n, 2^n)`, where :math:`n` is the number of wires
+
+        **Example:**
+
+        >>> coeffs = [1, -0.45]
+        >>> obs = [qml.PauliZ(0) @ qml.PauliZ(1), qml.PauliY(0) @ qml.PauliZ(1)]
+        >>> H = qml.Hamiltonian(coeffs, obs)
+        >>> H_sparse = H.sparse_matrix()
+        >>> H_sparse
+        <4x4 sparse matrix of type '<class 'numpy.complex128'>'
+                with 8 stored elements in Compressed Sparse Row format>
+
+        The resulting sparse matrix can be either used directly or transformed into a numpy array:
+
+        >>> H_sparse.toarray()
+        array([[ 1.+0.j  ,  0.+0.j  ,  0.+0.45j,  0.+0.j  ],
+               [ 0.+0.j  , -1.+0.j  ,  0.+0.j  ,  0.-0.45j],
+               [ 0.-0.45j,  0.+0.j  , -1.+0.j  ,  0.+0.j  ],
+               [ 0.+0.j  ,  0.+0.45j,  0.+0.j  ,  1.+0.j  ]])
+        """
+        if wire_order is None:
+            wires = self.wires
+        else:
+            wires = wire_order
+        n = len(wires)
+        matrix = scipy.sparse.csr_matrix((2**n, 2**n), dtype="complex128")
+
+        coeffs = qml.math.toarray(self.data)
+
+        temp_mats = []
+        for coeff, op in zip(coeffs, self.ops):
+            obs = []
+            for o in qml.operation.Tensor(op).obs:
+                if len(o.wires) > 1:
+                    # todo: deal with operations created from multi-qubit operations such as Hermitian
+                    raise ValueError(
+                        f"Can only sparsify Hamiltonians whose constituent observables consist of "
+                        f"(tensor products of) single-qubit operators; got {op}."
+                    )
+                obs.append(o.matrix())
+
+            # Array to store the single-wire observables which will be Kronecker producted together
+            mat = []
+            # i_count tracks the number of consecutive single-wire identity matrices encountered
+            # in order to avoid unnecessary Kronecker products, since I_n x I_m = I_{n+m}
+            i_count = 0
+            for wire_lab in wires:
+                if wire_lab in op.wires:
+                    if i_count > 0:
+                        mat.append(scipy.sparse.eye(2**i_count, format="coo"))
+                    i_count = 0
+                    idx = op.wires.index(wire_lab)
+                    # obs is an array storing the single-wire observables which
+                    # make up the full Hamiltonian term
+                    sp_obs = scipy.sparse.coo_matrix(obs[idx])
+                    mat.append(sp_obs)
+                else:
+                    i_count += 1
+
+            if i_count > 0:
+                mat.append(scipy.sparse.eye(2**i_count, format="coo"))
+
+            red_mat = (
+                functools.reduce(lambda i, j: scipy.sparse.kron(i, j, format="coo"), mat) * coeff
+            )
+
+            temp_mats.append(red_mat.tocsr())
+            # Value of 100 arrived at empirically to balance time savings vs memory use. At this point
+            # the `temp_mats` are summed into the final result and the temporary storage array is
+            # cleared.
+            if (len(temp_mats) % 100) == 0:
+                matrix += sum(temp_mats)
+                temp_mats = []
+
+        matrix += sum(temp_mats)
+        return matrix
+
     def simplify(self):
         r"""Simplifies the Hamiltonian by combining like-terms.
 
diff --git a/pennylane/optimize/lie_algebra.py b/pennylane/optimize/lie_algebra.py
index 252d66c13..9329ccbc2 100644
--- a/pennylane/optimize/lie_algebra.py
+++ b/pennylane/optimize/lie_algebra.py
@@ -59,7 +59,7 @@ def append_time_evolution(tape, riemannian_gradient, t, n, exact=False):
         qml.apply(obj)
     if exact:
         qml.QubitUnitary(
-            expm(-1j * t * qml.utils.sparse_hamiltonian(riemannian_gradient).toarray()),
+            expm(-1j * t * riemannian_gradient.sparse_matrix().toarray()),
             wires=range(len(riemannian_gradient.wires)),
         )
     else:
diff --git a/pennylane/utils.py b/pennylane/utils.py
index 5268b8796..4ec976899 100644
--- a/pennylane/utils.py
+++ b/pennylane/utils.py
@@ -20,21 +20,24 @@ from collections.abc import Iterable
 import functools
 import inspect
 import numbers
+import warnings
+
 
 import numpy as np
-import scipy
 
 import pennylane as qml
 
 
 def sparse_hamiltonian(H, wires=None):
-    r"""Computes the sparse matrix representation a Hamiltonian in the computational basis.
+    r"""Warning: This method is deprecated. Use :meth:~.Hamiltonian.sparse_matrix` instead.
+
+    Computes the sparse matrix representation a Hamiltonian in the computational basis.
 
     Args:
         H (~.Hamiltonian): Hamiltonian operator for which the matrix representation should be
-         computed
+            computed
         wires (Iterable): Wire labels that indicate the order of wires according to which the matrix
-         is constructed. If not profided, ``H.wires`` is used.
+            is constructed. If not profided, ``H.wires`` is used.
 
     Returns:
         csr_matrix: a sparse matrix in scipy Compressed Sparse Row (CSR) format with dimension
@@ -59,65 +62,22 @@ def sparse_hamiltonian(H, wires=None):
            [ 0.+0.j  , -1.+0.j  ,  0.+0.j  ,  0.-0.45j],
            [ 0.-0.45j,  0.+0.j  , -1.+0.j  ,  0.+0.j  ],
            [ 0.+0.j  ,  0.+0.45j,  0.+0.j  ,  1.+0.j  ]])
+
     """
+    warnings.warn(
+        "The method sparse_hamiltonian is deprecated. Please use the method "
+        "sparse_matrix of the Hamiltonian operator instead.",
+        UserWarning,
+    )
     if not isinstance(H, qml.Hamiltonian):
         raise TypeError("Passed Hamiltonian must be of type `qml.Hamiltonian`")
 
-    if wires is None:
+    if wires is None:  # not sure if this if-else is still necessary
         wires = H.wires
     else:
         wires = qml.wires.Wires(wires)
 
-    n = len(wires)
-    matrix = scipy.sparse.csr_matrix((2**n, 2**n), dtype="complex128")
-
-    coeffs = qml.math.toarray(H.data)
-
-    temp_mats = []
-    for coeff, op in zip(coeffs, H.ops):
-        obs = []
-        for o in qml.operation.Tensor(op).obs:
-            if len(o.wires) > 1:
-                # todo: deal with operations created from multi-qubit operations such as Hermitian
-                raise ValueError(
-                    f"Can only sparsify Hamiltonians whose constituent observables consist of "
-                    f"(tensor products of) single-qubit operators; got {op}."
-                )
-            obs.append(o.matrix())
-
-        # Array to store the single-wire observables which will be Kronecker producted together
-        mat = []
-        # i_count tracks the number of consecutive single-wire identity matrices encountered
-        # in order to avoid unnecessary Kronecker products, since I_n x I_m = I_{n+m}
-        i_count = 0
-        for wire_lab in wires:
-            if wire_lab in op.wires:
-                if i_count > 0:
-                    mat.append(scipy.sparse.eye(2**i_count, format="coo"))
-                i_count = 0
-                idx = op.wires.index(wire_lab)
-                # obs is an array storing the single-wire observables which
-                # make up the full Hamiltonian term
-                sp_obs = scipy.sparse.coo_matrix(obs[idx])
-                mat.append(sp_obs)
-            else:
-                i_count += 1
-
-        if i_count > 0:
-            mat.append(scipy.sparse.eye(2**i_count, format="coo"))
-
-        red_mat = functools.reduce(lambda i, j: scipy.sparse.kron(i, j, format="coo"), mat) * coeff
-
-        temp_mats.append(red_mat.tocsr())
-        # Value of 100 arrived at empirically to balance time savings vs memory use. At this point
-        # the `temp_mats` are summed into the final result and the temporary storage array is
-        # cleared.
-        if (len(temp_mats) % 100) == 0:
-            matrix += sum(temp_mats)
-            temp_mats = []
-
-    matrix += sum(temp_mats)
-    return matrix
+    return H.sparse_matrix(wire_order=wires)
 
 
 def _flatten(x):


Test patch:
diff --git a/tests/ops/qubit/test_hamiltonian.py b/tests/ops/qubit/test_hamiltonian.py
index 49d6ef74e..ae061e894 100644
--- a/tests/ops/qubit/test_hamiltonian.py
+++ b/tests/ops/qubit/test_hamiltonian.py
@@ -19,6 +19,7 @@ from unittest.mock import patch
 
 import numpy as np
 import pytest
+import scipy
 
 import pennylane as qml
 from pennylane import numpy as pnp
@@ -1217,6 +1218,165 @@ class TestHamiltonianArithmeticAutograd:
         assert H.compare(H1 @ H2)
 
 
+class TestHamiltonianSparseMatrix:
+    """Tests for sparse matrix representation."""
+
+    @pytest.mark.parametrize(
+        ("coeffs", "obs", "wires", "ref_matrix"),
+        [
+            (
+                [1, -0.45],
+                [qml.PauliZ(0) @ qml.PauliZ(1), qml.PauliY(0) @ qml.PauliZ(1)],
+                None,
+                np.array(
+                    [
+                        [1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.45j, 0.0 + 0.0j],
+                        [0.0 + 0.0j, -1.0 + 0.0j, 0.0 + 0.0j, 0.0 - 0.45j],
+                        [0.0 - 0.45j, 0.0 + 0.0j, -1.0 + 0.0j, 0.0 + 0.0j],
+                        [0.0 + 0.0j, 0.0 + 0.45j, 0.0 + 0.0j, 1.0 + 0.0j],
+                    ]
+                ),
+            ),
+            (
+                [0.1],
+                [qml.PauliZ("b") @ qml.PauliX("a")],
+                ["a", "c", "b"],
+                np.array(
+                    [
+                        [
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.1 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                        ],
+                        [
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            -0.1 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                        ],
+                        [
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.1 + 0.0j,
+                            0.0 + 0.0j,
+                        ],
+                        [
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            -0.1 + 0.0j,
+                        ],
+                        [
+                            0.1 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                        ],
+                        [
+                            0.0 + 0.0j,
+                            -0.1 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                        ],
+                        [
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.1 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                        ],
+                        [
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            -0.1 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                            0.0 + 0.0j,
+                        ],
+                    ]
+                ),
+            ),
+            (
+                [0.21, -0.78, 0.52],
+                [
+                    qml.PauliZ(0) @ qml.PauliZ(1),
+                    qml.PauliX(0) @ qml.PauliZ(1),
+                    qml.PauliY(0) @ qml.PauliZ(1),
+                ],
+                None,
+                np.array(
+                    [
+                        [0.21 + 0.0j, 0.0 + 0.0j, -0.78 - 0.52j, 0.0 + 0.0j],
+                        [0.0 + 0.0j, -0.21 + 0.0j, 0.0 + 0.0j, 0.78 + 0.52j],
+                        [-0.78 + 0.52j, 0.0 + 0.0j, -0.21 + 0.0j, 0.0 + 0.0j],
+                        [0.0 + 0.0j, 0.78 - 0.52j, 0.0 + 0.0j, 0.21 + 0.0j],
+                    ]
+                ),
+            ),
+        ],
+    )
+    def test_sparse_matrix(self, coeffs, obs, wires, ref_matrix):
+        """Tests that sparse_hamiltonian returns a correct sparse matrix"""
+        H = qml.Hamiltonian(coeffs, obs)
+
+        sparse_matrix = H.sparse_matrix(wire_order=wires)
+
+        assert np.allclose(sparse_matrix.toarray(), ref_matrix)
+
+    def test_sparse_format(self):
+        """Tests that sparse_hamiltonian returns a scipy.sparse.csr_matrix object"""
+
+        coeffs = [-0.25, 0.75]
+        obs = [
+            qml.PauliX(wires=[0]) @ qml.PauliZ(wires=[1]),
+            qml.PauliY(wires=[0]) @ qml.PauliZ(wires=[1]),
+        ]
+        H = qml.Hamiltonian(coeffs, obs)
+
+        sparse_matrix = H.sparse_matrix()
+
+        assert isinstance(sparse_matrix, scipy.sparse.csr_matrix)
+
+    def test_observable_error(self):
+        """Tests that an error is thrown if the observables are themselves constructed from multi-qubit
+        operations."""
+        with pytest.raises(ValueError, match="Can only sparsify Hamiltonians"):
+            H = qml.Hamiltonian(
+                [0.1], [qml.PauliZ("c") @ qml.Hermitian(np.eye(4), wires=["a", "b"])]
+            )
+            H.sparse_matrix(wire_order=["a", "c", "b"])
+
+
 class TestHamiltonianArithmeticJax:
     """Tests creation of Hamiltonians using arithmetic
     operations with jax tensor coefficients."""
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 9567bb506..abfd82ace 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -199,6 +199,16 @@ class TestSparse:
 
         assert np.allclose(sparse_matrix.toarray(), ref_matrix)
 
+    def test_warning(self):
+        """Tests that a warning is raised if sparse_hamiltonian function is called"""
+
+        coeffs = [0.2, -0.543]
+        obs = [qml.PauliX(0) @ qml.PauliZ(1), qml.PauliZ(0) @ qml.Hadamard(2)]
+        H = qml.Hamiltonian(coeffs, obs)
+
+        with pytest.warns(UserWarning, match="sparse_hamiltonian is deprecated"):
+            qml.utils.sparse_hamiltonian(H)
+
     def test_sparse_format(self):
         """Tests that sparse_hamiltonian returns a scipy.sparse.csr_matrix object"""