Repository: qiboteam/qibo

Issue:
Energy callback for sparse hamiltonians
I am currently working in a field theory project where we must compute the energy of a ~20 qubit hamiltonian. The problem is that the statevector is 20-qubit long, but the hamiltonian takes the same memory as a 40-qubit state, which does not fit in any RAM. I have observed that `cirq` and `openfermion` are capable to deal with this kind of hamiltonians by letting the matrix be sparse.  This in fact can be done by qibo obtaining the expectation value of a sparse matrix as 

```
E = openfermion.expectation(sparse_matrix, state)
```

or equivalently
```
E = np.conj(state) @ sparse_matrix @ state
```


Now, the problem is to configure a callback for the energy from a sparse matrix as our hamiltonian, for example

```
from qibo import hamiltonians, callbacks
ham = hamiltonians.Hamiltonian(qubits, sparse_matrix) 
energy_callback = callbacks.Energy(ham)

C = Circuit(qubits)
C.add(operations_1)
adiabatic_circuit.add(gates.CallbackGate(energy_callback))
C.add(operations_2)

C.execute()
```


Is it possible? When I run a similar code my error is
`  TypeError: Matrix of invalid type <class 'scipy.sparse.csc.csc_matrix'> given during Hamiltonian initialization `


Discussion:

Patch:
diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst
index cbe909297..ae7441cbd 100644
--- a/doc/source/api-reference/qibo.rst
+++ b/doc/source/api-reference/qibo.rst
@@ -581,9 +581,19 @@ of the Hamiltonian operator in the computational basis. This matrix has size
 ``(2 ** nqubits, 2 ** nqubits)`` and therefore its construction is feasible
 only when number of qubits is small.
 
-.. autoclass:: qibo.abstractions.hamiltonians.MatrixHamiltonian
-    :members:
-    :member-order: bysource
+Alternatively, the user can construct this Hamiltonian using a sparse matrices.
+Sparse matrices from the
+`scipy.sparse <https://docs.scipy.org/doc/scipy/reference/sparse.html>`_
+module are supported by the numpy and qibojit backends while the
+`tf.sparse <https://www.tensorflow.org/api_docs/python/tf/sparse>_` can be
+used for tensorflow and qibotf. Scipy sparse matrices support algebraic
+operations (addition, subtraction, scalar multiplication), linear algebra
+operations (eigenvalues, eigenvectors, matrix exponentiation) and
+multiplication to dense or other sparse matrices. All these properties are
+inherited by :class:`qibo.core.hamiltonians.Hamiltonian` objects created
+using sparse matrices. Tensorflow sparse matrices support only multiplication
+to dense matrices. Both backends support calculating Hamiltonian expectation
+values using a sparse Hamiltonian matrix.
 
 .. autoclass:: qibo.core.hamiltonians.Hamiltonian
     :members:
diff --git a/src/qibo/abstractions/hamiltonians.py b/src/qibo/abstractions/hamiltonians.py
index 58f058036..4e4e5b2c2 100644
--- a/src/qibo/abstractions/hamiltonians.py
+++ b/src/qibo/abstractions/hamiltonians.py
@@ -23,20 +23,36 @@ class AbstractHamiltonian(ABC):
         self._nqubits = n
 
     @abstractmethod
-    def eigenvalues(self): # pragma: no cover
-        """Computes the eigenvalues for the Hamiltonian."""
+    def eigenvalues(self, k=6): # pragma: no cover
+        """Computes the eigenvalues for the Hamiltonian.
+
+        Args:
+            k (int): Number of eigenvalues to calculate if the Hamiltonian
+                was created using a sparse matrix. This argument is ignored
+                if the Hamiltonian was created using a dense matrix.
+                See :meth:`qibo.backends.abstract.AbstractBackend.eigvalsh` for 
+                more details.
+        """
         raise_error(NotImplementedError)
 
     @abstractmethod
-    def eigenvectors(self): # pragma: no cover
-        """Computes a tensor with the eigenvectors for the Hamiltonian."""
+    def eigenvectors(self, k=6): # pragma: no cover
+        """Computes a tensor with the eigenvectors for the Hamiltonian.
+
+        Args:
+            k (int): Number of eigenvalues to calculate if the Hamiltonian
+                was created using a sparse matrix. This argument is ignored
+                if the Hamiltonian was created using a dense matrix.
+                See :meth:`qibo.backends.abstract.AbstractBackend.eigh` for
+                more details.
+        """
         raise_error(NotImplementedError)
 
     def ground_state(self):
         """Computes the ground state of the Hamiltonian.
 
-        Uses the ``eigenvectors`` method and returns the lowest energy
-        eigenvector.
+        Uses :meth:`qibo.abstractions.hamiltonians.AbstractHamiltonian.eigenvectors`
+        and returns eigenvector corresponding to the lowest energy.
         """
         return self.eigenvectors()[:, 0]
 
@@ -99,16 +115,7 @@ class AbstractHamiltonian(ABC):
 
 
 class MatrixHamiltonian(AbstractHamiltonian):
-    """Abstract Hamiltonian based on full matrix representation.
-
-    Args:
-        nqubits (int): number of quantum bits.
-        matrix (np.ndarray): Matrix representation of the Hamiltonian in the
-            computational basis as an array of shape ``(2 ** nqubits, 2 ** nqubits)``.
-        numpy (bool): If ``True`` the Hamiltonian is created using numpy as the
-            calculation backend, otherwise the selected backend is used.
-            Default option is ``numpy = False``.
-    """
+    """Abstract Hamiltonian based on full matrix representation."""
 
     def __init__(self, nqubits, matrix=None):
         super().__init__()
@@ -120,7 +127,11 @@ class MatrixHamiltonian(AbstractHamiltonian):
 
     @property
     def matrix(self):
-        """Returns the full ``(2 ** nqubits, 2 ** nqubits)`` matrix representation."""
+        """Returns the full matrix representation.
+
+        Can be a dense ``(2 ** nqubits, 2 ** nqubits)`` array or a sparse
+        matrix, depending on how the Hamiltonian was created.
+        """
         return self._matrix
 
     @matrix.setter
@@ -174,11 +185,11 @@ class SymbolicHamiltonian(AbstractHamiltonian):
         """Returns the full ``(2 ** nqubits, 2 ** nqubits)`` matrix representation."""
         return self.dense.matrix
 
-    def eigenvalues(self):
-        return self.dense.eigenvalues()
+    def eigenvalues(self, k=6):
+        return self.dense.eigenvalues(k)
 
-    def eigenvectors(self):
-        return self.dense.eigenvectors()
+    def eigenvectors(self, k=6):
+        return self.dense.eigenvectors(k)
 
     def ground_state(self):
         if self._ground_state is None:
diff --git a/src/qibo/backends/abstract.py b/src/qibo/backends/abstract.py
index b1d2b3fff..5c5c3a7d4 100644
--- a/src/qibo/backends/abstract.py
+++ b/src/qibo/backends/abstract.py
@@ -34,6 +34,7 @@ class AbstractBackend(ABC):
         self.Tensor = None
         self.random = None
         self.newaxis = None
+        self.sparse = None
         self.oom_error = None
         self.optimization = None
         self.supports_multigpu = False
@@ -198,6 +199,11 @@ class AbstractBackend(ABC):
         """Casts tensor to the given dtype."""
         raise_error(NotImplementedError)
 
+    @abstractmethod
+    def issparse(self, x): # pragma: no cover
+        """Checks if the given tensor is sparse."""
+        raise_error(NotImplementedError)
+
     @abstractmethod
     def diag(self, x, dtype='DTYPECPX'): # pragma: no cover
         raise_error(NotImplementedError)
@@ -337,6 +343,11 @@ class AbstractBackend(ABC):
         """Sum of tensor elements."""
         raise_error(NotImplementedError)
 
+    @abstractmethod
+    def dot(self, x, y): # pragma: no cover
+        """Dot product of two tensors."""
+        raise_error(NotImplementedError)
+
     @abstractmethod
     def matmul(self, x, y): # pragma: no cover
         """Matrix multiplication of two tensors."""
@@ -373,13 +384,30 @@ class AbstractBackend(ABC):
         raise_error(NotImplementedError)
 
     @abstractmethod
-    def eigh(self, x): # pragma: no cover
-        """Hermitian matrix eigenvalues and eigenvectors."""
+    def eigh(self, x, k=6): # pragma: no cover
+        """Hermitian matrix eigenvalues and eigenvectors.
+
+        Args:
+            x: Tensor to calculate the eigenvectors of.
+            k (int): Number of eigenvectors to calculate if a sparse matrix is
+                given. The eigenvectors corresponding to the k-th algebraically
+                smallest eigenvalues are calculated.
+                This argument is ignored if the given tensor is not sparse and
+                all eigenvectors are calculated.
+        """
         raise_error(NotImplementedError)
 
     @abstractmethod
-    def eigvalsh(self, x): # pragma: no cover
-        """Hermitian matrix eigenvalues."""
+    def eigvalsh(self, x, k=6): # pragma: no cover
+        """Hermitian matrix eigenvalues.
+
+        Args:
+            x: Tensor to calculate the eigenvalues.
+            k (int): Number of eigenvalues to calculate if a sparse matrix is
+                given. The k-th algebraically smallest eigenvalues are calculated.
+                This argument is ignored if the given tensor is not sparse and
+                all eigenvalues are calculated.
+        """
         raise_error(NotImplementedError)
 
     @abstractmethod
diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index f8e5e2cfd..5670b84c3 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -38,6 +38,8 @@ class NumpyBackend(abstract.AbstractBackend):
         self.Tensor = np.ndarray
         self.random = np.random
         self.newaxis = np.newaxis
+        from scipy import sparse
+        self.sparse = sparse
         self.oom_error = MemoryError
         self.optimization = None
         self.cpu_devices = ["/CPU:0"]
@@ -60,6 +62,8 @@ class NumpyBackend(abstract.AbstractBackend):
         return None
 
     def to_numpy(self, x):
+        if self.issparse(x):
+            return x.toarray()
         return x
 
     def to_complex(self, re, img): # pragma: no cover
@@ -70,8 +74,13 @@ class NumpyBackend(abstract.AbstractBackend):
             dtype = self.dtypes(dtype)
         if isinstance(x, self.backend.ndarray):
             return x.astype(dtype, copy=False)
+        elif self.issparse(x):
+            return x.astype(dtype=dtype, copy=False)
         return self.backend.array(x, dtype=dtype)
 
+    def issparse(self, x):
+        return self.sparse.issparse(x)
+
     def diag(self, x, dtype='DTYPECPX'):
         if isinstance(dtype, str):
             dtype = self.dtypes(dtype)
@@ -161,12 +170,18 @@ class NumpyBackend(abstract.AbstractBackend):
         return self.backend.trace(x)
 
     def expm(self, x):
-        from scipy import linalg
-        return linalg.expm(x)
+        if self.issparse(x):
+            from scipy.sparse.linalg import expm
+        else:
+            from scipy.linalg import expm
+        return expm(x)
 
     def sum(self, x, axis=None):
         return self.backend.sum(x, axis=axis)
 
+    def dot(self, x, y):
+        return x.dot(y)
+
     def matmul(self, x, y):
         return self.backend.matmul(x, y)
 
@@ -188,10 +203,19 @@ class NumpyBackend(abstract.AbstractBackend):
     def inv(self, x):
         return self.backend.linalg.inv(x)
 
-    def eigh(self, x):
+    def eigh(self, x, k=6):
+        if self.issparse(x):
+            if k < x.shape[0]:
+                from scipy.sparse.linalg import eigsh
+                return eigsh(x, k=k, which='SA')
+            x = self.to_numpy(x)
         return self.backend.linalg.eigh(x)
 
-    def eigvalsh(self, x):
+    def eigvalsh(self, x, k=6):
+        if self.issparse(x):
+            log.warning("Calculating sparse matrix eigenvectors because "
+                        "sparse modules do not provide ``eigvals`` method.")
+            return self.eigh(x, k=k)[0]
         return self.backend.linalg.eigvalsh(x)
 
     def less(self, x, y):
@@ -428,4 +452,6 @@ class NumpyBackend(abstract.AbstractBackend):
         return self.reshape(self.transpose(pieces, order), state.shape)
 
     def assert_allclose(self, value, target, rtol=1e-7, atol=0.0):
+        value = self.to_numpy(value)
+        target = self.to_numpy(target)
         self.np.testing.assert_allclose(value, target, rtol=rtol, atol=atol)
diff --git a/src/qibo/backends/tensorflow.py b/src/qibo/backends/tensorflow.py
index 21a5f7b54..e8f9ee920 100644
--- a/src/qibo/backends/tensorflow.py
+++ b/src/qibo/backends/tensorflow.py
@@ -70,6 +70,7 @@ class TensorflowBackend(NumpyBackend):
         self.Tensor = tf.Tensor
         self.random = tf.random
         self.newaxis = tf.newaxis
+        self.sparse = tf.sparse
         from tensorflow.python.framework import errors_impl  # pylint: disable=E0611,E0401
         self.oom_error = errors_impl.ResourceExhaustedError
         self.optimization = Optimization()
@@ -77,7 +78,6 @@ class TensorflowBackend(NumpyBackend):
         # seed to use in the measurement frequency custom op
         self._seed = None
         # seed can be modified using ``K.set_seed``
-
         self.supports_gradients = True
 
     def test_regressions(self, name):
@@ -99,9 +99,15 @@ class TensorflowBackend(NumpyBackend):
         AbstractBackend.set_threads(self, nthreads)
 
     def to_numpy(self, x):
-        if isinstance(x, self.np.ndarray):
+        if isinstance(x, self.numeric_types):
             return x
-        return x.numpy()
+        elif self.issparse(x):
+            from scipy.sparse import coo_matrix
+            idx = self.np.array(x.indices)
+            values = self.np.array(x.values)
+            shape = tuple(x.shape)
+            return coo_matrix((values, (idx[:, 0], idx[:, 1])), shape=shape)
+        return self.np.array(x, copy=False)
 
     def to_complex(self, re, img):
         return self.backend.complex(re, img)
@@ -112,8 +118,13 @@ class TensorflowBackend(NumpyBackend):
         if isinstance(x, self.np.ndarray):
             dtypestr = dtype.__repr__().split(".")[1]
             x = x.astype(getattr(self.np, dtypestr))
+        elif self.issparse(x):
+            return x
         return self.backend.cast(x, dtype=dtype)
 
+    def issparse(self, x):
+        return isinstance(x, self.sparse.SparseTensor)
+
     def diag(self, x, dtype='DTYPECPX'):
         if isinstance(dtype, str):
             dtype = self.dtypes(dtype)
@@ -166,6 +177,15 @@ class TensorflowBackend(NumpyBackend):
     def sum(self, x, axis=None):
         return self.backend.reduce_sum(x, axis=axis)
 
+    def dot(self, x, y):
+        if self.issparse(x) or self.issparse(y):
+            if self.issparse(x) and self.issparse(y):
+                raise_error(NotImplementedError, "Tensorflow does not support "
+                                                 "multiplication of two sparse "
+                                                 "matrices.")
+            return self.sparse.sparse_dense_matmul(x, y)
+        return self.tensordot(x, y, axes=1)
+
     def outer(self, x, y):
         return self.tensordot(x, y, axes=0)
 
diff --git a/src/qibo/core/callbacks.py b/src/qibo/core/callbacks.py
index 5c3fbd476..c04379dbb 100644
--- a/src/qibo/core/callbacks.py
+++ b/src/qibo/core/callbacks.py
@@ -149,7 +149,7 @@ class Gap(BackendCallback, abstract_callbacks.Gap):
         if self.evolution is None:
             raise_error(ValueError, "Gap callback can only be used in "
                                     "adiabatic evolution models.")
-        hamiltonian = self.evolution.solver.current_hamiltonian
+        hamiltonian = self.evolution.solver.current_hamiltonian  # pylint: disable=E1101
         # Call the eigenvectors so that they are cached for the ``exp`` call
         hamiltonian.eigenvectors()
         eigvals = hamiltonian.eigenvalues()
diff --git a/src/qibo/core/hamiltonians.py b/src/qibo/core/hamiltonians.py
index 5c5eaa81c..76c629643 100644
--- a/src/qibo/core/hamiltonians.py
+++ b/src/qibo/core/hamiltonians.py
@@ -6,16 +6,19 @@ from qibo.abstractions import hamiltonians, states
 
 
 class Hamiltonian(hamiltonians.MatrixHamiltonian):
-    """Backend implementation of :class:`qibo.abstractions.hamiltonians.MatrixHamiltonian`.
+    """Hamiltonian based on a dense or sparse matrix representation.
 
     Args:
         nqubits (int): number of quantum bits.
         matrix (np.ndarray): Matrix representation of the Hamiltonian in the
             computational basis as an array of shape ``(2 ** nqubits, 2 ** nqubits)``.
+            Sparse matrices based on ``scipy.sparse`` for numpy/qibojit backends
+            or on ``tf.sparse`` for tensorflow/qibotf backends are also
+            supported.
     """
 
     def __init__(self, nqubits, matrix):
-        if not isinstance(matrix, K.tensor_types):
+        if not (isinstance(matrix, K.tensor_types) or K.issparse(matrix)):
             raise_error(TypeError, "Matrix of invalid type {} given during "
                                    "Hamiltonian initialization"
                                    "".format(type(matrix)))
@@ -46,20 +49,20 @@ class Hamiltonian(hamiltonians.MatrixHamiltonian):
                     "to construct Hamiltonians using symbols.")
         return SymbolicHamiltonian(symbolic_hamiltonian, symbol_map)
 
-    def eigenvalues(self):
+    def eigenvalues(self, k=6):
         if self._eigenvalues is None:
-            self._eigenvalues = self.K.eigvalsh(self.matrix)
+            self._eigenvalues = self.K.eigvalsh(self.matrix, k)
         return self._eigenvalues
 
-    def eigenvectors(self):
+    def eigenvectors(self, k=6):
         if self._eigenvectors is None:
-            self._eigenvalues, self._eigenvectors = self.K.eigh(self.matrix)
+            self._eigenvalues, self._eigenvectors = self.K.eigh(self.matrix, k)
         return self._eigenvectors
 
     def exp(self, a):
         if self._exp.get("a") != a:
             self._exp["a"] = a
-            if self._eigenvectors is None:
+            if self._eigenvectors is None or self.K.issparse(self.matrix):
                 self._exp["result"] = self.K.expm(-1j * a * self.matrix)
             else:
                 expd = self.K.diag(self.K.exp(-1j * a * self._eigenvalues))
@@ -146,7 +149,7 @@ class Hamiltonian(hamiltonians.MatrixHamiltonian):
         if self._eigenvalues is not None:
             if K.qnp.cast(o).real >= 0:
                 r._eigenvalues = o * self._eigenvalues
-            else:
+            elif not K.issparse(self.matrix):
                 r._eigenvalues = o * self._eigenvalues[::-1]
         if self._eigenvectors is not None:
             if K.qnp.cast(o).real > 0:
@@ -157,7 +160,7 @@ class Hamiltonian(hamiltonians.MatrixHamiltonian):
 
     def __matmul__(self, o):
         if isinstance(o, self.__class__):
-            new_matrix = self.K.matmul(self.matrix, o.matrix)
+            new_matrix = self.K.dot(self.matrix, o.matrix)
             return self.__class__(self.nqubits, new_matrix)
 
         if isinstance(o, states.AbstractState):
@@ -165,9 +168,9 @@ class Hamiltonian(hamiltonians.MatrixHamiltonian):
         if isinstance(o, K.tensor_types):
             rank = len(tuple(o.shape))
             if rank == 1: # vector
-                return self.K.matmul(self.matrix, o[:, self.K.newaxis])[:, 0]
+                return self.K.dot(self.matrix, o[:, self.K.newaxis])[:, 0]
             elif rank == 2: # matrix
-                return self.K.matmul(self.matrix, o)
+                return self.K.dot(self.matrix, o)
             else:
                 raise_error(ValueError, "Cannot multiply Hamiltonian with "
                                         "rank-{} tensor.".format(rank))


Test patch:
diff --git a/src/qibo/tests/test_backends_agreement.py b/src/qibo/tests/test_backends_agreement.py
index 19870e69b..2cc87fa90 100644
--- a/src/qibo/tests/test_backends_agreement.py
+++ b/src/qibo/tests/test_backends_agreement.py
@@ -101,14 +101,46 @@ def test_backend_stack(tested_backend, target_backend):
     tested_backend.assert_allclose(tested_result, target_result)
 
 
-def test_backend_eigh(tested_backend, target_backend):
+@pytest.mark.parametrize("sparse_type", [None, "coo", "csr", "csc", "dia"])
+def test_backend_eigh(tested_backend, target_backend, sparse_type):
+    if tested_backend in ("qibotf", "tensorflow") and sparse_type is not None:
+        pytest.skip("Temporary skip.")
     tested_backend = K.construct_backend(tested_backend)
     target_backend = K.construct_backend(target_backend)
-    m = rand((5, 5))
-    eigvals2, eigvecs2 = target_backend.eigh(m)
+    if sparse_type is None:
+        m = rand((5, 5))
+    else:
+        from scipy import sparse
+        m = sparse.rand(16, 16, format=sparse_type, density=0.1)
+    eigvals1, eigvecs1 = tested_backend.eigh(tested_backend.cast(m), k=m.shape[0])
+    eigvals2, eigvecs2 = target_backend.eigh(target_backend.cast(m), k=m.shape[0])
+    tested_backend.assert_allclose(eigvals1, eigvals2, atol=1e-10)
+    m1 = eigvecs1 @ np.diag(eigvals1) @ tested_backend.transpose(eigvecs1)
+    m2 = eigvecs2 @ np.diag(eigvals2) @ target_backend.transpose(eigvecs2)
+    tested_backend.assert_allclose(m1, m2, atol=1e-10)
+
+
+@pytest.mark.parametrize("sparse_type", ["coo", "csr", "csc", "dia"])
+def test_backend_sparse_eigh(tested_backend, target_backend, sparse_type):
+    if tested_backend in ("qibotf", "tensorflow"):
+        pytest.skip("Temporary skip.")
+    tested_backend = K.construct_backend(tested_backend)
+    target_backend = K.construct_backend(target_backend)
+    from scipy import sparse
+    from qibo import hamiltonians
+    ham = hamiltonians.TFIM(6, h=1.0)
+    m = getattr(sparse, f"{sparse_type}_matrix")(K.to_numpy(ham.matrix))
     eigvals1, eigvecs1 = tested_backend.eigh(tested_backend.cast(m))
+    eigvals2, eigvecs2 = target_backend.eigh(target_backend.cast(m))
+    eigvals1 = sorted(K.to_numpy(eigvals1))
+    eigvals2 = sorted(K.to_numpy(eigvals2))
+    tested_backend.assert_allclose(eigvals1, eigvals2)
+
+    eigvals1 = tested_backend.eigvalsh(tested_backend.cast(m))
+    eigvals2 = target_backend.eigvalsh(target_backend.cast(m))
+    eigvals1 = sorted(K.to_numpy(eigvals1))
+    eigvals2 = sorted(K.to_numpy(eigvals2))
     tested_backend.assert_allclose(eigvals1, eigvals2)
-    tested_backend.assert_allclose(np.abs(eigvecs1), np.abs(eigvecs2))
 
 
 def test_backend_compile(tested_backend, target_backend):
diff --git a/src/qibo/tests/test_core_hamiltonians.py b/src/qibo/tests/test_core_hamiltonians.py
index 89825d4d7..393450fc2 100644
--- a/src/qibo/tests/test_core_hamiltonians.py
+++ b/src/qibo/tests/test_core_hamiltonians.py
@@ -1,10 +1,24 @@
 """Test methods in `qibo/core/hamiltonians.py`."""
 import pytest
 import numpy as np
+from scipy import sparse
 from qibo import hamiltonians, K
 from qibo.tests.utils import random_complex
 
 
+def random_sparse_matrix(n, sparse_type=None):
+    if K.name in ("qibotf", "tensorflow"):
+        nonzero = int(0.1 * n * n)
+        indices = np.random.randint(0, n, size=(nonzero, 2))
+        data = np.random.random(nonzero) + 1j * np.random.random(nonzero)
+        data = K.cast(data)
+        return K.sparse.SparseTensor(indices, data, (n, n))
+    else:
+        re = sparse.rand(n, n, format=sparse_type)
+        im = sparse.rand(n, n, format=sparse_type)
+        return re + 1j * im
+
+
 def test_hamiltonian_init():
     with pytest.raises(TypeError):
         H = hamiltonians.Hamiltonian(2, "test")
@@ -21,7 +35,8 @@ def test_hamiltonian_init():
 
 
 @pytest.mark.parametrize("dtype", K.numeric_types)
-def test_hamiltonian_algebraic_operations(dtype):
+@pytest.mark.parametrize("sparse_type", [None, "coo", "csr", "csc", "dia"])
+def test_hamiltonian_algebraic_operations(dtype, sparse_type):
     """Test basic hamiltonian overloading."""
 
     def transformation_a(a, b):
@@ -45,10 +60,16 @@ def test_hamiltonian_algebraic_operations(dtype):
         else:
             return c1 - a + c2 * b
 
-    H1 = hamiltonians.XXZ(nqubits=2, delta=0.5)
-    H2 = hamiltonians.XXZ(nqubits=2, delta=1)
+    if sparse_type is None:
+        H1 = hamiltonians.XXZ(nqubits=2, delta=0.5)
+        H2 = hamiltonians.XXZ(nqubits=2, delta=1)
+        mH1, mH2 = K.to_numpy(H1.matrix), K.to_numpy(H2.matrix)
+    else:
+        mH1 = sparse.rand(64, 64, format=sparse_type)
+        mH2 = sparse.rand(64, 64, format=sparse_type)
+        H1 = hamiltonians.Hamiltonian(6, mH1)
+        H2 = hamiltonians.Hamiltonian(6, mH2)
 
-    mH1, mH2 = K.to_numpy(H1.matrix), K.to_numpy(H2.matrix)
     hH1 = transformation_a(mH1, mH2)
     hH2 = transformation_b(mH1, mH2)
     hH3 = transformation_c(mH1, mH2, use_eye=True)
@@ -65,9 +86,15 @@ def test_hamiltonian_algebraic_operations(dtype):
     K.assert_allclose(hH4, HT4.matrix)
 
 
-def test_hamiltonian_addition():
-    H1 = hamiltonians.Y(nqubits=3)
-    H2 = hamiltonians.TFIM(nqubits=3, h=1.0)
+@pytest.mark.parametrize("sparse_type", [None, "coo", "csr", "csc", "dia"])
+def test_hamiltonian_addition(sparse_type):
+    if sparse_type is None:
+        H1 = hamiltonians.Y(nqubits=3)
+        H2 = hamiltonians.TFIM(nqubits=3, h=1.0)
+    else:
+        H1 = hamiltonians.Hamiltonian(6, sparse.rand(64, 64, format=sparse_type))
+        H2 = hamiltonians.Hamiltonian(6, sparse.rand(64, 64, format=sparse_type))
+
     H = H1 + H2
     matrix = H1.matrix + H2.matrix
     K.assert_allclose(H.matrix, matrix)
@@ -98,58 +125,78 @@ def test_hamiltonian_operation_errors():
         R = [3] - H1
 
 
-def test_hamiltonian_matmul():
-    """Test matrix multiplication between Hamiltonians and state vectors."""
-    H1 = hamiltonians.TFIM(nqubits=3, h=1.0)
-    H2 = hamiltonians.Y(nqubits=3)
+@pytest.mark.parametrize("sparse_type", [None, "coo", "csr", "csc", "dia"])
+def test_hamiltonian_matmul(backend, sparse_type):
+    """Test matrix multiplication between Hamiltonians."""
+    if sparse_type is None:
+        nqubits = 3
+        H1 = hamiltonians.TFIM(nqubits, h=1.0)
+        H2 = hamiltonians.Y(nqubits)
+    else:
+        nqubits = 5
+        nstates = 2 ** nqubits
+        H1 = hamiltonians.Hamiltonian(nqubits, random_sparse_matrix(nstates, sparse_type))
+        H2 = hamiltonians.Hamiltonian(nqubits, random_sparse_matrix(nstates, sparse_type))
+
     m1 = K.to_numpy(H1.matrix)
     m2 = K.to_numpy(H2.matrix)
-
-    K.assert_allclose((H1 @ H2).matrix, m1 @ m2)
-    K.assert_allclose((H2 @ H1).matrix, m2 @ m1)
-
-    v = random_complex(8, dtype=m1.dtype)
-    m = random_complex((8, 8), dtype=m1.dtype)
-    H1v = H1 @ K.cast(v)
-    H1m = H1 @ K.cast(m)
-    K.assert_allclose(H1v, m1.dot(v))
-    K.assert_allclose(H1m, m1 @ m)
-
-    from qibo.core.states import VectorState
-    H1state = H1 @ VectorState.from_tensor(K.cast(v))
-    K.assert_allclose(H1state, m1.dot(v))
+    if K.name in ("qibotf", "tensorflow") and sparse_type is not None:
+        with pytest.raises(NotImplementedError):
+            _ = H1 @ H2
+    else:
+        K.assert_allclose((H1 @ H2).matrix, m1 @ m2)
+        K.assert_allclose((H2 @ H1).matrix, m2 @ m1)
 
     with pytest.raises(ValueError):
-        H1 @ np.zeros((8, 8, 8), dtype=m1.dtype)
+        H1 @ np.zeros(3 * (2 ** nqubits,), dtype=m1.dtype)
     with pytest.raises(NotImplementedError):
         H1 @ 2
 
 
-@pytest.mark.parametrize("dense", [True, False])
-def test_hamiltonian_exponentiation(dense):
-    from scipy.linalg import expm
-    H = hamiltonians.XXZ(nqubits=2, delta=0.5, dense=dense)
-    target_matrix = expm(-0.5j * K.to_numpy(H.matrix))
-    K.assert_allclose(H.exp(0.5), target_matrix)
+@pytest.mark.parametrize("sparse_type", [None, "coo", "csr", "csc", "dia"])
+def test_hamiltonian_matmul_states(backend, sparse_type):
+    """Test matrix multiplication between Hamiltonian and states."""
+    if sparse_type is None:
+        nqubits = 3
+        H = hamiltonians.TFIM(nqubits, h=1.0)
+    else:
+        nqubits = 5
+        nstates = 2 ** nqubits
+        H = hamiltonians.Hamiltonian(nqubits, random_sparse_matrix(nstates, sparse_type))
+
+    hm = K.to_numpy(H.matrix)
+    v = random_complex(2 ** nqubits, dtype=hm.dtype)
+    m = random_complex((2 ** nqubits, 2 ** nqubits), dtype=hm.dtype)
+    Hv = H @ K.cast(v)
+    Hm = H @ K.cast(m)
+    K.assert_allclose(Hv, hm.dot(v))
+    K.assert_allclose(Hm, hm @ m)
 
-    H = hamiltonians.XXZ(nqubits=2, delta=0.5, dense=dense)
-    _ = H.eigenvectors()
-    K.assert_allclose(H.exp(0.5), target_matrix)
+    from qibo.core.states import VectorState
+    Hstate = H @ VectorState.from_tensor(K.cast(v))
+    K.assert_allclose(Hstate, hm.dot(v))
 
 
-@pytest.mark.parametrize("dense", [True, False])
 @pytest.mark.parametrize("density_matrix", [True, False])
-def test_hamiltonian_expectation(dense, density_matrix):
-    h = hamiltonians.XXZ(nqubits=3, delta=0.5, dense=dense)
-    matrix = K.to_numpy(h.matrix)
+@pytest.mark.parametrize("sparse_type,dense",
+                         [(None, True), (None, False),
+                          ("coo", True), ("csr", True),
+                          ("csc", True), ("dia", True)])
+def test_hamiltonian_expectation(backend, dense, density_matrix, sparse_type):
+    """Test Hamiltonian expectation value calculation."""
+    if sparse_type is None:
+        h = hamiltonians.XXZ(nqubits=3, delta=0.5, dense=dense)
+    else:
+        h = hamiltonians.Hamiltonian(6, random_sparse_matrix(64, sparse_type))
 
+    matrix = K.to_numpy(h.matrix)
     if density_matrix:
-        state = random_complex((8, 8))
+        state = random_complex((2 ** h.nqubits, 2 ** h.nqubits))
         state = state + state.T.conj()
         norm = np.trace(state)
         target_ev = np.trace(matrix.dot(state)).real
     else:
-        state = random_complex(8)
+        state = random_complex(2 ** h.nqubits)
         norm = np.sum(np.abs(state) ** 2)
         target_ev = np.sum(state.conj() * matrix.dot(state)).real
 
@@ -167,24 +214,38 @@ def test_hamiltonian_expectation_errors():
 
 
 @pytest.mark.parametrize("dtype", K.numeric_types)
-@pytest.mark.parametrize("dense", [True, False])
-def test_hamiltonian_eigenvalues(dtype, dense):
+@pytest.mark.parametrize("sparse_type,dense",
+                         [(None, True), (None, False),
+                          ("coo", True), ("csr", True),
+                          ("csc", True), ("dia", True)])
+def test_hamiltonian_eigenvalues(dtype, sparse_type, dense):
     """Testing hamiltonian eigenvalues scaling."""
-    H1 = hamiltonians.XXZ(nqubits=2, delta=0.5, dense=dense)
+    if sparse_type is None:
+        H1 = hamiltonians.XXZ(nqubits=2, delta=0.5, dense=dense)
+    else:
+        from scipy import sparse
+        H1 = hamiltonians.XXZ(nqubits=5, delta=0.5)
+        m = getattr(sparse, f"{sparse_type}_matrix")(K.to_numpy(H1.matrix))
+        H1 = hamiltonians.Hamiltonian(5, m)
 
-    H1_eigen = H1.eigenvalues()
-    hH1_eigen = K.eigvalsh(H1.matrix)
-    K.assert_allclose(H1_eigen, hH1_eigen)
+    H1_eigen = sorted(K.to_numpy(H1.eigenvalues()))
+    hH1_eigen = sorted(K.to_numpy(K.eigvalsh(H1.matrix)))
+    K.assert_allclose(sorted(H1_eigen), hH1_eigen)
 
     c1 = dtype(2.5)
     H2 = c1 * H1
-    hH2_eigen = K.eigvalsh(c1 * H1.matrix)
-    K.assert_allclose(H2._eigenvalues, hH2_eigen)
+    H2_eigen = sorted(K.to_numpy(H2._eigenvalues))
+    hH2_eigen = sorted(K.to_numpy(K.eigvalsh(c1 * H1.matrix)))
+    K.assert_allclose(H2_eigen, hH2_eigen)
 
     c2 = dtype(-11.1)
     H3 = H1 * c2
-    hH3_eigen = K.eigvalsh(H1.matrix * c2)
-    K.assert_allclose(H3._eigenvalues, hH3_eigen)
+    if sparse_type is None:
+        H3_eigen = sorted(K.to_numpy(H3._eigenvalues))
+        hH3_eigen = sorted(K.to_numpy(K.eigvalsh(H1.matrix * c2)))
+        K.assert_allclose(H3_eigen, hH3_eigen)
+    else:
+        assert H3._eigenvalues is None
 
 
 @pytest.mark.parametrize("dtype", K.numeric_types)
@@ -196,8 +257,6 @@ def test_hamiltonian_eigenvectors(dtype, dense):
     V1 = K.to_numpy(H1.eigenvectors())
     U1 = K.to_numpy(H1.eigenvalues())
     K.assert_allclose(H1.matrix, V1 @ np.diag(U1) @ V1.T)
-    # Check ground state
-    K.assert_allclose(H1.ground_state(), V1[:, 0])
 
     c1 = dtype(2.5)
     H2 = c1 * H1
@@ -216,3 +275,44 @@ def test_hamiltonian_eigenvectors(dtype, dense):
     V4 = K.to_numpy(H4._eigenvectors)
     U4 = K.to_numpy(H4._eigenvalues)
     K.assert_allclose(H4.matrix, V4 @ np.diag(U4) @ V4.T)
+
+
+@pytest.mark.parametrize("sparse_type,dense",
+                         [(None, True), (None, False),
+                          ("coo", True), ("csr", True),
+                          ("csc", True), ("dia", True)])
+def test_hamiltonian_ground_state(sparse_type, dense):
+    """Test Hamiltonian ground state."""
+    if sparse_type is None:
+        H = hamiltonians.XXZ(nqubits=2, delta=0.5, dense=dense)
+    else:
+        from scipy import sparse
+        H = hamiltonians.XXZ(nqubits=5, delta=0.5)
+        m = getattr(sparse, f"{sparse_type}_matrix")(K.to_numpy(H.matrix))
+        H = hamiltonians.Hamiltonian(5, m)
+    V = K.to_numpy(H.eigenvectors())
+    K.assert_allclose(H.ground_state(), V[:, 0])
+
+
+@pytest.mark.parametrize("sparse_type,dense",
+                         [(None, True), (None, False),
+                          ("coo", True), ("csr", True),
+                          ("csc", True), ("dia", True)])
+def test_hamiltonian_exponentiation(sparse_type, dense):
+    """Test matrix exponentiation of Hamiltonians ``exp(1j * t * H)``."""
+    from scipy.linalg import expm
+    def construct_hamiltonian():
+        if sparse_type is None:
+            return hamiltonians.XXZ(nqubits=2, delta=0.5, dense=dense)
+        else:
+            ham = hamiltonians.XXZ(nqubits=5, delta=0.5)
+            m = getattr(sparse, f"{sparse_type}_matrix")(K.to_numpy(ham.matrix))
+            return hamiltonians.Hamiltonian(5, m)
+
+    H = construct_hamiltonian()
+    target_matrix = expm(-0.5j * K.to_numpy(H.matrix))
+    K.assert_allclose(H.exp(0.5), target_matrix)
+
+    H = construct_hamiltonian()
+    _ = H.eigenvectors()
+    K.assert_allclose(H.exp(0.5), target_matrix)