Repository: PennyLaneAI/pennylane (python)

Issue:
[BUG] Issue with expval of SparseH not building correct sized matrix
### Expected behavior

The following example should run to completion:
```python
import pennylane as qml
dev = qml.device("default.qubit", wires=4)

@qml.qnode(dev, diff_method="parameter-shift")
def circuit(p):
    qml.PauliX(0)
    return qml.expval(qml.SparseHamiltonian( qml.utils.sparse_hamiltonian( 1.1 * qml.PauliX(wires=0) @ qml.PauliZ(wires=2) ), wires=range(4),))
circuit()
```

### Actual behavior

The internal representation does not build a sparse H of the correct size for application to the dense statevector:

```python
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [1], line 8
      6     qml.PauliX(0)
      7     return qml.expval(qml.SparseHamiltonian( qml.utils.sparse_hamiltonian( 1.1 * qml.PauliX(wires=0) @ qml.PauliZ(wires=2) ), wires=range(4),))
----> 8 circuit()

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/pennylane/qnode.py:661, in QNode.__call__(self, *args, **kwargs)
    657             res = tuple(res)
    659     return res
--> 661 res = qml.execute(
    662     [self.tape],
    663     device=self.device,
    664     gradient_fn=self.gradient_fn,
    665     interface=self.interface,
    666     gradient_kwargs=self.gradient_kwargs,
    667     override_shots=override_shots,
    668     **self.execute_kwargs,
    669 )
    671 if autograd.isinstance(res, (tuple, list)) and len(res) == 1:
    672     # If a device batch transform was applied, we need to 'unpack'
    673     # the returned tuple/list to a float.
   (...)
    680     # TODO: find a more explicit way of determining that a batch transform
    681     # was applied.
    683     res = res[0]

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/pennylane/interfaces/execution.py:443, in execute(tapes, device, gradient_fn, interface, mode, gradient_kwargs, cache, cachesize, max_diff, override_shots, expand_fn, max_expansion, device_batch_transform)
    437 except ImportError as e:
    438     raise qml.QuantumFunctionError(
    439         f"{mapped_interface} not found. Please install the latest "
    440         f"version of {mapped_interface} to enable the '{mapped_interface}' interface."
    441     ) from e
--> 443 res = _execute(
    444     tapes, device, execute_fn, gradient_fn, gradient_kwargs, _n=1, max_diff=max_diff, mode=_mode
    445 )
    447 return batch_fn(res)

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/pennylane/interfaces/autograd.py:66, in execute(tapes, device, execute_fn, gradient_fn, gradient_kwargs, _n, max_diff, mode)
     60 # pylint misidentifies autograd.builtins as a dict
     61 # pylint: disable=no-member
     62 parameters = autograd.builtins.tuple(
     63     [autograd.builtins.list(t.get_parameters()) for t in tapes]
     64 )
---> 66 return _execute(
     67     parameters,
     68     tapes=tapes,
     69     device=device,
     70     execute_fn=execute_fn,
     71     gradient_fn=gradient_fn,
     72     gradient_kwargs=gradient_kwargs,
     73     _n=_n,
     74     max_diff=max_diff,
     75 )[0]

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/autograd/tracer.py:48, in primitive.<locals>.f_wrapped(*args, **kwargs)
     46     return new_box(ans, trace, node)
     47 else:
---> 48     return f_raw(*args, **kwargs)

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/pennylane/interfaces/autograd.py:110, in _execute(parameters, tapes, device, execute_fn, gradient_fn, gradient_kwargs, _n, max_diff)
     89 """Autodifferentiable wrapper around ``Device.batch_execute``.
     90 
     91 The signature of this function is designed to work around Autograd restrictions.
   (...)
    107 understand the consequences!
    108 """
    109 with qml.tape.Unwrap(*tapes):
--> 110     res, jacs = execute_fn(tapes, **gradient_kwargs)
    112 for i, r in enumerate(res):
    114     if any(
    115         m.return_type in (qml.measurements.Counts, qml.measurements.AllCounts)
    116         for m in tapes[i].measurements
    117     ):

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/pennylane/interfaces/execution.py:197, in cache_execute.<locals>.wrapper(tapes, **kwargs)
    193         return (res, []) if return_tuple else res
    195 else:
    196     # execute all unique tapes that do not exist in the cache
--> 197     res = fn(execution_tapes.values(), **kwargs)
    199 final_res = []
    201 for i, tape in enumerate(tapes):

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/pennylane/interfaces/execution.py:122, in cache_execute.<locals>.fn(tapes, **kwargs)
    120 def fn(tapes, **kwargs):  # pylint: disable=function-redefined
    121     tapes = [expand_fn(tape) for tape in tapes]
--> 122     return original_fn(tapes, **kwargs)

File /usr/lib/python3.10/contextlib.py:79, in ContextDecorator.__call__.<locals>.inner(*args, **kwds)
     76 @wraps(func)
     77 def inner(*args, **kwds):
     78     with self._recreate_cm():
---> 79         return func(*args, **kwds)

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/pennylane/_qubit_device.py:586, in QubitDevice.batch_execute(self, circuits)
    583     self.reset()
    585     # TODO: Insert control on value here
--> 586     res = self.execute(circuit)
    587     results.append(res)
    589 if self.tracker.active:

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/pennylane/_qubit_device.py:330, in QubitDevice.execute(self, circuit, **kwargs)
    328     results = self._collect_shotvector_results(circuit, counts_exist)
    329 else:
--> 330     results = self.statistics(circuit.observables, circuit=circuit)
    332 if not circuit.is_sampled:
    334     if len(circuit.measurements) == 1:

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/pennylane/_qubit_device.py:740, in QubitDevice.statistics(self, observables, shot_range, bin_size, circuit)
    736 for obs in observables:
    737     # Pass instances directly
    738     if obs.return_type is Expectation:
    739         # Appends a result of shape (num_bins,) if bin_size is not None, else a scalar
--> 740         results.append(self.expval(obs, shot_range=shot_range, bin_size=bin_size))
    742     elif obs.return_type is Variance:
    743         # Appends a result of shape (num_bins,) if bin_size is not None, else a scalar
    744         results.append(self.var(obs, shot_range=shot_range, bin_size=bin_size))

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/pennylane/devices/default_qubit.py:595, in DefaultQubit.expval(self, observable, shot_range, bin_size)
    583         res = qml.math.array(
    584             [
    585                 csr_matrix.dot(
   (...)
    590             ]
    591         )
    592     else:
    593         res = csr_matrix.dot(
    594             csr_matrix(self._conj(state)),
--> 595             csr_matrix.dot(Hmat, csr_matrix(state[..., None])),
    596         ).toarray()[0]
    598 if observable.name == "Hamiltonian":
    599     res = qml.math.squeeze(res)

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/scipy/sparse/_base.py:416, in spmatrix.dot(self, other)
    414     return self * other
    415 else:
--> 416     return self @ other

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/scipy/sparse/_base.py:630, in spmatrix.__matmul__(self, other)
    627 if isscalarlike(other):
    628     raise ValueError("Scalar operands are not allowed, "
    629                      "use '*' instead")
--> 630 return self._mul_dispatch(other)

File /tmp/pennylane-lightning-gpu/pyenv/lib/python3.10/site-packages/scipy/sparse/_base.py:540, in spmatrix._mul_dispatch(self, other)
    538 if issparse(other):
    539     if self.shape[1] != other.shape[0]:
--> 540         raise ValueError('dimension mismatch')
    541     return self._mul_sparse_matrix(other)
    543 # If it's a list or whatever, treat it like a matrix

ValueError: dimension mismatch
```

### Additional information

This is in the most recent released and master branches of PennyLane.

### Source code

_No response_

### Tracebacks

_No response_

### 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: /tmp/pennylane-lightning-gpu/pyenv/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.19.13-arch1-1-x86_64-with-glibc2.36
Python version:          3.10.7
Numpy version:           1.23.4
Scipy version:           1.9.3
Installed devices:
- lightning.qubit (PennyLane-Lightning-0.26.1)
- 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:
Add shape validation to `SparseHamiltonian`
Fixes #3272 .

Currently, `SparseHamiltonian` does not perform any validation on the shape of the matrix.  This leads to confusing errors downstream if the shape does not match the wires.

This PR just adds a simple check on the shape, and raises an error if the sparse matrix is not the right shape.

This is a rapid check, so I do not see any problems adding it in.

Patch:
diff --git a/doc/releases/changelog-0.27.0.md b/doc/releases/changelog-0.27.0.md
index e7f2db581..853c0da7a 100644
--- a/doc/releases/changelog-0.27.0.md
+++ b/doc/releases/changelog-0.27.0.md
@@ -592,6 +592,9 @@
 
 <h3>Bug fixes</h3>
 
+* `qml.SparseHamiltonian` now validates the size of the input matrix.
+  [(#3278)](https://github.com/PennyLaneAI/pennylane/pull/3278)
+
 * Fixed a bug where `qml.sample()` and `qml.counts()` would output incorrect results when mixed with measurements whose
   operators do not qubit-wise commute with computational basis projectors.
   [(#3207)](https://github.com/PennyLaneAI/pennylane/pull/3207)
diff --git a/pennylane/ops/qubit/observables.py b/pennylane/ops/qubit/observables.py
index 7d1e50ca8..5138d7ec3 100644
--- a/pennylane/ops/qubit/observables.py
+++ b/pennylane/ops/qubit/observables.py
@@ -252,6 +252,11 @@ class SparseHamiltonian(Observable):
         if not isinstance(H, csr_matrix):
             raise TypeError("Observable must be a scipy sparse csr_matrix.")
         super().__init__(H, wires=wires, do_queue=do_queue, id=id)
+        mat_len = 2 ** len(self.wires)
+        if H.shape != (mat_len, mat_len):
+            raise ValueError(
+                f"Sparse Matrix must be of shape ({mat_len}, {mat_len}). Got {H.shape}."
+            )
 
     def label(self, decimals=None, base_label=None, cache=None):
         return super().label(decimals=decimals, base_label=base_label or "𝓗", cache=cache)


Test patch:
diff --git a/tests/ops/functions/test_eigvals.py b/tests/ops/functions/test_eigvals.py
index a0f4a9074..732b38e6c 100644
--- a/tests/ops/functions/test_eigvals.py
+++ b/tests/ops/functions/test_eigvals.py
@@ -203,23 +203,24 @@ class TestSingleOperation:
     )
     def test_sparse_hamiltonian(self, row, col, dat, val_ref):
         """Test that the eigenvalues of a sparse Hamiltonian are correctly returned"""
-        # N x N matrix with N = 20
-        h_sparse = scipy.sparse.csr_matrix((dat, (row, col)), shape=(len(row), len(col)))
-        h_sparse = qml.SparseHamiltonian(h_sparse, wires=all)
+        # N x N matrix with N = 16
+        h_mat = scipy.sparse.csr_matrix((dat, (row, col)), shape=(16, 16))
+        h_sparse = qml.SparseHamiltonian(h_mat, wires=range(4))
+
+        dense_mat = h_mat.todense()
+        dense_eigvals = np.sort(np.linalg.eigvals(dense_mat))
 
         # k = 1  (< N-1) scipy.sparse.linalg is used:
         val_groundstate = qml.eigvals(h_sparse, k=1)
-        # k = 18  (< N-1) scipy.sparse.linalg is used:
-        val_n_sparse = qml.eigvals(h_sparse, k=len(val_ref) - 2)
-        # k = 20 (> N-1) qml.math.linalg is used:
-        val_all = qml.eigvals(h_sparse, k=len(val_ref))
-        # k = 19 (= N-1) qml.math.linalg is used:
-        val_n_dense = qml.eigvals(h_sparse, k=len(val_ref) - 1)
-
-        assert np.allclose(val_groundstate, val_ref[0])
-        assert np.allclose(np.sort(val_n_sparse), val_ref[0:18])
-        assert np.allclose(np.sort(val_all), val_ref)
-        assert np.allclose(np.sort(val_n_dense), val_ref)
+        assert np.allclose(val_groundstate, dense_eigvals[0])
+
+        # k = 14  (< N-1) scipy.sparse.linalg is used:
+        val_n_sparse = np.sort(qml.eigvals(h_sparse, k=14))
+        assert np.allclose(val_n_sparse, dense_eigvals[0:-2])
+
+        # k = 16 (> N-1) qml.math.linalg is used:
+        val_all = np.sort(qml.eigvals(h_sparse, k=16))
+        assert np.allclose(val_all, dense_eigvals)
 
 
 class TestMultipleOperations:
diff --git a/tests/ops/qubit/test_sparse.py b/tests/ops/qubit/test_sparse.py
index 8ee3eb514..da13e17b9 100644
--- a/tests/ops/qubit/test_sparse.py
+++ b/tests/ops/qubit/test_sparse.py
@@ -62,6 +62,14 @@ class TestSparse:
         H = qml.SparseHamiltonian(csr_matrix(np.array([[1, 0], [-1.5, 0]])), 1)
         assert H.label() == "𝓗"
 
+    def test_valueeror_shape(self):
+        """Test that iniatilization raises a ValueError if the shape is not correct for the number of wires."""
+        mat = csr_matrix([[1, 0], [0, 1]])
+        with pytest.raises(
+            ValueError, match=r"Sparse Matrix must be of shape \(4, 4\). Got \(2, 2\)."
+        ):
+            qml.SparseHamiltonian(mat, wires=(0, 1))
+
     @pytest.mark.parametrize("sparse_hamiltonian", SPARSE_HAMILTONIAN_TEST_DATA)
     def test_sparse_typeerror(self, sparse_hamiltonian):
         """Test that the SparseHamiltonian class raises a TypeError on incorrect inputs."""
@@ -69,12 +77,12 @@ class TestSparse:
         sparse_hamiltonian = coo_matrix(mat)
 
         with pytest.raises(TypeError, match="Observable must be a scipy sparse csr_matrix"):
-            qml.SparseHamiltonian(sparse_hamiltonian)
+            qml.SparseHamiltonian(sparse_hamiltonian, wires=0)
 
     @pytest.mark.parametrize("sparse_hamiltonian", SPARSE_HAMILTONIAN_TEST_DATA)
     def test_sparse_matrix(self, sparse_hamiltonian, tol):
         """Test that the matrix property of the SparseHamiltonian class returns the correct matrix."""
-        num_wires = len(sparse_hamiltonian[0])
+        num_wires = int(np.log2(len(sparse_hamiltonian[0])))
         sparse_hamiltonian_csr = csr_matrix(sparse_hamiltonian)
         res_dynamic = qml.SparseHamiltonian(
             sparse_hamiltonian_csr, range(num_wires)
@@ -88,7 +96,7 @@ class TestSparse:
     @pytest.mark.parametrize("sparse_hamiltonian", SPARSE_HAMILTONIAN_TEST_DATA)
     def test_matrix(self, sparse_hamiltonian, tol):
         """Test that the matrix property of the SparseHamiltonian class returns the correct matrix."""
-        num_wires = len(sparse_hamiltonian[0])
+        num_wires = int(np.log2(len(sparse_hamiltonian[0])))
         sparse_hamiltonian_csr = csr_matrix(sparse_hamiltonian)
         res_dynamic = qml.SparseHamiltonian(sparse_hamiltonian_csr, range(num_wires)).matrix()
         res_static = qml.SparseHamiltonian.compute_matrix(sparse_hamiltonian_csr)