Repository: qiboteam/qibo

Issue:
sympy vs qibo symbols performance
We should verify the impact of the `expand` call in qibo.
It looks like this call deteriorates our performance in comparison to a pure sympy implementation.

Discussion:

Patch:
diff --git a/doc/source/advancedexamples.rst b/doc/source/advancedexamples.rst
index 994e224b1..1374722c6 100644
--- a/doc/source/advancedexamples.rst
+++ b/doc/source/advancedexamples.rst
@@ -1427,3 +1427,20 @@ form is very close to the form of the Hamiltonian on paper. Note that when a
 :class:`qibo.core.hamiltonians.SymbolicHamiltonian` is used for time evolution,
 Qibo handles automatically automatically the Trotter decomposition by splitting
 to the appropriate terms.
+
+Qibo symbols support an additional ``commutative`` argument which is set to
+``False`` by default since quantum operators are non-commuting objects.
+When the user knows that the Hamiltonian consists of commuting terms only, such
+as products of Z operators, switching ``commutative`` to ``True`` may speed-up
+some symbolic calculations, such as the ``sympy.expand`` used when calculating
+the Trotter decomposition for the Hamiltonian. This option can be used when
+constructing each symbol:
+
+
+.. code-block::  python
+
+    from qibo import hamiltonians
+    from qibo.symbols import Z
+
+    form = Z(0, commutative=True) * Z(1, commutative=True) + Z(1, commutative=True) * Z(2, commutative=True)
+    ham = hamiltonians.SymbolicHamiltonian(form)
diff --git a/src/qibo/core/hamiltonians.py b/src/qibo/core/hamiltonians.py
index f3a4056ec..5c5eaa81c 100644
--- a/src/qibo/core/hamiltonians.py
+++ b/src/qibo/core/hamiltonians.py
@@ -1,4 +1,5 @@
 import itertools
+import sympy
 from qibo import K, gates
 from qibo.config import log, raise_error, EINSUM_CHARS
 from qibo.abstractions import hamiltonians, states
@@ -212,56 +213,187 @@ class TrotterCircuit:
 class SymbolicHamiltonian(hamiltonians.SymbolicHamiltonian):
     """Backend implementation of :class:`qibo.abstractions.hamiltonians.SymbolicHamiltonian`.
 
+    Calculations using symbolic Hamiltonians are either done directly using
+    the given ``sympy`` expression as it is (``form``) or by parsing the
+    corresponding ``terms`` (which are :class:`qibo.core.terms.SymbolicTerm`
+    objects). The latter approach is more computationally costly as it uses
+    a ``sympy.expand`` call on the given form before parsing the terms.
+    For this reason the ``terms`` are calculated only when needed, for example
+    during Trotterization.
+    The dense matrix of the symbolic Hamiltonian can be calculated directly
+    from ``form`` without requiring ``terms`` calculation (see
+    :meth:`qibo.core.hamiltonians.SymbolicHamiltonian.calculate_dense` for details).
+
     Args:
-        form (sympy.Expr): Hamiltonian form as a ``sympy.Expr``. The Hamiltonian
-            should be created using Qibo symbols.
+        form (sympy.Expr): Hamiltonian form as a ``sympy.Expr``. Ideally the
+            Hamiltonian should be written using Qibo symbols.
             See :ref:`How to define custom Hamiltonians using symbols? <symbolicham-example>`
             example for more details.
-        symbol_map (dict): Dictionary that maps each ``sympy.Symbol`` to the
-            corresponding target qubit and matrix representation. This feature
-            is kept for compatibility with older versions where Qibo symbols
-            were not available and will be deprecated in the future.
+        symbol_map (dict): Dictionary that maps each ``sympy.Symbol`` to a tuple
+            of (target qubit, matrix representation). This feature is kept for
+            compatibility with older versions where Qibo symbols were not available
+            and may be deprecated in the future.
             It is not required if the Hamiltonian is constructed using Qibo symbols.
+            The symbol_map can also be used to pass non-quantum operator arguments
+            to the symbolic Hamiltonian, such as the parameters in the
+            :meth:`qibo.hamiltonians.MaxCut` Hamiltonian.
+        ground_state (Callable): Function with no arguments that returns the
+            ground state of this Hamiltonian. This is useful in cases where
+            the ground state is trivial and is used for initialization,
+            for example the easy Hamiltonian in adiabatic evolution,
+            however we would like to avoid constructing and diagonalizing the
+            full Hamiltonian matrix only to find the ground state.
     """
 
-    def __init__(self, form=None, symbol_map=None, ground_state=None):
+    def __init__(self, form=None, symbol_map={}, ground_state=None):
         super().__init__(ground_state)
+        self._form = None
+        self._terms = None
+        self.constant = 0 # used only when we perform calculations using ``_terms``
         self._dense = None
-        self.terms = None # list of :class:`qibo.core.terms.HamiltonianTerm` objects
-        self.constant = 0
-        self.form = None
+        self.symbol_map = symbol_map
+        # if a symbol in the given form is not a Qibo symbol it must be
+        # included in the ``symbol_map``
         self.trotter_circuit = None
+        from qibo.symbols import Symbol
+        self._qiboSymbol = Symbol # also used in ``self._get_symbol_matrix``
         if form is not None:
-            self.set_form(form, symbol_map)
+            self.form = form
 
-    @classmethod
-    def from_terms(cls, terms, ground_state=None):
-        """Constructs a symbolic Hamiltonian directly from a list of terms."""
-        ham = cls(ground_state=ground_state)
-        ham.set_terms(terms)
-        return ham
-
-    def set_terms(self, terms):
-        self.terms = terms
-        self.nqubits = max(q for term in self.terms for q in term.target_qubits) + 1
-
-    def set_form(self, form, symbol_map=None):
-        import sympy
-        from qibo.core.terms import SymbolicTerm
-        if not issubclass(form.__class__, sympy.Expr):
+    @property
+    def form(self):
+        return self._form
+
+    @form.setter
+    def form(self, form):
+        # Check that given form is a ``sympy`` expression
+        if not isinstance(form, sympy.Expr):
             raise_error(TypeError, "Symbolic Hamiltonian should be a ``sympy`` "
                                    "expression but is {}.".format(type(form)))
-        self.form = sympy.expand(form)
-        terms = []
-        for f, c in self.form.as_coefficients_dict().items():
-            term = SymbolicTerm(c, f, symbol_map)
-            if term.target_qubits:
-                terms.append(term)
+        # Calculate number of qubits in the system described by the given
+        # Hamiltonian formula
+        nqubits = 0
+        for symbol in form.free_symbols:
+            if isinstance(symbol, self._qiboSymbol):
+                q = symbol.target_qubit
+            elif isinstance(symbol, sympy.Expr):
+                if symbol not in self.symbol_map:
+                    raise_error(ValueError, "Symbol {} is not in symbol "
+                                            "map.".format(symbol))
+                q, matrix = self.symbol_map.get(symbol)
+                if not isinstance(matrix, K.tensor_types):
+                    # ignore symbols that do not correspond to quantum operators
+                    # for example parameters in the MaxCut Hamiltonian
+                    q = 0
+            if q > nqubits:
+                nqubits = q
+
+        self._form = form
+        self.nqubits = nqubits + 1
+
+    @property
+    def terms(self):
+        """List of :class:`qibo.core.terms.HamiltonianTerm` objects of which the Hamiltonian is a sum of."""
+        if self._terms is None:
+            # Calculate terms based on ``self.form``
+            from qibo.core.terms import SymbolicTerm
+            form = sympy.expand(self.form)
+            terms = []
+            for f, c in form.as_coefficients_dict().items():
+                term = SymbolicTerm(c, f, self.symbol_map)
+                if term.target_qubits:
+                    terms.append(term)
+                else:
+                    self.constant += term.coefficient
+            assert self.nqubits == max(q for term in terms for q in term.target_qubits) + 1
+            self._terms = terms
+        return self._terms
+
+    @terms.setter
+    def terms(self, terms):
+        self._terms = terms
+        self.nqubits = max(q for term in self._terms for q in term.target_qubits) + 1
+
+    def _get_symbol_matrix(self, term):
+        """Calculates numerical matrix corresponding to symbolic expression.
+
+        This is partly equivalent to sympy's ``.subs``, which does not work
+        in our case as it does not allow us to substitute ``sympy.Symbol``
+        with numpy arrays and there are different complication when switching
+        to ``sympy.MatrixSymbol``. Here we calculate the full numerical matrix
+        given the symbolic expression using recursion.
+        Helper method for ``_calculate_dense_from_form``.
+
+        Args:
+            term (sympy.Expr): Symbolic expression containing local operators.
+
+        Returns:
+            Numerical matrix corresponding to the given expression as a numpy
+            array of size ``(2 ** self.nqubits, 2 ** self.nqubits).
+        """
+        if isinstance(term, sympy.Add):
+            # symbolic op for addition
+            result = sum(self._get_symbol_matrix(subterm)
+                         for subterm in term.as_ordered_terms())
+
+        elif isinstance(term, sympy.Mul):
+            # symbolic op for multiplication
+            # note that we need to use matrix multiplication even though
+            # we use scalar symbols for convenience
+            factors = term.as_ordered_factors()
+            result = self._get_symbol_matrix(factors[0])
+            for subterm in factors[1:]:
+                result = result @ self._get_symbol_matrix(subterm)
+
+        elif isinstance(term, sympy.Pow):
+            # symbolic op for power
+            base, exponent = term.as_base_exp()
+            matrix = self._get_symbol_matrix(base)
+            # multiply ``base`` matrix ``exponent`` times to itself
+            result = matrix
+            for _ in range(exponent - 1):
+                result = result @ matrix
+
+        elif isinstance(term, sympy.Symbol):
+            # if the term is a ``Symbol`` then it corresponds to a quantum
+            # operator for which we can construct the full matrix directly
+            if isinstance(term, self._qiboSymbol):
+                # if we have a Qibo symbol the matrix construction is
+                # implemented in :meth:`qibo.core.terms.SymbolicTerm.full_matrix`.
+                result = term.full_matrix(self.nqubits)
             else:
-                self.constant += term.coefficient
-        self.set_terms(terms)
+                q, matrix = self.symbol_map.get(term)
+                if not isinstance(matrix, K.tensor_types):
+                    # symbols that do not correspond to quantum operators
+                    # for example parameters in the MaxCut Hamiltonian
+                    result = complex(matrix) * K.qnp.eye(2 ** self.nqubits)
+                else:
+                    # if we do not have a Qibo symbol we construct one and use
+                    # :meth:`qibo.core.terms.SymbolicTerm.full_matrix`.
+                    result = self._qiboSymbol(q, matrix).full_matrix(self.nqubits)
+
+        elif term.is_number:
+            # if the term is number we should return in the form of identity
+            # matrix because in expressions like `1 + Z`, `1` is not correspond
+            # to the float 1 but the identity operator (matrix)
+            result = complex(term) * K.qnp.eye(2 ** self.nqubits)
 
-    def calculate_dense(self):
+        else:
+            raise_error(TypeError, "Cannot calculate matrix for symbolic term "
+                                   "of type {}.".format(type(term)))
+
+        return result
+
+    def _calculate_dense_from_form(self):
+        """Calculates equivalent :class:`qibo.core.hamiltonians.Hamiltonian` using symbolic form.
+
+        Useful when the term representation is not available.
+        """
+        matrix = self._get_symbol_matrix(self.form)
+        return Hamiltonian(self.nqubits, matrix)
+
+    def _calculate_dense_from_terms(self):
+        """Calculates equivalent :class:`qibo.core.hamiltonians.Hamiltonian` using the term representation."""
         if 2 * self.nqubits > len(EINSUM_CHARS): # pragma: no cover
             # case not tested because it only happens in large examples
             raise_error(NotImplementedError, "Not enough einsum characters.")
@@ -280,6 +412,13 @@ class SymbolicHamiltonian(hamiltonians.SymbolicHamiltonian):
         matrix = K.np.reshape(matrix, 2 * (2 ** self.nqubits,))
         return Hamiltonian(self.nqubits, matrix) + self.constant
 
+    def calculate_dense(self):
+        if self._terms is None:
+            # calculate dense matrix directly using the form to avoid the
+            # costly ``sympy.expand`` call
+            return self._calculate_dense_from_form()
+        return self._calculate_dense_from_terms()
+
     def expectation(self, state, normalize=False):
         return Hamiltonian.expectation(self, state, normalize)
 
@@ -288,17 +427,23 @@ class SymbolicHamiltonian(hamiltonians.SymbolicHamiltonian):
             if self.nqubits != o.nqubits:
                 raise_error(RuntimeError, "Only hamiltonians with the same "
                                           "number of qubits can be added.")
-            new_ham = self.__class__.from_terms(self.terms + o.terms)
-            if self.form is not None and o.form is not None:
+            new_ham = self.__class__(symbol_map=dict(self.symbol_map))
+            if self._form is not None and o._form is not None:
                 new_ham.form = self.form + o.form
+                new_ham.symbol_map.update(o.symbol_map)
+            if self._terms is not None and o._terms is not None:
+                new_ham.terms = self.terms + o.terms
+                new_ham.constant = self.constant + o.constant
             if self._dense is not None and o._dense is not None:
                 new_ham.dense = self.dense + o.dense
 
         elif isinstance(o, K.numeric_types):
-            new_ham = self.__class__.from_terms(self.terms)
-            new_ham.constant = self.constant + o
-            if self.form is not None:
+            new_ham = self.__class__(symbol_map=dict(self.symbol_map))
+            if self._form is not None:
                 new_ham.form = self.form + o
+            if self._terms is not None:
+                new_ham.terms = self.terms
+                new_ham.constant = self.constant + o
             if self._dense is not None:
                 new_ham.dense = self.dense + o
 
@@ -312,18 +457,23 @@ class SymbolicHamiltonian(hamiltonians.SymbolicHamiltonian):
             if self.nqubits != o.nqubits:
                 raise_error(RuntimeError, "Only hamiltonians with the same "
                                           "number of qubits can be subtracted.")
-            new_terms = self.terms + [-1 * x for x in o.terms]
-            new_ham = self.__class__.from_terms(new_terms)
-            if self.form is not None and o.form is not None:
+            new_ham = self.__class__(symbol_map=dict(self.symbol_map))
+            if self._form is not None and o._form is not None:
                 new_ham.form = self.form - o.form
+                new_ham.symbol_map.update(o.symbol_map)
+            if self._terms is not None and o._terms is not None:
+                new_ham.terms = self.terms + [-1 * x for x in o.terms]
+                new_ham.constant = self.constant - o.constant
             if self._dense is not None and o._dense is not None:
                 new_ham.dense = self.dense - o.dense
 
         elif isinstance(o, K.numeric_types):
-            new_ham = self.__class__.from_terms(self.terms)
-            new_ham.constant = self.constant - o
-            if self.form is not None:
+            new_ham = self.__class__(symbol_map=dict(self.symbol_map))
+            if self._form is not None:
                 new_ham.form = self.form - o
+            if self._terms is not None:
+                new_ham.terms = self.terms
+                new_ham.constant = self.constant - o
             if self._dense is not None:
                 new_ham.dense = self.dense - o
 
@@ -334,10 +484,12 @@ class SymbolicHamiltonian(hamiltonians.SymbolicHamiltonian):
 
     def __rsub__(self, o):
         if isinstance(o, K.numeric_types):
-            new_ham = self.__class__.from_terms([-1 * x for x in self.terms])
-            new_ham.constant = o - self.constant
-            if self.form is not None:
+            new_ham = self.__class__(symbol_map=dict(self.symbol_map))
+            if self._form is not None:
                 new_ham.form = o - self.form
+            if self._terms is not None:
+                new_ham.terms = [-1 * x for x in self.terms]
+                new_ham.constant = o - self.constant
             if self._dense is not None:
                 new_ham.dense = o - self.dense
         else:
@@ -350,10 +502,12 @@ class SymbolicHamiltonian(hamiltonians.SymbolicHamiltonian):
             raise_error(NotImplementedError, "Hamiltonian multiplication to {} "
                                              "not implemented.".format(type(o)))
         o = complex(o)
-        new_ham = self.__class__.from_terms([o * x for x in self.terms])
-        new_ham.constant = self.constant * o
-        if self.form is not None:
+        new_ham = self.__class__(symbol_map=dict(self.symbol_map))
+        if self._form is not None:
             new_ham.form = o * self.form
+        if self._terms is not None:
+            new_ham.terms = [o * x for x in self.terms]
+            new_ham.constant = self.constant * o
         if self._dense is not None:
             new_ham.dense = o * self._dense
         return new_ham
@@ -373,11 +527,13 @@ class SymbolicHamiltonian(hamiltonians.SymbolicHamiltonian):
     def __matmul__(self, o):
         """Matrix multiplication with other Hamiltonians or state vectors."""
         if isinstance(o, self.__class__):
-            if self.form is None or o.form is None:
+            if self._form is None or o._form is None:
                 raise_error(NotImplementedError, "Multiplication of symbolic Hamiltonians "
                                                  "without symbolic form is not implemented.")
             new_form = self.form * o.form
-            new_ham = self.__class__(new_form)
+            new_symbol_map = dict(self.symbol_map)
+            new_symbol_map.update(o.symbol_map)
+            new_ham = self.__class__(new_form, symbol_map=new_symbol_map)
             if self._dense is not None and o._dense is not None:
                 new_ham.dense = self.dense @ o.dense
             return new_ham
diff --git a/src/qibo/core/terms.py b/src/qibo/core/terms.py
index 32306801b..7f4370d75 100644
--- a/src/qibo/core/terms.py
+++ b/src/qibo/core/terms.py
@@ -124,7 +124,7 @@ class SymbolicTerm(HamiltonianTerm):
             symbols were not available.
     """
 
-    def __init__(self, coefficient, factors=1, symbol_map=None):
+    def __init__(self, coefficient, factors=1, symbol_map={}):
         self.coefficient = complex(coefficient)
         self._matrix = None
         self._gate = None
@@ -148,7 +148,7 @@ class SymbolicTerm(HamiltonianTerm):
 
                 # if the user is using ``symbol_map`` instead of qibo symbols,
                 # create the corresponding symbols
-                if symbol_map is not None and factor in symbol_map:
+                if factor in symbol_map:
                     from qibo.symbols import Symbol
                     q, matrix = symbol_map.get(factor)
                     factor = Symbol(q, matrix, name=factor.name)
@@ -171,6 +171,8 @@ class SymbolicTerm(HamiltonianTerm):
                         self.coefficient *= factor.matrix
                 elif factor == sympy.I:
                     self.coefficient *= 1j
+                elif factor.is_number:
+                    self.coefficient *= complex(factor)
                 else: # pragma: no cover
                     raise_error(TypeError, "Cannot parse factor {}.".format(factor))
 
diff --git a/src/qibo/hamiltonians.py b/src/qibo/hamiltonians.py
index c4da00f01..ac56463f3 100644
--- a/src/qibo/hamiltonians.py
+++ b/src/qibo/hamiltonians.py
@@ -61,7 +61,9 @@ def XXZ(nqubits, delta=0.5, dense=True):
     matrix = hx + hy + delta * hz
     terms = [HamiltonianTerm(matrix, i, i + 1) for i in range(nqubits - 1)]
     terms.append(HamiltonianTerm(matrix, nqubits - 1, 0))
-    return SymbolicHamiltonian.from_terms(terms)
+    ham = SymbolicHamiltonian()
+    ham.terms = terms
+    return ham
 
 
 def _OneBodyPauli(nqubits, matrix, dense=True, ground_state=None):
@@ -73,7 +75,9 @@ def _OneBodyPauli(nqubits, matrix, dense=True, ground_state=None):
 
     matrix = - matrix
     terms = [HamiltonianTerm(matrix, i) for i in range(nqubits)]
-    return SymbolicHamiltonian.from_terms(terms, ground_state)
+    ham = SymbolicHamiltonian(ground_state=ground_state)
+    ham.terms = terms
+    return ham
 
 
 def X(nqubits, dense=True):
@@ -150,7 +154,9 @@ def TFIM(nqubits, h=0.0, dense=True):
     matrix = -(K.np.kron(matrices.Z, matrices.Z) + h * K.np.kron(matrices.X, matrices.I))
     terms = [HamiltonianTerm(matrix, i, i + 1) for i in range(nqubits - 1)]
     terms.append(HamiltonianTerm(matrix, nqubits - 1, 0))
-    return SymbolicHamiltonian.from_terms(terms)
+    ham = SymbolicHamiltonian()
+    ham.terms = terms
+    return ham
 
 
 def MaxCut(nqubits, dense=True):
diff --git a/src/qibo/symbols.py b/src/qibo/symbols.py
index 2ecdc6658..7e9eeacc3 100644
--- a/src/qibo/symbols.py
+++ b/src/qibo/symbols.py
@@ -7,19 +7,38 @@ class Symbol(sympy.Symbol):
     """Qibo specialization for ``sympy`` symbols.
 
     These symbols can be used to create :class:`qibo.core.hamiltonians.SymbolicHamiltonian`.
+    See :ref:`How to define custom Hamiltonians using symbols? <symbolicham-example>`
+    for more details.
+
+    Example:
+        ::
+
+            from qibo import hamiltonians
+            from qibo.symbols import X, Y, Z
+            # construct a XYZ Hamiltonian on two qubits using Qibo symbols
+            form = X(0) * X(1) + Y(0) * Y(1) + Z(0) * Z(1)
+            ham = hamiltonians.SymbolicHamiltonian(form)
 
     Args:
         q (int): Target qubit id.
         matrix (np.ndarray): 2x2 matrix represented by this symbol.
         name (str): Name of the symbol which defines how it is represented in
             symbolic expressions.
+        commutative (bool): If ``True`` the constructed symbols commute with
+            each other. Default is ``False``.
+            This argument should be used with caution because quantum operators
+            are not commutative objects and therefore switching this to ``True``
+            may lead to wrong results. It is useful for improving performance
+            in symbolic calculations in cases where the user is sure that
+            the operators participating in the Hamiltonian form are commuting
+            (for example when the Hamiltonian consists of Z terms only).
     """
 
-    def __new__(cls, q, matrix=None, name="Symbol"):
+    def __new__(cls, q, matrix=None, name="Symbol", commutative=False):
         name = "{}{}".format(name, q)
-        return super().__new__(cls=cls, name=name, commutative=False)
+        return super().__new__(cls=cls, name=name, commutative=commutative)
 
-    def __init__(self, q, matrix=None, name="Symbol"):
+    def __init__(self, q, matrix=None, name="Symbol", commutative=False):
         self.target_qubit = q
         self._gate = None
         if not (matrix is None or isinstance(matrix, K.qnp.numeric_types) or
@@ -38,16 +57,34 @@ class Symbol(sympy.Symbol):
     def calculate_gate(self):
         return gates.Unitary(self.matrix, self.target_qubit)
 
+    def full_matrix(self, nqubits):
+        """Calculates the full dense matrix corresponding to the symbol as part of a bigger system.
+
+        Args:
+            nqubits (int): Total number of qubits in the system.
+
+        Returns:
+            Matrix of dimension (2^nqubits, 2^nqubits) composed of the Kronecker
+            product between identities and the symbol's single-qubit matrix.
+        """
+        from qibo.hamiltonians import multikron
+        matrix_list = self.target_qubit * [matrices.I]
+        matrix_list.append(self.matrix)
+        n = nqubits - self.target_qubit - 1
+        matrix_list.extend(matrices.I for _ in range(n))
+        return multikron(matrix_list)
+
 
 class PauliSymbol(Symbol):
 
-    def __new__(cls, q):
-        return super().__new__(cls=cls, q=q, name=cls.__name__)
+    def __new__(cls, q, commutative=False):
+        matrix = getattr(matrices, cls.__name__)
+        return super().__new__(cls, q, matrix, cls.__name__, commutative)
 
-    def __init__(self, q):
-        self.target_qubit = q
-        self._gate = None
-        self.matrix = getattr(matrices, self.__class__.__name__)
+    def __init__(self, q, commutative=False):
+        name = self.__class__.__name__
+        matrix = getattr(matrices, name)
+        super().__init__(q, matrix, name, commutative)
 
     def calculate_gate(self):
         name = self.__class__.__name__


Test patch:
diff --git a/src/qibo/tests/test_core_hamiltonians_from_symbols.py b/src/qibo/tests/test_core_hamiltonians_from_symbols.py
index 973e9bc9c..06e076cc9 100644
--- a/src/qibo/tests/test_core_hamiltonians_from_symbols.py
+++ b/src/qibo/tests/test_core_hamiltonians_from_symbols.py
@@ -8,7 +8,8 @@ from qibo.tests.utils import random_hermitian
 
 @pytest.mark.parametrize("nqubits", [4, 5])
 @pytest.mark.parametrize("hamtype", ["normal", "symbolic"])
-def test_tfim_hamiltonian_from_symbols(nqubits, hamtype):
+@pytest.mark.parametrize("calcterms", [False, True])
+def test_tfim_hamiltonian_from_symbols(nqubits, hamtype, calcterms):
     """Check creating TFIM Hamiltonian using sympy."""
     if hamtype == "symbolic":
         from qibo.symbols import X, Z
@@ -29,13 +30,16 @@ def test_tfim_hamiltonian_from_symbols(nqubits, hamtype):
         symmap.update({x: (i, matrices.X) for i, x in enumerate(x_symbols)})
         ham = hamiltonians.Hamiltonian.from_symbolic(-symham, symmap)
 
+    if calcterms:
+        _ = ham.terms
     final_matrix = ham.matrix
     target_matrix = hamiltonians.TFIM(nqubits, h=h).matrix
     K.assert_allclose(final_matrix, target_matrix)
 
 
 @pytest.mark.parametrize("hamtype", ["normal", "symbolic"])
-def test_from_symbolic_with_power(hamtype):
+@pytest.mark.parametrize("calcterms", [False, True])
+def test_from_symbolic_with_power(hamtype, calcterms):
     """Check ``from_symbolic`` when the expression contains powers."""
     if hamtype == "symbolic":
         from qibo.symbols import Symbol
@@ -50,6 +54,8 @@ def test_from_symbolic_with_power(hamtype):
         symmap = {x: (i, matrix) for i, x in enumerate(z)}
         ham = hamiltonians.Hamiltonian.from_symbolic(symham, symmap)
 
+    if calcterms:
+        _ = ham.terms
     final_matrix = ham.matrix
     matrix2 = matrix.dot(matrix)
     eye = np.eye(2, dtype=matrix.dtype)
@@ -62,7 +68,8 @@ def test_from_symbolic_with_power(hamtype):
 
 
 @pytest.mark.parametrize("hamtype", ["normal", "symbolic"])
-def test_from_symbolic_with_complex_numbers(hamtype):
+@pytest.mark.parametrize("calcterms", [False, True])
+def test_from_symbolic_with_complex_numbers(hamtype, calcterms):
     """Check ``from_symbolic`` when the expression contains imaginary unit."""
     if hamtype == "symbolic":
         from qibo.symbols import X, Y
@@ -76,6 +83,8 @@ def test_from_symbolic_with_complex_numbers(hamtype):
         symmap.update({s: (i, matrices.Y) for i, s in enumerate(y)})
         ham = hamiltonians.Hamiltonian.from_symbolic(symham, symmap)
 
+    if calcterms:
+        _ = ham.terms
     final_matrix = ham.matrix
     target_matrix = (1 + 2j) * np.kron(matrices.X, matrices.X)
     target_matrix += 2 * np.kron(matrices.Y, matrices.Y)
@@ -84,7 +93,8 @@ def test_from_symbolic_with_complex_numbers(hamtype):
     K.assert_allclose(final_matrix, target_matrix)
 
 
-def test_from_symbolic_application_hamiltonian():
+@pytest.mark.parametrize("calcterms", [False, True])
+def test_from_symbolic_application_hamiltonian(calcterms):
     """Check ``from_symbolic`` for a specific four-qubit Hamiltonian."""
     z1, z2, z3, z4 = sympy.symbols("z1 z2 z3 z4")
     symmap = {z: (i, matrices.Z) for i, z in enumerate([z1, z2, z3, z4])}
@@ -96,12 +106,15 @@ def test_from_symbolic_application_hamiltonian():
     symham = (Z(0) * Z(1) - 0.5 * Z(0) * Z(2) + 2 * Z(1) * Z(2) + 0.35 * Z(1)
               + 0.25 * Z(2) * Z(3) + 0.5 * Z(2) + Z(3) - Z(0))
     sham = hamiltonians.SymbolicHamiltonian(symham)
+    if calcterms:
+        _ = sham.terms
     K.assert_allclose(sham.matrix, fham.matrix)
 
 
 @pytest.mark.parametrize("nqubits", [4, 5])
 @pytest.mark.parametrize("hamtype", ["normal", "symbolic"])
-def test_x_hamiltonian_from_symbols(nqubits, hamtype):
+@pytest.mark.parametrize("calcterms", [False, True])
+def test_x_hamiltonian_from_symbols(nqubits, hamtype, calcterms):
     """Check creating sum(X) Hamiltonian using sympy."""
     if hamtype == "symbolic":
         from qibo.symbols import X
@@ -112,13 +125,16 @@ def test_x_hamiltonian_from_symbols(nqubits, hamtype):
         symham =  -sum(x_symbols)
         symmap = {x: (i, matrices.X) for i, x in enumerate(x_symbols)}
         ham = hamiltonians.Hamiltonian.from_symbolic(symham, symmap)
+    if calcterms:
+        _ = ham.terms
     final_matrix = ham.matrix
     target_matrix = hamiltonians.X(nqubits).matrix
     K.assert_allclose(final_matrix, target_matrix)
 
 
 @pytest.mark.parametrize("hamtype", ["normal", "symbolic"])
-def test_three_qubit_term_hamiltonian_from_symbols(hamtype):
+@pytest.mark.parametrize("calcterms", [False, True])
+def test_three_qubit_term_hamiltonian_from_symbols(hamtype, calcterms):
     """Check creating Hamiltonian with three-qubit interaction using sympy."""
     if hamtype == "symbolic":
         from qibo.symbols import X, Y, Z
@@ -142,8 +158,9 @@ def test_three_qubit_term_hamiltonian_from_symbols(hamtype):
         symham -= 2
         ham = hamiltonians.Hamiltonian.from_symbolic(symham, symmap)
 
+    if calcterms:
+        _ = ham.terms
     final_matrix = ham.matrix
-
     target_matrix = np.kron(np.kron(matrices.X, matrices.Y),
                             np.kron(matrices.Z, matrices.I))
     target_matrix += 0.5 * np.kron(np.kron(matrices.Y, matrices.Z),
diff --git a/src/qibo/tests/test_core_hamiltonians_symbolic.py b/src/qibo/tests/test_core_hamiltonians_symbolic.py
index 56304d4de..c136bc155 100644
--- a/src/qibo/tests/test_core_hamiltonians_symbolic.py
+++ b/src/qibo/tests/test_core_hamiltonians_symbolic.py
@@ -16,40 +16,39 @@ def symbolic_tfim(nqubits, h=1.0):
     return sham
 
 
-def test_symbolic_hamiltonian_init():
-    # Wrong type of symbolic expression
-    with pytest.raises(TypeError):
-        ham = hamiltonians.SymbolicHamiltonian("test")
+def test_symbolic_hamiltonian_errors():
     # Wrong type of Symbol matrix
     from qibo.symbols import Symbol
     with pytest.raises(TypeError):
         s = Symbol(0, "test")
-    # Wrong HamiltonianTerm matrix
-    from qibo.core.terms import HamiltonianTerm
+    # Wrong type of symbolic expression
     with pytest.raises(TypeError):
-        t = HamiltonianTerm("test", 0, 1)
-    # Passing negative target qubits in HamiltonianTerm
-    with pytest.raises(ValueError):
-        t = HamiltonianTerm("test", 0, -1)
-    # Passing matrix shape incompatible with number of qubits
-    with pytest.raises(ValueError):
-        t = HamiltonianTerm(np.random.random((4, 4)), 0, 1, 2)
-    # Merging terms with invalid qubits
-    t1 = HamiltonianTerm(np.random.random((4, 4)), 0, 1)
-    t2 = HamiltonianTerm(np.random.random((4, 4)), 1, 2)
+        ham = hamiltonians.SymbolicHamiltonian("test")
+    # Passing form with symbol that is not in ``symbol_map``
+    from qibo import matrices
+    Z, X = sympy.Symbol("Z"), sympy.Symbol("X")
+    symbol_map = {Z: (0, matrices.Z)}
     with pytest.raises(ValueError):
-        t = t1.merge(t2)
+        ham = hamiltonians.SymbolicHamiltonian(Z * X, symbol_map)
+    # Invalid operation in Hamiltonian expresion
+    ham = hamiltonians.SymbolicHamiltonian(sympy.cos(Z), symbol_map)
+    with pytest.raises(TypeError):
+        dense = ham.dense
 
 
 @pytest.mark.parametrize("nqubits", [3, 4])
-def test_symbolictfim_hamiltonian_to_dense(backend, nqubits):
+@pytest.mark.parametrize("calcterms", [False, True])
+def test_symbolictfim_hamiltonian_to_dense(backend, nqubits, calcterms):
     final_ham = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=1))
     target_ham = hamiltonians.TFIM(nqubits, h=1)
+    if calcterms:
+        _ = final_ham.terms
     K.assert_allclose(final_ham.matrix, target_ham.matrix, atol=1e-15)
 
 
 @pytest.mark.parametrize("nqubits", [3, 4])
-def test_symbolicxxz_hamiltonian_to_dense(backend, nqubits):
+@pytest.mark.parametrize("calcterms", [False, True])
+def test_symbolicxxz_hamiltonian_to_dense(backend, nqubits, calcterms):
     from qibo.symbols import X, Y, Z
     sham = sum(X(i) * X(i + 1) for i in range(nqubits - 1))
     sham += sum(Y(i) * Y(i + 1) for i in range(nqubits - 1))
@@ -57,48 +56,66 @@ def test_symbolicxxz_hamiltonian_to_dense(backend, nqubits):
     sham += X(0) * X(nqubits - 1) + Y(0) * Y(nqubits - 1) + 0.5 * Z(0) * Z(nqubits - 1)
     final_ham = hamiltonians.SymbolicHamiltonian(sham)
     target_ham = hamiltonians.XXZ(nqubits)
+    if calcterms:
+        _ = final_ham.terms
     K.assert_allclose(final_ham.matrix, target_ham.matrix, atol=1e-15)
 
 
+@pytest.mark.parametrize("nqubits", [3])
+@pytest.mark.parametrize("calcterms", [False, True])
 @pytest.mark.parametrize("calcdense", [False, True])
-def test_symbolic_hamiltonian_scalar_mul(backend, calcdense, nqubits=3):
+def test_symbolic_hamiltonian_scalar_mul(backend, nqubits, calcterms, calcdense):
     """Test multiplication of Trotter Hamiltonian with scalar."""
     local_ham = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=1.0))
     target_ham = 2 * hamiltonians.TFIM(nqubits, h=1.0)
+    if calcterms:
+        _ = local_ham.terms
     if calcdense:
         _ = local_ham.dense
     local_dense = (2 * local_ham).dense
     K.assert_allclose(local_dense.matrix, target_ham.matrix)
 
     local_ham = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=1.0))
+    if calcterms:
+        _ = local_ham.terms
     if calcdense:
         _ = local_ham.dense
     local_dense = (local_ham * 2).dense
     K.assert_allclose(local_dense.matrix, target_ham.matrix)
 
 
+@pytest.mark.parametrize("nqubits", [4])
+@pytest.mark.parametrize("calcterms", [False, True])
 @pytest.mark.parametrize("calcdense", [False, True])
-def test_symbolic_hamiltonian_scalar_add(backend, calcdense, nqubits=4):
+def test_symbolic_hamiltonian_scalar_add(backend, nqubits, calcterms, calcdense):
     """Test addition of Trotter Hamiltonian with scalar."""
     local_ham = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=1.0))
     target_ham = 2 + hamiltonians.TFIM(nqubits, h=1.0)
+    if calcterms:
+        _ = local_ham.terms
     if calcdense:
         _ = local_ham.dense
     local_dense = (2 + local_ham).dense
     K.assert_allclose(local_dense.matrix, target_ham.matrix)
 
     local_ham = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=1.0))
+    if calcterms:
+        _ = local_ham.terms
     if calcdense:
         _ = local_ham.dense
     local_dense = (local_ham + 2).dense
     K.assert_allclose(local_dense.matrix, target_ham.matrix)
 
 
+@pytest.mark.parametrize("nqubits", [3])
+@pytest.mark.parametrize("calcterms", [False, True])
 @pytest.mark.parametrize("calcdense", [False, True])
-def test_symbolic_hamiltonian_scalar_sub(backend, calcdense, nqubits=3):
+def test_symbolic_hamiltonian_scalar_sub(backend, nqubits, calcterms, calcdense):
     """Test subtraction of Trotter Hamiltonian with scalar."""
     local_ham = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=1.0))
     target_ham = 2 - hamiltonians.TFIM(nqubits, h=1.0)
+    if calcterms:
+        _ = local_ham.terms
     if calcdense:
         _ = local_ham.dense
     local_dense = (2 - local_ham).dense
@@ -106,17 +123,24 @@ def test_symbolic_hamiltonian_scalar_sub(backend, calcdense, nqubits=3):
 
     target_ham = hamiltonians.TFIM(nqubits, h=1.0) - 2
     local_ham = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=1.0))
+    if calcterms:
+        _ = local_ham.terms
     if calcdense:
         _ = local_ham.dense
     local_dense = (local_ham - 2).dense
     K.assert_allclose(local_dense.matrix, target_ham.matrix)
 
 
+@pytest.mark.parametrize("nqubits", [3])
+@pytest.mark.parametrize("calcterms", [False, True])
 @pytest.mark.parametrize("calcdense", [False, True])
-def test_symbolic_hamiltonian_operator_add_and_sub(backend, calcdense, nqubits=3):
+def test_symbolic_hamiltonian_operator_add_and_sub(backend, nqubits, calcterms, calcdense):
     """Test addition and subtraction between Trotter Hamiltonians."""
     local_ham1 = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=1.0))
     local_ham2 = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=0.5))
+    if calcterms:
+        _ = local_ham1.terms
+        _ = local_ham2.terms
     if calcdense:
         _ = local_ham1.dense
         _ = local_ham2.dense
@@ -128,6 +152,9 @@ def test_symbolic_hamiltonian_operator_add_and_sub(backend, calcdense, nqubits=3
 
     local_ham1 = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=1.0))
     local_ham2 = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=0.5))
+    if calcterms:
+        _ = local_ham1.terms
+        _ = local_ham2.terms
     if calcdense:
         _ = local_ham1.dense
         _ = local_ham2.dense
@@ -138,12 +165,17 @@ def test_symbolic_hamiltonian_operator_add_and_sub(backend, calcdense, nqubits=3
     K.assert_allclose(dense.matrix, target_ham.matrix)
 
 
+@pytest.mark.parametrize("nqubits", [5])
+@pytest.mark.parametrize("calcterms", [False, True])
 @pytest.mark.parametrize("calcdense", [False, True])
-def test_symbolic_hamiltonian_hamiltonianmatmul(backend, calcdense, nqubits=5):
+def test_symbolic_hamiltonian_hamiltonianmatmul(backend, nqubits, calcterms, calcdense):
     local_ham1 = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=1.0))
     local_ham2 = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=0.5))
     dense_ham1 = hamiltonians.TFIM(nqubits, h=1.0)
     dense_ham2 = hamiltonians.TFIM(nqubits, h=0.5)
+    if calcterms:
+        _ = local_ham1.terms
+        _ = local_ham2.terms
     if calcdense:
         _ = local_ham1.dense
         _ = local_ham2.dense
@@ -152,9 +184,10 @@ def test_symbolic_hamiltonian_hamiltonianmatmul(backend, calcdense, nqubits=5):
     K.assert_allclose(local_matmul.matrix, target_matmul.matrix)
 
 
-@pytest.mark.parametrize("density_matrix", [False, True])
 @pytest.mark.parametrize("nqubits", [3, 4])
-def test_symbolic_hamiltonian_matmul(backend, density_matrix, nqubits):
+@pytest.mark.parametrize("density_matrix", [False, True])
+@pytest.mark.parametrize("calcterms", [False, True])
+def test_symbolic_hamiltonian_matmul(backend, nqubits, density_matrix, calcterms):
     if density_matrix:
         from qibo.core.states import MatrixState
         shape = (2 ** nqubits, 2 ** nqubits)
@@ -165,15 +198,20 @@ def test_symbolic_hamiltonian_matmul(backend, density_matrix, nqubits):
         state = VectorState.from_tensor(random_complex(shape))
     local_ham = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=1.0))
     dense_ham = hamiltonians.TFIM(nqubits, h=1.0)
+    if calcterms:
+        _ = local_ham.terms
     local_matmul = local_ham @ state
     target_matmul = dense_ham @ state
     K.assert_allclose(local_matmul, target_matmul)
 
 
-@pytest.mark.parametrize("calcdense", [False, True])
 @pytest.mark.parametrize("nqubits,normalize", [(3, False), (4, False)])
-def test_symbolic_hamiltonian_state_ev(backend, calcdense, nqubits, normalize):
+@pytest.mark.parametrize("calcterms", [False, True])
+@pytest.mark.parametrize("calcdense", [False, True])
+def test_symbolic_hamiltonian_state_ev(backend, nqubits, normalize, calcterms, calcdense):
     local_ham = hamiltonians.SymbolicHamiltonian(symbolic_tfim(nqubits, h=1.0)) + 2
+    if calcterms:
+        _ = local_ham.terms
     if calcdense:
         _ = local_ham.dense
     dense_ham = hamiltonians.TFIM(nqubits, h=1.0) + 2
@@ -190,11 +228,14 @@ def test_symbolic_hamiltonian_state_ev(backend, calcdense, nqubits, normalize):
 
 
 @pytest.mark.parametrize("density_matrix", [False, True])
-def test_symbolic_hamiltonian_abstract_symbol_ev(backend, density_matrix):
+@pytest.mark.parametrize("calcterms", [False, True])
+def test_symbolic_hamiltonian_abstract_symbol_ev(backend, density_matrix, calcterms):
     from qibo.symbols import X, Symbol
     matrix = np.random.random((2, 2))
     form = X(0) * Symbol(1, matrix) + Symbol(0, matrix) * X(1)
     local_ham = hamiltonians.SymbolicHamiltonian(form)
+    if calcterms:
+        _ = local_ham.terms
     if density_matrix:
         state = K.cast(random_complex((4, 4)))
     else:
diff --git a/src/qibo/tests/test_core_hamiltonians_trotter.py b/src/qibo/tests/test_core_hamiltonians_trotter.py
index 09a9ee85a..031df3297 100644
--- a/src/qibo/tests/test_core_hamiltonians_trotter.py
+++ b/src/qibo/tests/test_core_hamiltonians_trotter.py
@@ -104,7 +104,8 @@ def test_trotter_hamiltonian_three_qubit_term(backend):
 
     terms = [HamiltonianTerm(m1, 0, 1, 2), HamiltonianTerm(m2, 2, 3),
              HamiltonianTerm(m3, 1)]
-    ham = hamiltonians.SymbolicHamiltonian.from_terms(terms)
+    ham = hamiltonians.SymbolicHamiltonian()
+    ham.terms = terms
 
     # Test that the `TrotterHamiltonian` dense matrix is correct
     eye = np.eye(2, dtype=m1.dtype)
diff --git a/src/qibo/tests/test_core_terms.py b/src/qibo/tests/test_core_terms.py
index fbc9e677a..a3b707a47 100644
--- a/src/qibo/tests/test_core_terms.py
+++ b/src/qibo/tests/test_core_terms.py
@@ -18,6 +18,24 @@ def test_hamiltonian_term_initialization(backend):
     K.assert_allclose(term.matrix, matrix)
 
 
+def test_hamiltonian_term_initialization_errors():
+    """Test initializing ``HamiltonianTerm`` with wrong parameters."""
+    # Wrong HamiltonianTerm matrix
+    with pytest.raises(TypeError):
+        t = terms.HamiltonianTerm("test", 0, 1)
+    # Passing negative target qubits in HamiltonianTerm
+    with pytest.raises(ValueError):
+        t = terms.HamiltonianTerm("test", 0, -1)
+    # Passing matrix shape incompatible with number of qubits
+    with pytest.raises(ValueError):
+        t = terms.HamiltonianTerm(np.random.random((4, 4)), 0, 1, 2)
+    # Merging terms with invalid qubits
+    t1 = terms.HamiltonianTerm(np.random.random((4, 4)), 0, 1)
+    t2 = terms.HamiltonianTerm(np.random.random((4, 4)), 1, 2)
+    with pytest.raises(ValueError):
+        t = t1.merge(t2)
+
+
 def test_hamiltonian_term_gates(backend):
     """Test gate application of ``HamiltonianTerm``."""
     matrix = np.random.random((4, 4))
@@ -81,7 +99,7 @@ def test_symbolic_term_creation(use_symbols):
     if use_symbols:
         from qibo.symbols import X, Y
         expression = X(0) * Y(1) * X(1)
-        symbol_map = None
+        symbol_map = {}
     else:
         import sympy
         x0, x1, y1 = sympy.symbols("X0 X1 Y1", commutative=False)
@@ -96,6 +114,28 @@ def test_symbolic_term_creation(use_symbols):
     K.assert_allclose(term.matrix_map.get(1)[1], matrices.X)
 
 
+def test_symbolic_term_with_power_creation():
+    """Test creating ``SymbolicTerm`` from sympy expression that contains powers."""
+    from qibo.symbols import X, Z
+    expression = X(0) ** 4 * Z(1) ** 2 * X(2)
+    term = terms.SymbolicTerm(2, expression)
+    assert term.target_qubits == (0, 1, 2)
+    assert len(term.matrix_map) == 3
+    assert term.coefficient == 2
+    K.assert_allclose(term.matrix_map.get(0), 4 * [matrices.X])
+    K.assert_allclose(term.matrix_map.get(1), 2 * [matrices.Z])
+    K.assert_allclose(term.matrix_map.get(2), [matrices.X])
+
+
+def test_symbolic_term_with_imag_creation():
+    """Test creating ``SymbolicTerm`` from sympy expression that contains imaginary coefficients."""
+    from qibo.symbols import Y
+    expression = 3j * Y(0)
+    term = terms.SymbolicTerm(2, expression)
+    assert term.target_qubits == (0,)
+    assert term.coefficient == 6j
+
+
 def test_symbolic_term_matrix(backend):
     """Test matrix calculation of ``SymbolicTerm``."""
     from qibo.symbols import X, Y, Z
diff --git a/src/qibo/tests/test_hamiltonians.py b/src/qibo/tests/test_hamiltonians.py
index abe646b39..1aff1ffa8 100644
--- a/src/qibo/tests/test_hamiltonians.py
+++ b/src/qibo/tests/test_hamiltonians.py
@@ -19,7 +19,7 @@ models_config = [
     ("MaxCut", {"nqubits": 5}, "maxcut_N5.out"),
 ]
 @pytest.mark.parametrize(("model", "kwargs", "filename"), models_config)
-def test_tfim_model_hamiltonian(model, kwargs, filename):
+def test_hamiltonian_models(model, kwargs, filename):
     """Test pre-coded Hamiltonian models generate the proper matrices."""
     from qibo.tests.test_models_variational import assert_regression_fixture
     H = getattr(hamiltonians, model)(**kwargs)
@@ -28,7 +28,8 @@ def test_tfim_model_hamiltonian(model, kwargs, filename):
 
 
 @pytest.mark.parametrize("nqubits", [3, 4])
-def test_maxcut(nqubits):
+@pytest.mark.parametrize("dense,calcterms", [(True, False), (False, False), (False, True)])
+def test_maxcut(nqubits, dense, calcterms):
     size = 2 ** nqubits
     ham = np.zeros(shape=(size, size), dtype=np.complex128)
     for i in range(nqubits):
@@ -42,5 +43,7 @@ def test_maxcut(nqubits):
             M = np.eye(2**nqubits) - h
             ham += M
     target_ham = K.cast(- ham / 2)
-    final_ham = hamiltonians.MaxCut(nqubits)
+    final_ham = hamiltonians.MaxCut(nqubits, dense)
+    if (not dense) and calcterms:
+        _ = final_ham.terms
     K.assert_allclose(final_ham.matrix, target_ham)