Repository: PennyLaneAI/pennylane (python)

Issue:
[BUG] Error when multiplying certain Hamiltonian observables
### Expected behavior

I've been playing around with the different observable types in Pennylane and noticed we do support multiplying Hamiltonians together, both from tensor products of Paulis and from Hermitian observables:

```py
import pennylane as qml

tensor = qml.PauliX(0) @ qml.PauliX(1)
herm = qml.Hermitian([[1,0], [0,1]], wires=4)

qml.Hamiltonian([1.0], [tensor]) @ qml.Hamiltonian([1.0], [herm])
```

### Actual behavior

However, for this case I get an unexpected error:

```py
import pennylane as qml

tensor = qml.PauliX(0) @ qml.PauliX(1)
herm = qml.Hermitian([[1,0], [0,1]], wires=4)

qml.Hamiltonian([1.0, 1.0], [tensor, qml.PauliX(2)]) @ qml.Hamiltonian([1.0], [herm])
```

### Additional information

_No response_

### Source code

_No response_

### Tracebacks

```shell
AttributeError                            Traceback (most recent call last)
Untitled-1.ipynb Cell 5 in <cell line: 6>()
      3 tensor = qml.PauliX(0) @ qml.PauliX(1)
      4 herm = qml.Hermitian([[1,0], [0,1]], wires=4)
----> 6 qml.Hamiltonian([1.0, 1.0], [tensor, qml.PauliX(2)]) @ qml.Hamiltonian([1.0], [herm])

File ~/.local/lib/python3.10/site-packages/pennylane/ops/qubit/hamiltonian.py:543, in Hamiltonian.__matmul__(self, H)
    540     ops_list = itertools.product(ops1, ops2)
    541     terms = [qml.operation.Tensor(t[0], t[1]) for t in ops_list]
--> 543     return qml.Hamiltonian(coeffs, terms, simplify=True)
    545 if isinstance(H, (Tensor, Observable)):
    546     terms = [op @ H for op in ops1]

File ~/.local/lib/python3.10/site-packages/pennylane/ops/qubit/hamiltonian.py:203, in Hamiltonian.__init__(self, coeffs, observables, simplify, grouping_type, method, id, do_queue)
    200 self._grouping_indices = None
    202 if simplify:
--> 203     self.simplify()
    204 if grouping_type is not None:
    205     with qml.tape.stop_recording():

File ~/.local/lib/python3.10/site-packages/pennylane/ops/qubit/hamiltonian.py:374, in Hamiltonian.simplify(self)
    371 c = self.coeffs[i]
    372 op = op if isinstance(op, Tensor) else Tensor(op)
--> 374 ind = next((j for j, o in enumerate(new_ops) if op.compare(o)), None)
...
-> 1589     parameters = tuple(param.tobytes() for param in ob.parameters)
   1590     tensor.add((ob.name, ob.wires, parameters))
   1592 return tensor

AttributeError: 'list' object has no attribute 'tobytes'
```


### System information

```shell
Name: PennyLane
Version: 0.26.0
Summary: PennyLane is a Python quantum machine learning library by Xanadu Inc.
Home-page: https://github.com/XanaduAI/pennylane
Author: 
Author-email: 
License: Apache License 2.0
Location: /home/david/.local/lib/python3.10/site-packages
Requires: appdirs, autograd, autoray, cachetools, networkx, numpy, pennylane-lightning, retworkx, scipy, semantic-version, toml
Required-by: PennyLane-Lightning

Platform info:           Linux-5.15.0-1020-aws-x86_64-with-glibc2.31
Python version:          3.10.7
Numpy version:           1.22.4
Scipy version:           1.8.1
Installed devices:
- lightning.qubit (PennyLane-Lightning-0.26.0)
- default.gaussian (PennyLane-0.26.0)
- default.mixed (PennyLane-0.26.0)
- default.qubit (PennyLane-0.26.0)
- default.qubit.autograd (PennyLane-0.26.0)
- default.qubit.jax (PennyLane-0.26.0)
- default.qubit.tf (PennyLane-0.26.0)
- default.qubit.torch (PennyLane-0.26.0)
- default.qutrit (PennyLane-0.26.0)
```


### Existing GitHub issues

- [X] I have searched existing GitHub issues to make sure the issue does not already exist.

Pull request description:
Cast input of `qml.Hermitian` to array
If not an array, but a sequence is inputted to `qml.Hermitian` (contrary to its docstring), users may see non-intuitive behaviour.

This PR considers the input as an array and adds a validation check when creating such observables.

Closes #3164.

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 21471d126..158ac49c3 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -86,6 +86,9 @@ keyword argument when using `GellMann`, which determines which of the 8 Gell-Man
 
 <h3>Improvements</h3>
 
+* The matrix passed to `qml.Hermitian` is validated when creating the observable if the input is not abstract.
+  [(#3181)](https://github.com/PennyLaneAI/pennylane/pull/3181)
+
 * Added a new `pennylane.tape.QuantumScript` class that contains all the non-queuing behavior of `QuantumTape`. Now `QuantumTape` inherits from `QuantumScript` as well
   as `AnnotatedQueue`.
   This is a developer-facing change, and users should not manipulate `QuantumScript` directly.  Instead, they
@@ -242,6 +245,9 @@ keyword argument when using `GellMann`, which determines which of the 8 Gell-Man
 
 <h3>Bug fixes</h3>
 
+* Users no longer see unintuitive errors when inputing sequences to `qml.Hermitian`.
+  [(#3181)](https://github.com/PennyLaneAI/pennylane/pull/3181)
+
 * `ControlledQubitUnitary.pow` now copies over the `control_values`.
   [(#3206)](https://github.com/PennyLaneAI/pennylane/pull/3206)
 
diff --git a/pennylane/ops/qubit/observables.py b/pennylane/ops/qubit/observables.py
index e57540400..7d1e50ca8 100644
--- a/pennylane/ops/qubit/observables.py
+++ b/pennylane/ops/qubit/observables.py
@@ -18,6 +18,7 @@ excepting the Pauli gates and Hadamard gate in ``non_parametric_ops.py``.
 
 from copy import copy
 
+from collections.abc import Sequence
 import numpy as np
 from scipy.sparse import csr_matrix
 
@@ -49,7 +50,7 @@ class Hermitian(Observable):
     * Gradient recipe: None
 
     Args:
-        A (array): square hermitian matrix
+        A (array or Sequence): square hermitian matrix
         wires (Sequence[int] or int): the wire(s) the operation acts on
         do_queue (bool): Indicates whether the operator should be
             immediately pushed into the Operator queue (optional)
@@ -60,11 +61,43 @@ class Hermitian(Observable):
     """int: Number of trainable parameters that the operator depends on."""
 
     grad_method = "F"
+
+    # Qubit case
+    _num_basis_states = 2
     _eigs = {}
 
     def __init__(self, A, wires, do_queue=True, id=None):
+        A = qml.math.asarray(A)
+        if not qml.math.is_abstract(A):
+            if isinstance(wires, Sequence) and not isinstance(wires, str):
+                if len(wires) == 0:
+                    raise ValueError(
+                        "Hermitian: wrong number of wires. At least one wire has to be given."
+                    )
+                expected_mx_shape = self._num_basis_states ** len(wires)
+            else:
+                # Assumably wires is an int; further validation checks are performed by calling super().__init__
+                expected_mx_shape = self._num_basis_states
+
+            Hermitian._validate_input(A, expected_mx_shape)
+
         super().__init__(A, wires=wires, do_queue=do_queue, id=id)
 
+    @staticmethod
+    def _validate_input(A, expected_mx_shape=None):
+        """Validate the input matrix."""
+        if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
+            raise ValueError("Observable must be a square matrix.")
+
+        if expected_mx_shape is not None and A.shape[0] != expected_mx_shape:
+            raise ValueError(
+                f"Expected input matrix to have shape {expected_mx_shape}x{expected_mx_shape}, but "
+                f"a matrix with shape {A.shape[0]}x{A.shape[0]} was passed."
+            )
+
+        if not qml.math.allclose(A, qml.math.T(qml.math.conj(A))):
+            raise ValueError("Observable must be Hermitian.")
+
     def label(self, decimals=None, base_label=None, cache=None):
         return super().label(decimals=decimals, base_label=base_label or "𝓗", cache=cache)
 
@@ -91,13 +124,7 @@ class Hermitian(Observable):
          [ 1.+2.j -1.+0.j]]
         """
         A = qml.math.asarray(A)
-
-        if A.shape[0] != A.shape[1]:
-            raise ValueError("Observable must be a square matrix.")
-
-        if not qml.math.allclose(A, qml.math.T(qml.math.conj(A))):
-            raise ValueError("Observable must be Hermitian.")
-
+        Hermitian._validate_input(A)
         return A
 
     @property
diff --git a/pennylane/ops/qutrit/observables.py b/pennylane/ops/qutrit/observables.py
index 2e2223a04..45c05a48f 100644
--- a/pennylane/ops/qutrit/observables.py
+++ b/pennylane/ops/qutrit/observables.py
@@ -69,6 +69,8 @@ class THermitian(Hermitian):
 
     """
 
+    # Qutrit case
+    _num_basis_states = 3
     _eigs = {}
 
     # This method is overridden to update the docstring.


Test patch:
diff --git a/tests/ops/qubit/test_hamiltonian.py b/tests/ops/qubit/test_hamiltonian.py
index ec7b88028..f953ae3d6 100644
--- a/tests/ops/qubit/test_hamiltonian.py
+++ b/tests/ops/qubit/test_hamiltonian.py
@@ -957,6 +957,14 @@ class TestHamiltonian:
         for coeff1, coeff2 in zip(mapped_h.coeffs, h.coeffs):
             assert coeff1 == coeff2
 
+    def test_hermitian_tensor_prod(self):
+        """Test that the tensor product of a Hamiltonian with Hermitian observable works."""
+        tensor = qml.PauliX(0) @ qml.PauliX(1)
+        herm = qml.Hermitian([[1, 0], [0, 1]], wires=4)
+
+        ham = qml.Hamiltonian([1.0, 1.0], [tensor, qml.PauliX(2)]) @ qml.Hamiltonian([1.0], [herm])
+        assert isinstance(ham, qml.Hamiltonian)
+
 
 class TestHamiltonianCoefficients:
     """Test the creation of a Hamiltonian"""
diff --git a/tests/ops/qubit/test_observables.py b/tests/ops/qubit/test_observables.py
index 8b52c7108..e41dd7790 100644
--- a/tests/ops/qubit/test_observables.py
+++ b/tests/ops/qubit/test_observables.py
@@ -148,6 +148,36 @@ class TestSimpleObservables:
 class TestHermitian:
     """Test the Hermitian observable"""
 
+    def test_hermitian_creation_exceptions(self):
+        """Tests that the hermitian matrix method raises the proper errors."""
+        H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
+
+        # test non-square matrix
+        with pytest.raises(ValueError, match="must be a square matrix"):
+            qml.Hermitian(H[1:], wires=0)
+
+        H1 = np.array([[1]]) / np.sqrt(2)
+
+        # test matrix with incorrect dimensions
+        with pytest.raises(ValueError, match="Expected input matrix to have shape"):
+            qml.Hermitian(H1, wires=[0])
+
+        # test non-Hermitian matrix
+        H2 = H.copy()
+        H2[0, 1] = 2
+        with pytest.raises(ValueError, match="must be Hermitian"):
+            qml.Hermitian(H2, wires=0)
+
+    def test_ragged_input_raises(self):
+        """Tests that an error is raised if the input to Hermitian is ragged."""
+        H = [[1, 0], [0, 1, 2]]
+
+        with pytest.warns(
+            np.VisibleDeprecationWarning, match="Creating an ndarray from ragged nested sequences"
+        ):
+            with pytest.raises(ValueError, match="must be a square matrix"):
+                qml.Hermitian(H, wires=0)
+
     @pytest.mark.parametrize("observable, eigvals, eigvecs", EIGVALS_TEST_DATA)
     def test_hermitian_eigegendecomposition_single_wire(self, observable, eigvals, eigvecs, tol):
         """Tests that the eigendecomposition property of the Hermitian class returns the correct results
@@ -361,7 +391,7 @@ class TestHermitian:
 
     def test_hermitian_empty_wire_list_error(self):
         """Tests that the hermitian operator raises an error when instantiated with wires=[]."""
-        herm_mat = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
+        herm_mat = np.array([]).reshape((0, 0))
 
         with pytest.raises(ValueError, match="wrong number of wires"):
             qml.Hermitian(herm_mat, wires=[])
diff --git a/tests/test_measurements.py b/tests/test_measurements.py
index c1e735a8c..ac68d7bd6 100644
--- a/tests/test_measurements.py
+++ b/tests/test_measurements.py
@@ -1158,7 +1158,7 @@ class TestProperties:
     def test_wires_match_observable(self):
         """Test that the wires of the measurement process
         match an internal observable"""
-        obs = qml.Hermitian(np.diag([1, 2, 3]), wires=["a", "b", "c"])
+        obs = qml.Hermitian(np.diag([1, 2, 3, 4]), wires=["a", "b"])
         m = MeasurementProcess(Expectation, obs=obs)
 
         assert np.all(m.wires == obs.wires)
@@ -1166,19 +1166,19 @@ class TestProperties:
     def test_eigvals_match_observable(self):
         """Test that the eigenvalues of the measurement process
         match an internal observable"""
-        obs = qml.Hermitian(np.diag([1, 2, 3]), wires=[0, 1, 2])
+        obs = qml.Hermitian(np.diag([1, 2, 3, 4]), wires=[0, 1])
         m = MeasurementProcess(Expectation, obs=obs)
 
-        assert np.all(m.eigvals() == np.array([1, 2, 3]))
+        assert np.all(m.eigvals() == np.array([1, 2, 3, 4]))
 
         # changing the observable data should be reflected
-        obs.data = [np.diag([5, 6, 7])]
-        assert np.all(m.eigvals() == np.array([5, 6, 7]))
+        obs.data = [np.diag([5, 6, 7, 8])]
+        assert np.all(m.eigvals() == np.array([5, 6, 7, 8]))
 
     def test_error_obs_and_eigvals(self):
         """Test that providing both eigenvalues and an observable
         results in an error"""
-        obs = qml.Hermitian(np.diag([1, 2, 3]), wires=[0, 1, 2])
+        obs = qml.Hermitian(np.diag([1, 2, 3, 4]), wires=[0, 1])
 
         with pytest.raises(ValueError, match="Cannot set the eigenvalues"):
             MeasurementProcess(Expectation, obs=obs, eigvals=[0, 1])
@@ -1186,7 +1186,7 @@ class TestProperties:
     def test_error_obs_and_wires(self):
         """Test that providing both wires and an observable
         results in an error"""
-        obs = qml.Hermitian(np.diag([1, 2, 3]), wires=[0, 1, 2])
+        obs = qml.Hermitian(np.diag([1, 2, 3, 4]), wires=[0, 1])
 
         with pytest.raises(ValueError, match="Cannot set the wires"):
             MeasurementProcess(Expectation, obs=obs, wires=qml.wires.Wires([0, 1]))
diff --git a/tests/test_operation.py b/tests/test_operation.py
index d8f876252..17d0247f2 100644
--- a/tests/test_operation.py
+++ b/tests/test_operation.py
@@ -1204,7 +1204,7 @@ class TestTensor:
 
     def test_num_wires(self):
         """Test that the correct number of wires is returned"""
-        p = np.array([0.5])
+        p = np.eye(4)
         X = qml.PauliX(0)
         Y = qml.Hermitian(p, wires=[1, 2])
         t = Tensor(X, Y)
@@ -1212,7 +1212,7 @@ class TestTensor:
 
     def test_wires(self):
         """Test that the correct nested list of wires is returned"""
-        p = np.array([0.5])
+        p = np.eye(4)
         X = qml.PauliX(0)
         Y = qml.Hermitian(p, wires=[1, 2])
         t = Tensor(X, Y)
@@ -1220,7 +1220,7 @@ class TestTensor:
 
     def test_params(self):
         """Test that the correct flattened list of parameters is returned"""
-        p = np.array([0.5])
+        p = np.eye(4)
         X = qml.PauliX(0)
         Y = qml.Hermitian(p, wires=[1, 2])
         t = Tensor(X, Y)
@@ -1228,7 +1228,7 @@ class TestTensor:
 
     def test_num_params(self):
         """Test that the correct number of parameters is returned"""
-        p = np.array([0.5])
+        p = np.eye(4)
         X = qml.PauliX(0)
         Y = qml.Hermitian(p, wires=[1, 2])
         Z = qml.Hermitian(p, wires=[1, 2])
@@ -1237,7 +1237,7 @@ class TestTensor:
 
     def test_parameters(self):
         """Test that the correct nested list of parameters is returned"""
-        p = np.array([0.5])
+        p = np.eye(4)
         X = qml.PauliX(0)
         Y = qml.Hermitian(p, wires=[1, 2])
         t = Tensor(X, Y)