Expectation values from samples

src/qibo/hamiltonians/abstract.py

@@ -83,6 +83,21 @@ def expectation(self, state, normalize=False):  # pragma: no cover
         """
         raise_error(NotImplementedError)
 
+    @abstractmethod
+    def expectation_from_samples(self, freq, qubit_map=None):  # pragma: no cover
+        """Computes the real expectation value of a diagonal observable given the frequencies when measuring in the computational basis.
+
+        Args:
+            freq (collections.Counter): the keys are the observed values in binary form
+            and the values the corresponding frequencies, that is the number
+            of times each measured value/bitstring appears.
+            qubit_map (tuple): Mapping between frequencies and qubits. If None, [1,...,len(key)]
+
+        Returns:
+            Real number corresponding to the expectation value.
+        """
+        raise_error(NotImplementedError)
+
     @abstractmethod
     def __add__(self, o):  # pragma: no cover
         """Add operator."""

src/qibo/hamiltonians/hamiltonians.py

@@ -2,6 +2,7 @@
 
 from qibo.config import EINSUM_CHARS, log, raise_error
 from qibo.hamiltonians.abstract import AbstractHamiltonian
+from qibo.symbols import Z
 
 
 class Hamiltonian(AbstractHamiltonian):
@@ -133,6 +134,25 @@ def expectation(self, state, normalize=False):
                 "".format(type(state)),
             )
 
+    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.")
+        keys = list(freq.keys())
+        if qubit_map is None:
+            qubit_map = list(range(int(np.log2(len(obs)))))
+        counts = np.array(list(freq.values())) / sum(freq.values())
+        expval = 0
+        kl = len(qubit_map)
+        for j, k in enumerate(keys):
+            index = 0
+            for i in qubit_map:
+                index += int(k[qubit_map.index(i)]) * 2 ** (kl - 1 - i)
+            expval += obs[index, index] * counts[j]
+        return expval
+
     def eye(self, n=None):
         if n is None:
             n = int(self.matrix.shape[0])
@@ -534,6 +554,42 @@ 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):
+        import numpy as np
+
+        terms = self.terms
+        for term in terms:
+            for factor in term.factors:
+                if isinstance(factor, Z) == False:
+                    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.")
+        keys = list(freq.keys())
+        counts = np.array(list(freq.values())) / sum(freq.values())
+        coeff = list(self.form.as_coefficients_dict().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 * float(coeff[j])
+        return expval
+
     def __add__(self, o):
         if isinstance(o, self.__class__):
             if self.nqubits != o.nqubits:

src/qibo/states.py

@@ -231,3 +231,16 @@ def apply_bitflips(self, p0, p1=None):
             )
         noiseless_samples = self.samples()
         return self.backend.apply_bitflips(noiseless_samples, probs)
+
+    def expectation_from_samples(self, observable):
+        """Computes the real expectation value of a diagonal observable from frequencies.
+
+        Args:
+            observable (Hamiltonian/SymbolicHamiltonian): diagonal observable in the computational basis.
+
+        Returns:
+            Real number corresponding to the expectation value.
+        """
+        freq = self.frequencies(binary=True)
+        qubit_map = self.circuit.measurement_gate.qubits
+        return observable.expectation_from_samples(freq, qubit_map)

src/qibo/tests/test_hamiltonians.py

@@ -2,7 +2,9 @@
 import numpy as np
 import pytest
 
-from qibo import hamiltonians
+from qibo import gates, hamiltonians
+from qibo.models import Circuit
+from qibo.symbols import I, Z
 from qibo.tests.utils import random_complex, random_sparse_matrix
 
 
@@ -245,6 +247,37 @@ def test_hamiltonian_expectation_errors(backend):
         h.expectation("test")
 
 
+def test_hamiltonian_expectation_from_samples(backend):
+    """Test Hamiltonian expectation value calculation."""
+    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)
+    freq = result.frequencies(binary=True)
+    Obs0 = hamiltonians.Hamiltonian(3, matrix).expectation_from_samples(
+        freq, qubit_map=None
+    )
+    Obs1 = h1.expectation(result.state())
+
+    backend.assert_allclose(Obs0, Obs1, atol=10 / np.sqrt(nshots))
+
+
+def test_hamiltonian_expectation_from_samples_errors(backend):
+    obs = random_complex((4, 4))
+    h = hamiltonians.Hamiltonian(2, obs, backend=backend)
+    with pytest.raises(NotImplementedError):
+        h.expectation_from_samples(None, qubit_map=None)
+
+
 @pytest.mark.parametrize("dtype", [np.complex128, np.complex64])
 @pytest.mark.parametrize(
     "sparse_type,dense",

src/qibo/tests/test_hamiltonians_symbolic.py

@@ -3,7 +3,9 @@
 import pytest
 import sympy
 
-from qibo import hamiltonians
+from qibo import gates, hamiltonians
+from qibo.models import Circuit
+from qibo.symbols import I, Y, Z
 from qibo.tests.utils import random_complex
 
 
@@ -269,6 +271,36 @@ def test_symbolic_hamiltonian_state_ev(
     backend.assert_allclose(local_ev, target_ev)
 
 
+def test_hamiltonian_expectation_from_samples(backend):
+    """Test Hamiltonian expectation value calculation."""
+    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)
+    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)
+    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))
+
+
+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):

src/qibo/tests/test_states.py

@@ -1,7 +1,9 @@
+import numpy as np
 import pytest
 
-from qibo import gates
+from qibo import gates, hamiltonians
 from qibo.models import Circuit
+from qibo.symbols import I, Z
 
 
 @pytest.mark.parametrize("target", range(5))
@@ -52,3 +54,26 @@ def test_state_representation_max_terms(backend, density_matrix):
             result.symbolic(max_terms=5)
             == "(0.17678+0j)|00000> + (0.17678+0j)|00001> + (0.17678+0j)|00010> + (0.17678+0j)|00011> + (0.17678+0j)|00100> + ..."
         )
+
+
+@pytest.mark.parametrize("density_matrix", [False, True])
+def test_expectation_from_samples(backend, density_matrix):
+
+    obs0 = 2 * Z(0) * Z(1) + Z(0) * Z(2)
+    obs1 = 2 * Z(0) * Z(1) + Z(0) * Z(2) * I(3)
+    h_sym = hamiltonians.SymbolicHamiltonian(obs0, backend=backend)
+    h_dense = hamiltonians.Hamiltonian(3, h_sym.matrix, backend=backend)
+    h1 = hamiltonians.SymbolicHamiltonian(obs1, backend=backend)
+    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))
+    nshots = 10**5
+    result = c(nshots=nshots)
+    expval_sym = result.expectation_from_samples(h_sym)
+    expval_dense = result.expectation_from_samples(h_dense)
+    expval = h1.expectation(result.state())
+    backend.assert_allclose(expval_sym, expval_dense)
+    backend.assert_allclose(expval_sym, expval, atol=10 / np.sqrt(nshots))