Fix SymbolicHamiltonian.expectation for different nqubits

doc/source/api-reference/qibo.rst

@@ -851,7 +851,6 @@ values using a sparse Hamiltonian matrix.
 .. autoclass:: qibo.hamiltonians.Hamiltonian
     :members:
     :member-order: bysource
-    :noindex:
 
 
 Symbolic Hamiltonian
@@ -866,7 +865,6 @@ For more information on constructing Hamiltonians using symbols we refer to the
 .. autoclass:: qibo.hamiltonians.SymbolicHamiltonian
     :members:
     :member-order: bysource
-    :noindex:
 
 
 When a :class:`qibo.hamiltonians.SymbolicHamiltonian` is used for time

doc/source/code-examples/advancedexamples.rst

@@ -1897,3 +1897,75 @@ constructing each symbol:
 
     form = Z(0, commutative=True) * Z(1, commutative=True) + Z(1, commutative=True) * Z(2, commutative=True)
     ham = hamiltonians.SymbolicHamiltonian(form)
+
+
+.. _hamexpectation-example:
+
+How to calculate expectation values using samples?
+--------------------------------------------------
+
+It is possible to calculate the expectation value of a :class:`qibo.hamiltonians.Hamiltonian`
+on a given state using the :meth:`qibo.hamiltonians.Hamiltonian.expectation` method,
+which can be called on a state or density matrix. For example
+
+
+.. testcode::
+
+    from qibo import gates
+    from qibo.models import Circuit
+    from qibo.hamiltonians import XXZ
+
+    circuit = Circuit(4)
+    circuit.add(gates.H(i) for i in range(4))
+    circuit.add(gates.CNOT(0, 1))
+    circuit.add(gates.CNOT(1, 2))
+    circuit.add(gates.CNOT(2, 3))
+
+    hamiltonian = XXZ(4)
+
+    result = circuit()
+    expectation_value = hamiltonian.expectation(result.state())
+
+In this example, the circuit will be simulated to obtain the final state vector
+and the corresponding expectation value will be calculated through exact matrix
+multiplication with the Hamiltonian matrix.
+If a :class:`qibo.hamiltonians.SymbolicHamiltonian` is used instead, the expectation
+value will be calculated as a sum of expectation values of local terms, allowing
+calculations of more qubits with lower memory consumption. The calculation of each
+local term still requires the state vector.
+
+When executing a circuit on real hardware, usually only measurements of the state are
+available, not the state vector. Qibo provides :meth:`qibo.hamiltonians.Hamiltonian.expectation_from_samples`
+to allow calculation of expectation values directly from such samples:
+
+
+.. testcode::
+
+    from qibo import gates
+    from qibo.models import Circuit
+    from qibo.hamiltonians import Z
+
+    circuit = Circuit(4)
+    circuit.add(gates.H(i) for i in range(4))
+    circuit.add(gates.CNOT(0, 1))
+    circuit.add(gates.CNOT(1, 2))
+    circuit.add(gates.CNOT(2, 3))
+    circuit.add(gates.M(*range(4)))
+
+    hamiltonian = Z(4)
+
+    result = circuit(nshots=1024)
+    expectation_value = hamiltonian.expectation_from_samples(result.frequencies())
+
+
+This example simulates the circuit similarly to the previous one but calculates
+the expectation value using the frequencies of shots, instead of the exact state vector.
+This can also be invoked directly from the ``result`` object:
+
+.. testcode::
+
+    expectation_value = result.expectation_from_samples(hamiltonian)
+
+
+The expectation from samples currently works only for Hamiltonians that are diagonal in
+the computational basis.

src/qibo/hamiltonians/hamiltonians.py

@@ -1,3 +1,6 @@
+from itertools import chain
+
+import numpy as np
 import sympy
 
 from qibo.config import EINSUM_CHARS, log, raise_error
@@ -88,7 +91,9 @@ def from_symbolic(cls, symbolic_hamiltonian, symbol_map, backend=None):
             "deprecated. Please use `SymbolicHamiltonian` and Qibo symbols "
             "to construct Hamiltonians using symbols."
         )
-        return SymbolicHamiltonian(symbolic_hamiltonian, symbol_map, backend=backend)
+        return SymbolicHamiltonian(
+            symbolic_hamiltonian, symbol_map=symbol_map, backend=backend
+        )
 
     def eigenvalues(self, k=6):
         if self._eigenvalues is None:
@@ -135,8 +140,6 @@ def expectation(self, state, normalize=False):
             )
 
     def expectation_from_samples(self, freq, qubit_map=None):
-        import numpy as np
-
         obs = self.matrix
         if np.count_nonzero(obs - np.diag(np.diagonal(obs))) != 0:
             raise_error(NotImplementedError, "Observable is not diagonal.")
@@ -318,7 +321,7 @@ class SymbolicHamiltonian(AbstractHamiltonian):
             :meth:`qibo.hamiltonians.models.MaxCut` Hamiltonian.
     """
 
-    def __init__(self, form=None, symbol_map={}, backend=None):
+    def __init__(self, form=None, nqubits=None, symbol_map={}, backend=None):
         super().__init__()
         self._form = None
         self._terms = None
@@ -339,6 +342,8 @@ def __init__(self, form=None, symbol_map={}, backend=None):
             self.backend = backend
         if form is not None:
             self.form = form
+        if nqubits is not None:
+            self.nqubits = nqubits
 
     @property
     def dense(self):
@@ -409,9 +414,6 @@ def terms(self):
                     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
 
@@ -454,8 +456,6 @@ def _get_symbol_matrix(self, term):
             Numerical matrix corresponding to the given expression as a numpy
             array of size ``(2 ** self.nqubits, 2 ** self.nqubits).
         """
-        from numpy import eye
-
         if isinstance(term, sympy.Add):
             # symbolic op for addition
             result = sum(
@@ -492,7 +492,7 @@ def _get_symbol_matrix(self, term):
                 if not isinstance(matrix, self.backend.tensor_types):
                     # symbols that do not correspond to quantum operators
                     # for example parameters in the MaxCut Hamiltonian
-                    result = complex(matrix) * eye(2**self.nqubits)
+                    result = complex(matrix) * np.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`.
@@ -502,7 +502,7 @@ def _get_symbol_matrix(self, term):
             # 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) * eye(2**self.nqubits)
+            result = complex(term) * np.eye(2**self.nqubits)
 
         else:
             raise_error(
@@ -522,10 +522,6 @@ def _calculate_dense_from_form(self):
 
     def _calculate_dense_from_terms(self):
         """Calculates equivalent :class:`qibo.core.hamiltonians.Hamiltonian` using the term representation."""
-        from itertools import chain
-
-        import numpy as np
-
         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.")
@@ -555,8 +551,6 @@ def expectation(self, state, normalize=False):
         return Hamiltonian.expectation(self, state, normalize)
 
     def expectation_from_samples(self, freq, qubit_map=None):
-        import numpy as np
-
         terms = self.terms
         for term in terms:
             # pylint: disable=E1101
@@ -745,15 +739,22 @@ def __matmul__(self, o):
 
         if isinstance(o, self.backend.tensor_types):
             rank = len(tuple(o.shape))
-            if rank == 1:  # state vector
-                return self.apply_gates(o)
-            elif rank == 2:  # density matrix
-                return self.apply_gates(o, density_matrix=True)
-            else:
+            if rank not in (1, 2):
                 raise_error(
                     NotImplementedError,
                     "Cannot multiply Hamiltonian with " "rank-{} tensor.".format(rank),
                 )
+            state_qubits = int(np.log2(int(o.shape[0])))
+            if state_qubits != self.nqubits:
+                raise_error(
+                    ValueError,
+                    f"Cannot multiply Hamiltonian on {self.nqubits} qubits to "
+                    f"state of {state_qubits} qubits.",
+                )
+            if rank == 1:  # state vector
+                return self.apply_gates(o)
+            else:  # density matrix
+                return self.apply_gates(o, density_matrix=True)
 
         raise_error(
             NotImplementedError,

src/qibo/hamiltonians/models.py

@@ -191,7 +191,7 @@ def MaxCut(nqubits, dense=True, backend=None):
     smap = {s: (i, matrices.Z) for i, s in enumerate(Z)}
     smap.update({s: (i, v[i]) for i, s in enumerate(V)})
 
-    ham = SymbolicHamiltonian(sham, smap, backend=backend)
+    ham = SymbolicHamiltonian(sham, symbol_map=smap, backend=backend)
     if dense:
         return ham.dense
     return ham

tests/test_hamiltonians_symbolic.py

@@ -35,9 +35,13 @@ def test_symbolic_hamiltonian_errors(backend):
     Z, X = sympy.Symbol("Z"), sympy.Symbol("X")
     symbol_map = {Z: (0, matrices.Z)}
     with pytest.raises(ValueError):
-        ham = hamiltonians.SymbolicHamiltonian(Z * X, symbol_map, backend=backend)
+        ham = hamiltonians.SymbolicHamiltonian(
+            Z * X, symbol_map=symbol_map, backend=backend
+        )
     # Invalid operation in Hamiltonian expresion
-    ham = hamiltonians.SymbolicHamiltonian(sympy.cos(Z), symbol_map, backend=backend)
+    ham = hamiltonians.SymbolicHamiltonian(
+        sympy.cos(Z), symbol_map=symbol_map, backend=backend
+    )
     with pytest.raises(TypeError):
         dense = ham.dense
 
@@ -248,7 +252,7 @@ def test_symbolic_hamiltonian_matmul(backend, nqubits, density_matrix, calcterms
 @pytest.mark.parametrize("nqubits,normalize", [(3, False), (4, False)])
 @pytest.mark.parametrize("calcterms", [False, True])
 @pytest.mark.parametrize("calcdense", [False, True])
-def test_symbolic_hamiltonian_state_ev(
+def test_symbolic_hamiltonian_state_expectation(
     backend, nqubits, normalize, calcterms, calcdense
 ):
     local_ham = (
@@ -272,6 +276,45 @@ def test_symbolic_hamiltonian_state_ev(
     backend.assert_allclose(local_ev, target_ev)
 
 
+@pytest.mark.parametrize("give_nqubits", [False, True])
+@pytest.mark.parametrize("calcterms", [False, True])
+@pytest.mark.parametrize("calcdense", [False, True])
+def test_symbolic_hamiltonian_state_expectation_different_nqubits(
+    backend, give_nqubits, calcterms, calcdense
+):
+    expr = symbolic_tfim(3, h=1.0)
+    if give_nqubits:
+        local_ham = hamiltonians.SymbolicHamiltonian(expr, nqubits=5, backend=backend)
+    else:
+        local_ham = hamiltonians.SymbolicHamiltonian(expr, backend=backend)
+    if calcterms:
+        _ = local_ham.terms
+    if calcdense:
+        _ = local_ham.dense
+
+    dense_ham = hamiltonians.TFIM(3, h=1.0, backend=backend)
+    dense_matrix = np.kron(backend.to_numpy(dense_ham.matrix), np.eye(4))
+    dense_ham = hamiltonians.Hamiltonian(5, dense_matrix, backend=backend)
+
+    if give_nqubits:
+        state = backend.cast(random_complex((2**5,)))
+        local_ev = local_ham.expectation(state)
+        target_ev = dense_ham.expectation(state)
+        backend.assert_allclose(local_ev, target_ev)
+
+        state = random_complex((2**5,))
+        local_ev = local_ham.expectation(state)
+        target_ev = dense_ham.expectation(state)
+        backend.assert_allclose(local_ev, target_ev)
+    else:
+        state = backend.cast(random_complex((2**5,)))
+        with pytest.raises(ValueError):
+            local_ev = local_ham.expectation(state)
+        state = random_complex((2**5,))
+        with pytest.raises(ValueError):
+            local_ev = local_ham.expectation(state)
+
+
 def test_hamiltonian_expectation_from_samples(backend):
     """Test Hamiltonian expectation value calculation."""
     obs0 = 2 * Z(0) * Z(1) + Z(0) * Z(2)
@@ -285,11 +328,11 @@ def test_hamiltonian_expectation_from_samples(backend):
     c.add(gates.RX(3, np.random.rand()))
     c.add(gates.M(0, 1, 2, 3))
     nshots = 10**5
-    result = c(nshots=nshots)
+    result = backend.execute_circuit(c, nshots=nshots)
     freq = result.frequencies(binary=True)
-    Obs0 = h0.expectation_from_samples(freq, qubit_map=None)
-    Obs1 = h1.expectation(result.state())
-    backend.assert_allclose(Obs0, Obs1, atol=10 / np.sqrt(nshots))
+    ev0 = h0.expectation_from_samples(freq, qubit_map=None)
+    ev1 = h1.expectation(result.state())
+    backend.assert_allclose(ev0, ev1, atol=20 / np.sqrt(nshots))
 
 
 def test_hamiltonian_expectation_from_samples_errors(backend):