Expectation value of observables non diagonal in the computational basis

doc/source/code-examples/advancedexamples.rst

@@ -2015,7 +2015,7 @@ to allow calculation of expectation values directly from such samples:
 .. testcode::
 
     from qibo import Circuit, gates
-    from qibo.hamiltonians import Z
+    from qibo import hamiltonians
 
     circuit = Circuit(4)
     circuit.add(gates.H(i) for i in range(4))
@@ -2024,7 +2024,7 @@ to allow calculation of expectation values directly from such samples:
     circuit.add(gates.CNOT(2, 3))
     circuit.add(gates.M(*range(4)))
 
-    hamiltonian = Z(4)
+    hamiltonian = hamiltonians.Z(4)
 
     result = circuit(nshots=1024)
     expectation_value = hamiltonian.expectation_from_samples(result.frequencies())
@@ -2039,8 +2039,33 @@ This can also be invoked directly from the ``result`` object:
     expectation_value = result.expectation_from_samples(hamiltonian)
 
 
-The expectation from samples currently works only for Hamiltonians that are diagonal in
-the computational basis.
+For Hamiltonians that are not diagonal in the computational basis, or that are sum of terms that cannot be
+diagonalised simultaneously, one has to calculate the expectation value starting from the circuit:
+
+.. testcode::
+
+   from qibo.symbols import X, Y, Z
+   from qibo.hamiltonians import SymbolicHamiltonian
+
+   # build the circuit as before
+   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))
+   # but don't add any measurement at the end!
+   # they will be automatically added with the proper basis
+   # while calculating the expectation value
+
+   hamiltonian = SymbolicHamiltonian(3 * Z(2) * (1 - X(1)) ** 2 - (Y(0) * X(3)) / 2, nqubits=4)
+   expectation_value = hamiltonian.expectation_from_circuit(circuit)
+
+What is happening under the hood in this case, is that the expectation value is calculated for each term
+individually by measuring the circuit in the correct (rotated) basis. All the contributions are then
+summed to recover the global expectation value. This means, in particular, that several copies of the
+circuit are parallely executed, one for each term of the Hamiltonian. Note that, at the moment, no
+check is performed to verify whether a subset of the terms could be diagonalised simultaneously, but
+rather each term is treated separately every time.
 
 
 .. _tutorials_transpiler:

src/qibo/hamiltonians/hamiltonians.py

@@ -1,6 +1,7 @@
 """Module defining Hamiltonian classes."""
 
 from itertools import chain
+from math import prod
 from typing import Optional
 
 import numpy as np
@@ -152,7 +153,10 @@ def expectation_from_samples(self, freq, qubit_map=None):
             )
             != 0
         ):
-            raise_error(NotImplementedError, "Observable is not diagonal.")
+            raise_error(
+                NotImplementedError,
+                "Observable is not diagonal. Expectation of non diagonal observables starting from samples is currently supported for `qibo.hamiltonians.hamiltonians.SymbolicHamiltonian` only.",
+            )
         keys = list(freq.keys())
         if qubit_map is None:
             qubit_map = list(range(int(np.log2(len(obs)))))
@@ -552,41 +556,120 @@ def calculate_dense(self):
     def expectation(self, state, normalize=False):
         return Hamiltonian.expectation(self, state, normalize)
 
-    def expectation_from_samples(self, freq, qubit_map=None):
-        terms = self.terms
-        for term in terms:
+    def expectation_from_circuit(
+        self, circuit: "Circuit", qubit_map: dict = None, nshots: int = 1000
+    ) -> float:
+        """
+        Calculate the expectation value from a circuit.
+        This even works for observables not completely diagonal in the computational
+        basis, but only diagonal at a term level in a defined basis. Namely, for
+        an observable of the form :math:``H = \\sum_i H_i``, where each :math:``H_i``
+        consists in a `n`-qubits pauli operator :math:`P_0 \\otimes P_1 \\otimes \\cdots \\otimes P_n`,
+        the expectation value is computed by rotating the input circuit in the suitable
+        basis for each term :math:``H_i`` thus extracting the `term-wise` expectations
+        that are then summed to build the global expectation value.
+        Each term of the observable is treated separately, by measuring in the correct
+        basis and re-executing the circuit.
+
+        Args:
+            circuit (Circuit): input circuit.
+            qubit_map (dict): qubit map, defaults to ``{0: 0, 1: 1, ..., n: n}``
+            nshots (int): number of shots, defaults to 1000.
+
+        Returns:
+            (float): the calculated expectation value.
+        """
+        from qibo import Circuit, gates
+
+        if qubit_map is None:
+            qubit_map = list(range(self.nqubits))
+
+        rotated_circuits = []
+        coefficients = []
+        Z_observables = []
+        q_maps = []
+        for term in self.terms:
+            # store coefficient
+            coefficients.append(term.coefficient)
+            # build diagonal observable
+            Z_observables.append(
+                SymbolicHamiltonian(
+                    prod([Z(q) for q in term.target_qubits]),
+                    nqubits=circuit.nqubits,
+                    backend=self.backend,
+                )
+            )
+            # prepare the measurement basis and append it to the circuit
+            measurements = [
+                gates.M(factor.target_qubit, basis=factor.gate.__class__)
+                for factor in set(term.factors)
+            ]
+            circ_copy = circuit.copy(True)
+            circ_copy.add(measurements)
+            rotated_circuits.append(circ_copy)
+            # map the original qubits into 0, 1, 2, ...
+            q_maps.append(
+                dict(
+                    zip(
+                        [qubit_map[q] for q in term.target_qubits],
+                        range(len(term.target_qubits)),
+                    )
+                )
+            )
+        frequencies = [
+            result.frequencies()
+            for result in self.backend.execute_circuits(rotated_circuits, nshots=nshots)
+        ]
+        return sum(
+            [
+                coeff * obs.expectation_from_samples(freq, qmap)
+                for coeff, freq, obs, qmap in zip(
+                    coefficients, frequencies, Z_observables, q_maps
+                )
+            ]
+        )
+
+    def expectation_from_samples(self, freq: dict, qubit_map: dict = None) -> float:
+        """
+        Calculate the expectation value from the samples.
+        The observable has to be diagonal in the computational basis.
+
+        Args:
+            freq (dict): input frequencies of the samples.
+            qubit_map (dict): qubit map.
+
+        Returns:
+            (float): the calculated expectation value.
+        """
+        for term in self.terms:
             # pylint: disable=E1101
             for factor in term.factors:
                 if not isinstance(factor, Z):
                     raise_error(
                         NotImplementedError, "Observable is not a Z Pauli string."
                     )
-            if len(term.factors) != len(set(term.factors)):
-                raise_error(NotImplementedError, "Z^k is not implemented since Z^2=I.")
+
+        if qubit_map is None:
+            qubit_map = list(range(self.nqubits))
+
         keys = list(freq.keys())
         counts = self.backend.cast(list(freq.values()), self.backend.precision) / sum(
             freq.values()
         )
-        qubits = []
-        for term in terms:
-            qubits_term = []
-            for k in term.target_qubits:
-                qubits_term.append(k)
-            qubits.append(qubits_term)
-        if qubit_map is None:
-            qubit_map = list(range(len(keys[0])))
-        expval = 0
-        for j, q in enumerate(qubits):
-            subk = []
-            expval_q = 0
-            for i, k in enumerate(keys):
-                subk = [int(k[qubit_map.index(s)]) for s in q]
-                expval_k = 1
-                if subk.count(1) % 2 == 1:
-                    expval_k = -1
-                expval_q += expval_k * counts[i]
-            expval += expval_q * self.terms[j].coefficient.real
-        return expval + self.constant.real
+        expvals = []
+        for term in self.terms:
+            qubits = term.target_qubits
+            expvals.extend(
+                [
+                    term.coefficient.real
+                    * (-1) ** [state[qubit_map[q]] for q in qubits].count("1")
+                    for state in keys
+                ]
+            )
+        expvals = self.backend.cast(expvals, dtype=counts.dtype).reshape(
+            len(self.terms), len(freq)
+        )
+        return self.backend.np.sum(expvals @ counts.T) + self.constant.real
 
     def __add__(self, o):
         if isinstance(o, self.__class__):

src/qibo/hamiltonians/terms.py

@@ -1,8 +1,9 @@
 import numpy as np
 import sympy
 
-from qibo import gates
+from qibo import gates, symbols
 from qibo.config import raise_error
+from qibo.symbols import I, X, Y, Z
 
 
 class HamiltonianTerm:
@@ -152,7 +153,16 @@ def __init__(self, coefficient, factors=1, symbol_map={}):
                     factor, pow = factor.args
                     assert isinstance(pow, sympy.Integer)
                     assert isinstance(factor, sympy.Symbol)
-                    pow = int(pow)
+                    # if the symbol is a Pauli (i.e. a qibo symbol) and `pow` is even
+                    # the power is the identity, thus the factor vanishes. Otherwise,
+                    # for an odd exponent, it remains unchanged (i.e. `pow`=1)
+                    if factor.__class__ in (I, X, Y, Z):
+                        if not int(pow) % 2:
+                            factor = sympy.N(1)
+                        else:
+                            pow = 1
+                    else:
+                        pow = int(pow)
                 else:
                     pow = 1
 

tests/test_hamiltonians.py

@@ -4,8 +4,9 @@
 import pytest
 
 from qibo import Circuit, gates, hamiltonians
+from qibo.hamiltonians.hamiltonians import Hamiltonian, SymbolicHamiltonian
 from qibo.quantum_info.random_ensembles import random_density_matrix, random_statevector
-from qibo.symbols import I, Z
+from qibo.symbols import I, X, Y, Z
 
 from .utils import random_sparse_matrix
 
@@ -261,34 +262,41 @@ def test_hamiltonian_expectation_errors(backend):
         h.expectation("test")
 
 
+def non_exact_expectation_test_setup(backend, observable):
+
+    nqubits = 3
+    c = Circuit(nqubits)
+    for q in range(nqubits):
+        c.add(gates.RX(q, np.random.rand()))
+
+    H = hamiltonians.SymbolicHamiltonian(observable, nqubits=nqubits, backend=backend)
+    final_state = backend.execute_circuit(c.copy(True)).state()
+    exp = H.expectation(final_state)
+    return exp, H, c
+
+
 def test_hamiltonian_expectation_from_samples(backend):
     """Test Hamiltonian expectation value calculation."""
     backend.set_seed(12)
-    obs0 = 2 * Z(0) * Z(1) + Z(0) * Z(2)
-    obs1 = 2 * Z(0) * Z(1) + Z(0) * Z(2) * I(3)
-    h0 = hamiltonians.SymbolicHamiltonian(obs0, backend=backend)
-    h1 = hamiltonians.SymbolicHamiltonian(obs1, backend=backend)
-    matrix = backend.to_numpy(h0.matrix)
-    c = Circuit(4)
-    c.add(gates.RX(0, np.random.rand()))
-    c.add(gates.RX(1, np.random.rand()))
-    c.add(gates.RX(2, np.random.rand()))
-    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 = hamiltonians.Hamiltonian(
-        3, matrix, backend=backend
-    ).expectation_from_samples(freq, qubit_map=None)
-    Obs0 = backend.cast(Obs0, dtype=Obs0.dtype)
-
-    Obs1 = h1.expectation(result.state())
-    Obs1 = backend.cast(Obs1, dtype=Obs1.dtype)
-
-    backend.assert_allclose(Obs0, Obs1, atol=10 / np.sqrt(nshots))
+
+    nshots = 4 * 10**6
+    observable = 2 * Z(0) * (1 - Z(1)) ** 2 + Z(0) * Z(2)
+    exp, H, c = non_exact_expectation_test_setup(backend, observable)
+    c.add(gates.M(*range(c.nqubits)))
+    freq = backend.execute_circuit(c, nshots=nshots).frequencies()
+    exp_from_samples = H.expectation_from_samples(freq)
+    backend.assert_allclose(exp, exp_from_samples, atol=1e-2)
+
+
+def test_hamiltonian_expectation_from_circuit(backend):
+    """Test Hamiltonian expectation value calculation."""
+    backend.set_seed(12)
+
+    nshots = 4 * 10**6
+    observable = X(0) * Z(1) + Y(0) * X(2) / 2 - Z(0) * (1 - Y(1)) ** 3
+    exp, H, c = non_exact_expectation_test_setup(backend, observable)
+    exp_from_samples = H.expectation_from_circuit(c, nshots=nshots)
+    backend.assert_allclose(exp, exp_from_samples, atol=1e-2)
 
 
 def test_hamiltonian_expectation_from_samples_errors(backend):
@@ -441,8 +449,8 @@ def construct_hamiltonian():
     H1 = construct_hamiltonian()
     _ = H1.eigenvectors()
 
-    backend.assert_allclose(H.exp(0.5), target_matrix)
-    backend.assert_allclose(H1.exp(0.5), target_matrix)
+    backend.assert_allclose(H.exp(0.5), target_matrix, atol=1e-6)
+    backend.assert_allclose(H1.exp(0.5), target_matrix, atol=1e-6)
 
 
 def test_hamiltonian_energy_fluctuation(backend):

tests/test_hamiltonians_symbolic.py

@@ -330,16 +330,6 @@ def test_hamiltonian_expectation_from_samples(backend):
     backend.assert_allclose(ev0, ev1, atol=20 / np.sqrt(nshots))
 
 
-def test_hamiltonian_expectation_from_samples_errors(backend):
-    obs = [Z(0) * Y(1), Z(0) * Z(1) ** 3]
-    h1 = hamiltonians.SymbolicHamiltonian(obs[0], backend=backend)
-    h2 = hamiltonians.SymbolicHamiltonian(obs[1], backend=backend)
-    with pytest.raises(NotImplementedError):
-        h1.expectation_from_samples(None, qubit_map=None)
-    with pytest.raises(NotImplementedError):
-        h2.expectation_from_samples(None, qubit_map=None)
-
-
 @pytest.mark.parametrize("density_matrix", [False, True])
 @pytest.mark.parametrize("calcterms", [False, True])
 def test_symbolic_hamiltonian_abstract_symbol_ev(backend, density_matrix, calcterms):

tests/test_hamiltonians_terms.py

@@ -125,11 +125,9 @@ def test_symbolic_term_with_power_creation(backend):
 
     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.target_qubits == (2,)
+    assert len(term.matrix_map) == 1
     assert term.coefficient == 2
-    backend.assert_allclose(term.matrix_map.get(0), 4 * [matrices.X])
-    backend.assert_allclose(term.matrix_map.get(1), 2 * [matrices.Z])
     backend.assert_allclose(term.matrix_map.get(2), [matrices.X])