Add functions to support a sparse matrix observable

.github/CHANGELOG.md

@@ -23,6 +23,29 @@ Ashish Panigrahi
 
 <h3>New features since last release</h3>
 
+* Added a sparse hamiltonian observable and the functionality to support computing its expectation
+  value. [(#1398)](https://github.com/PennyLaneAI/pennylane/pull/1398)
+
+  For example, the expectation value of a sparse hamiltonian, an identity matrix in this example,
+  can be computed as:
+
+  ```python
+  dev = qml.device("default.qubit", wires=2)
+
+  @qml.qnode(dev, diff_method="parameter-shift")
+  def circuit(param, H):
+      qml.PauliX(0)
+      qml.SingleExcitation(param, wires = [0, 1])    
+      return qml.expval(qml.SparseHamiltonian(H, wires=[]))
+
+  >>> print(circuit([0.5], scipy.sparse.eye(4).tocoo()))
+  0.9999999999999999
+  ```
+
+  The expectation value of the sparse hamiltonian is computed directly, which leads to executions
+  that are faster by orders of magnitude. Note that "parameter-shift" is the only differentiation
+  method that is currently supported when the observable is a sparse hamiltonian.
+
 * Added functionality to compute the sparse matrix representation of a `qml.Hamiltonian` object.
   [(#1394)](https://github.com/PennyLaneAI/pennylane/pull/1394)
 

pennylane/devices/default_qubit.py

@@ -23,6 +23,7 @@
 from string import ascii_letters as ABC
 
 import numpy as np
+from scipy.sparse import coo_matrix
 
 from pennylane import QubitDevice, DeviceError, QubitStateVector, BasisState
 from pennylane.operation import DiagonalOperation
@@ -135,7 +136,16 @@ class DefaultQubit(QubitDevice):
         "QubitSum",
     }
 
-    observables = {"PauliX", "PauliY", "PauliZ", "Hadamard", "Hermitian", "Identity", "Projector"}
+    observables = {
+        "PauliX",
+        "PauliY",
+        "PauliZ",
+        "Hadamard",
+        "Hermitian",
+        "Identity",
+        "Projector",
+        "SparseHamiltonian",
+    }
 
     def __init__(self, wires, *, shots=None, cache=0, analytic=None):
         super().__init__(wires, shots, cache=cache, analytic=analytic)
@@ -437,6 +447,38 @@ def _apply_phase(self, state, axes, parameters, inverse=False):
         phase = self._conj(parameters) if inverse else parameters
         return self._stack([state[sl_0], phase * state[sl_1]], axis=axes[0])
 
+    def expval(self, observable, shot_range=None, bin_size=None):
+        """Returns the expectation value of a Hamiltonian observable. When the observable is a
+         ``SparseHamiltonian`` object, the expectation value is computed directly for the full
+         Hamiltonian, which leads to faster execution.
+
+        Args:
+            observable (~.Observable): a PennyLane observable
+            shot_range (tuple[int]): 2-tuple of integers specifying the range of samples
+                to use. If not specified, all samples are used.
+            bin_size (int): Divides the shot range into bins of size ``bin_size``, and
+                returns the measurement statistic separately over each bin. If not
+                provided, the entire shot range is treated as a single bin.
+
+        Returns:
+            float: returns the expectation value of the observable
+        """
+        if observable.name == "SparseHamiltonian" and self.shots is not None:
+            raise DeviceError("SparseHamiltonian must be used with shots=None")
+
+        if observable.name == "SparseHamiltonian" and self.shots is None:
+
+            ev = coo_matrix.dot(
+                coo_matrix(self.state),
+                coo_matrix.dot(
+                    observable.matrix, coo_matrix(self.state.reshape(len(self.state), 1))
+                ),
+            )
+
+            return np.real(ev.toarray()[0])
+
+        return super().expval(observable, shot_range=shot_range, bin_size=bin_size)
+
     def _get_unitary_matrix(self, unitary):  # pylint: disable=no-self-use
         """Return the matrix representing a unitary operation.
 

pennylane/ops/qubit.py

@@ -21,6 +21,7 @@
 # pylint:disable=abstract-method,arguments-differ,protected-access
 import math
 import numpy as np
+import scipy
 from scipy.linalg import block_diag
 
 import pennylane as qml
@@ -3085,6 +3086,41 @@ def diagonalizing_gates(self):
         return []
 
 
+class SparseHamiltonian(Observable):
+    r"""SparseHamiltonian(H)
+    A Hamiltonian represented directly as a sparse matrix in coordinate list (COO) format.
+
+    **Details:**
+
+    * Number of wires: Any
+    * Number of parameters: 1
+    * Gradient recipe: None
+
+    Args:
+        H (coo_matrix): a sparse matrix in scipy coordinate list (COO) format with
+            dimension :math:`(2^n, 2^n)`, where :math:`n` is the number of wires
+    """
+    num_wires = AnyWires
+    num_params = 1
+    par_domain = None
+    grad_method = None
+
+    wires = Wires([0])
+
+    def __init__(self, *params, wires=wires):
+        super().__init__(*params, wires=wires)
+
+    @classmethod
+    def _matrix(cls, *params):
+        A = params[0]
+        if not isinstance(A, scipy.sparse.coo_matrix):
+            raise TypeError("Observable must be a scipy sparse coo_matrix.")
+        return A
+
+    def diagonalizing_gates(self):
+        return []
+
+
 # =============================================================================
 # Arithmetic
 # =============================================================================
@@ -3331,7 +3367,7 @@ def adjoint(self):
 }
 
 
-obs = {"Hadamard", "PauliX", "PauliY", "PauliZ", "Hermitian", "Projector"}
+obs = {"Hadamard", "PauliX", "PauliY", "PauliZ", "Hermitian", "Projector", "SparseHamiltonian"}
 
 
 __all__ = list(ops | obs)

pennylane/qnode.py

@@ -538,13 +538,23 @@ def construct(self, args, kwargs):
             )
 
         for obj in self.qtape.operations + self.qtape.observables:
+
             if getattr(obj, "num_wires", None) is qml.operation.WiresEnum.AllWires:
                 # check here only if enough wires
                 if len(obj.wires) != self.device.num_wires:
                     raise qml.QuantumFunctionError(
                         "Operator {} must act on all wires".format(obj.name)
                     )
 
+            if (
+                isinstance(obj, qml.ops.qubit.SparseHamiltonian)
+                and self.diff_method != "parameter-shift"
+            ):
+                raise qml.QuantumFunctionError(
+                    "SparseHamiltonian observable must be used with the parameter-shift"
+                    " differentiation method"
+                )
+
         # pylint: disable=protected-access
         obs_on_same_wire = len(self.qtape._obs_sharing_wires) > 0
         ops_not_supported = any(

tests/devices/test_default_qubit.py

@@ -21,6 +21,7 @@
 
 import pytest
 import pennylane as qml
+from scipy.sparse import coo_matrix
 from pennylane import numpy as np, DeviceError
 from pennylane.devices.default_qubit import _get_slice, DefaultQubit
 from pennylane.wires import Wires, WireError
@@ -873,6 +874,84 @@ def circuit():
         # an estimated variance an an analytically calculated one
         assert expval != 0.0
 
+    H_row = np.array([0, 1, 2, 3, 3, 4, 5, 6, 6, 7, 8, 9, 9, 10, 11, 12, 12, 13, 14, 15])
+    H_col = np.array([0, 1, 2, 3, 12, 4, 5, 6, 9, 7, 8, 6, 9, 10, 11, 3, 12, 13, 14, 15])
+    H_data = np.array(
+        [
+            0.72004228 + 0.0j,
+            0.24819411 + 0.0j,
+            0.24819411 + 0.0j,
+            0.47493347 + 0.0j,
+            0.18092703 + 0.0j,
+            -0.5363422 + 0.0j,
+            -0.52452263 + 0.0j,
+            -0.34359561 + 0.0j,
+            -0.18092703 + 0.0j,
+            0.3668115 + 0.0j,
+            -0.5363422 + 0.0j,
+            -0.18092703 + 0.0j,
+            -0.34359561 + 0.0j,
+            -0.52452263 + 0.0j,
+            0.3668115 + 0.0j,
+            0.18092703 + 0.0j,
+            -1.11700225 + 0.0j,
+            -0.44058791 + 0.0j,
+            -0.44058791 + 0.0j,
+            0.93441396 + 0.0j,
+        ]
+    )
+    H_hydrogen = coo_matrix((H_data, (H_row, H_col)), shape=(16, 16)).toarray()
+
+    @pytest.mark.parametrize(
+        "qubits, operations, hamiltonian, expected_output",
+        [
+            (1, [], np.array([[1.0, 0.0], [0.0, 1.0]]), 1.0),
+            (
+                2,
+                [qml.PauliX(0), qml.PauliY(1)],
+                np.array(
+                    [
+                        [1.0, 0.0, 0.0, 0.0],
+                        [0.0, 1.0, 0.0, 0.0],
+                        [0.0, 0.0, 1.0, 0.0],
+                        [0.0, 0.0, 0.0, 1.0],
+                    ]
+                ),
+                -1.0,
+            ),
+            (
+                4,
+                [
+                    qml.PauliX(0),
+                    qml.PauliX(1),
+                    qml.DoubleExcitation(0.22350048065138242, wires=[0, 1, 2, 3]),
+                ],
+                H_hydrogen,
+                -1.1373060481,
+            ),
+        ],
+    )
+    def test_sparse_hamiltonian_expval(self, qubits, operations, hamiltonian, expected_output, tol):
+        """Test that expectation values of sparse hamiltonians are properly calculated."""
+
+        hamiltonian = coo_matrix(hamiltonian)
+
+        dev = qml.device("default.qubit", wires=qubits, shots=None)
+        dev.apply(operations)
+        expval = dev.expval(qml.SparseHamiltonian(hamiltonian))[0]
+
+        assert np.allclose(expval, expected_output, atol=tol, rtol=0)
+
+    def test_sparse_expval_error(self):
+        """Test that the expval function raises a DeviceError when the observable is
+        SparseHamiltonian and finite shots is requested."""
+        hamiltonian = coo_matrix(np.array([[1.0, 0.0], [0.0, 1.0]]))
+
+        dev = qml.device("default.qubit", wires=1, shots=1)
+
+        with pytest.raises(DeviceError, match="SparseHamiltonian must be used with shots=None"):
+            dev.expval(qml.SparseHamiltonian(hamiltonian))[0]
+
 
 class TestVar:
     """Tests that variances are properly calculated."""

tests/ops/test_qubit_ops.py

@@ -22,6 +22,7 @@
 from numpy.linalg import multi_dot
 from scipy.stats import unitary_group
 from scipy.linalg import expm
+from scipy.sparse import coo_matrix, csr_matrix
 
 import pennylane as qml
 from pennylane.wires import Wires
@@ -93,6 +94,9 @@
     (np.array([1, 0, 1])),
 ]
 
+# Testing SparseHamiltonian observable.
+SPARSE_HAMILTONIAN_TEST_DATA = [(np.array([[1, 0], [-1.5, 0]])), (np.eye(4))]
+
 
 @pytest.mark.usefixtures("tear_down_hermitian")
 class TestObservables:
@@ -394,6 +398,73 @@ def circuit(basis_state):
             circuit(basis_state)
 
 
+class TestSparse:
+    """Tests for sparse hamiltonian observable"""
+
+    @pytest.mark.parametrize("sparse_hamiltonian", SPARSE_HAMILTONIAN_TEST_DATA)
+    def test_sparse_diagonalization(self, sparse_hamiltonian):
+        """Test that the diagonalizing_gates property of the SparseHamiltonian class returns empty."""
+        num_wires = len(sparse_hamiltonian[0])
+        sparse_hamiltonian = coo_matrix(sparse_hamiltonian)
+        diag_gates = qml.SparseHamiltonian(sparse_hamiltonian).diagonalizing_gates()
+
+        assert diag_gates == []
+
+    @pytest.mark.parametrize("sparse_hamiltonian", SPARSE_HAMILTONIAN_TEST_DATA)
+    def test_sparse_typeerror(self, sparse_hamiltonian):
+        """Test that the matrix property of the SparseHamiltonian class raises a TypeError on incorrect inputs."""
+        num_wires = len(sparse_hamiltonian[0])
+        sparse_hamiltonian = csr_matrix(sparse_hamiltonian)
+
+        dev = qml.device("default.qubit", wires=num_wires)
+
+        @qml.qnode(dev, diff_method="parameter-shift")
+        def circuit(sparse_hamiltonian, num_wires):
+            obs = qml.SparseHamiltonian(sparse_hamiltonian)
+            return qml.expval(obs)
+
+        with pytest.raises(TypeError, match="Observable must be a scipy sparse coo_matrix"):
+            circuit(sparse_hamiltonian, num_wires)
+
+    @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])
+        sparse_hamiltonian = coo_matrix(sparse_hamiltonian)
+        returned_matrix = qml.SparseHamiltonian(sparse_hamiltonian).matrix
+        assert np.allclose(
+            returned_matrix.toarray(), sparse_hamiltonian.toarray(), atol=tol, rtol=0
+        )
+
+    def test_sparse_gradient(self, tol):
+        """Tests that gradients are computed correctly for a SparseHamiltonian observable."""
+        dev = qml.device("default.qubit", wires=2, shots=None)
+
+        @qml.qnode(dev, diff_method="parameter-shift")
+        def circuit(param):
+            qml.RX(param, wires=0)
+            return qml.expval(qml.SparseHamiltonian(coo_matrix(np.eye(4))))
+
+        assert np.allclose(qml.grad(circuit)([0.5]), -0.47942554, atol=tol, rtol=0)
+
+    def test_sparse_diffmethod_error(self):
+        """Test that an error is raised when the observable is SparseHamiltonian and the
+        differentiation method is not parameter-shift."""
+        dev = qml.device("default.qubit", wires=2, shots=None)
+
+        @qml.qnode(dev, diff_method="backprop")
+        def circuit(param):
+            qml.RX(param, wires=0)
+            return qml.expval(qml.SparseHamiltonian(coo_matrix(np.eye(4))))
+
+        with pytest.raises(
+            qml.QuantumFunctionError,
+            match="SparseHamiltonian observable must be"
+            " used with the parameter-shift differentiation method",
+        ):
+            qml.grad(circuit)([0.5])
+
+
 # Non-parametrized operations and their matrix representation
 NON_PARAMETRIZED_OPERATIONS = [
     (qml.CNOT, CNOT),