Merge pull request #1686 from qiboteam/exp_fix

Support simultaneously diagonalizable terms in ``exp_from_circuit``

Author: renatomello

src/qibo/hamiltonians/hamiltonians.py

@@ -135,7 +135,7 @@ def expectation_from_samples(self, freq, qubit_map=None):
         diag = self.backend.np.reshape(diag, self.nqubits * (2,))
         if qubit_map is None:
             qubit_map = range(self.nqubits)
-        diag = self.backend.np.transpose(diag, axes=qubit_map).ravel()
+        diag = self.backend.np.transpose(diag, qubit_map).ravel()
         # select only the elements with non-zero counts
         diag = diag[[int(state, 2) for state in freq.keys()]]
         counts = self.backend.cast(list(freq.values()), dtype=diag.dtype) / sum(
@@ -373,6 +373,35 @@ def terms(self):
                 self.constant += term.coefficient
         return terms
 
+    @cached_property
+    def diagonal_terms(self) -> list[list[SymbolicTerm]]:
+        """List of terms that can be diagonalized simultaneously, i.e. that
+        commute with each other. In detail each element of the list is a sublist
+        of commuting ``SymbolicTerm``s.
+        """
+        diagonal_terms = []
+        terms = self.terms
+        # loop until there are no more terms
+        while len(terms) > 0:
+            # take the first (remaining) term
+            t0 = terms[0]
+            diagonal_term = [t0]
+            removable_indices = {
+                0,
+            }
+            # look for all the following terms that commute with the
+            # first one and among themselves
+            for i, t1 in enumerate(terms[1:], 1):
+                # commutes with all the other terms -> append it
+                if all(term.commute(t1) for term in diagonal_term):
+                    diagonal_term.append(t1)
+                    removable_indices.add(i)
+            # append the new sublist of commuting terms found
+            diagonal_terms.append(diagonal_term)
+            # remove them from the original terms
+            terms = [term for i, term in enumerate(terms) if i not in removable_indices]
+        return diagonal_terms
+
     @property
     def matrix(self):
         """Returns the full matrix representation.
@@ -502,46 +531,56 @@ def expectation_from_circuit(self, circuit: "Circuit", nshots: int = 1000) -> fl
         from qibo import gates
 
         rotated_circuits = []
-        coefficients = []
         Z_observables = []
         qubit_maps = []
-        for term in self.terms:
-            # store coefficient
-            coefficients.append(term.coefficient)
-            # Only care about non-I terms
-            non_identity_factors = [
-                factor for factor in term.factors if factor.name[0] != "I"
-            ]
-            # build diagonal observable
-            Z_observables.append(
-                SymbolicHamiltonian(
-                    prod(Z(factor.target_qubit) for factor in non_identity_factors),
+        # loop over the terms that can be diagonalized simultaneously
+        for terms in self.diagonal_terms:
+            # for each term that can be diagonalized simultaneously
+            # extract the coefficient, Z observable and qubit map
+            # the basis rotation, instead, will be the same
+            tmp_obs = []
+            measurements = {}
+            for term in terms:
+                # Only care about non-I terms
+                non_identity_factors = []
+                # prepare the measurement basis and append it to the circuit
+                for factor in term.factors:
+                    if factor.name[0] != "I":
+                        non_identity_factors.append(factor)
+                        q = factor.target_qubit
+                        if q not in measurements:
+                            measurements[q] = gates.M(q, basis=factor.gate.__class__)
+
+                # build diagonal observable
+                # including the coefficient
+                symb_obs = prod(
+                    Z(factor.target_qubit) for factor in non_identity_factors
+                )
+                symb_obs = SymbolicHamiltonian(
+                    term.coefficient * symb_obs,
                     nqubits=circuit.nqubits,
                     backend=self.backend,
                 )
-            )
+                tmp_obs.append(symb_obs)
+
             # Get the qubits we want to measure for each term
-            qubit_map = sorted(factor.target_qubit for factor in non_identity_factors)
-            # prepare the measurement basis and append it to the circuit
-            measurements = [
-                gates.M(factor.target_qubit, basis=factor.gate.__class__)
-                for factor in non_identity_factors
-            ]
+            qubit_maps.append(sorted(measurements.keys()))
+            Z_observables.append(tmp_obs)
+
             circ_copy = circuit.copy(True)
-            circ_copy.add(measurements)
+            circ_copy.add(list(measurements.values()))
             rotated_circuits.append(circ_copy)
-            # for mapping the obtained sample frequencies to the original qubits
-            qubit_maps.append(qubit_map)
-        frequencies = [
-            result.frequencies()
-            for result in self.backend.execute_circuits(rotated_circuits, nshots=nshots)
-        ]
-        return sum(
-            coeff * obs.expectation_from_samples(freq, qubit_map)
-            for coeff, freq, obs, qubit_map in zip(
-                coefficients, frequencies, Z_observables, qubit_maps
-            )
-        )
+
+        # execute the circuits
+        results = self.backend.execute_circuits(rotated_circuits, nshots=nshots)
+
+        # construct the expectation value for each diagonal term
+        # and sum all together
+        expval = 0.0
+        for res, obs, qmap in zip(results, Z_observables, qubit_maps):
+            freq = res.frequencies()
+            expval += sum(o.expectation_from_samples(freq, qmap) for o in obs)
+        return self.constant + expval
 
     def expectation_from_samples(self, freq: dict, qubit_map: list = None) -> float:
         """

src/qibo/hamiltonians/terms.py

@@ -208,12 +208,20 @@ def matrix(self):
             where ``ntargets`` is the number of qubits included in the factors
             of this term.
         """
-        matrices = [
-            reduce(self.backend.np.matmul, self.matrix_map.get(q))
-            for q in self.target_qubits
-        ]
+        matrices = list(self.qubit_to_matrix_map.values())
         return complex(self.coefficient) * reduce(self.backend.np.kron, matrices)
 
+    @cached_property
+    def qubit_to_matrix_map(self) -> dict:
+        """Dictionary mapping each qubit to the corresponding global matrix
+        that acts on it, i.e. the product of the matrices of all the factors
+        acting on it.
+        """
+        return {
+            q: reduce(self.backend.np.matmul, self.matrix_map.get(q))
+            for q in self.target_qubits
+        }
+
     def copy(self):
         """Creates a shallow copy of the term with the same attributes."""
         new = self.__class__(self.coefficient)
@@ -236,6 +244,15 @@ def __call__(self, backend, state, nqubits, density_matrix=False):
             )
         return self.coefficient * state
 
+    def commute(self, term) -> bool:
+        """Check whether this term commutes with another term."""
+        for q in set(self.target_qubits).intersection(set(term.target_qubits)):
+            m1 = self.qubit_to_matrix_map.get(q)
+            m2 = term.qubit_to_matrix_map.get(q)
+            if not self.backend.np.all(m1 @ m2 == m2 @ m1):
+                return False
+        return True
+
 
 class TermGroup(list):
     """Collection of multiple :class:`qibo.hamiltonians.terms.HamiltonianTerm` objects.

tests/test_hamiltonians.py

@@ -305,10 +305,13 @@ def test_hamiltonian_expectation_from_samples_with_some_zero_counts(backend, qma
 
 def test_hamiltonian_expectation_from_circuit(backend):
     """Test Hamiltonian expectation value calculation."""
-    backend.set_seed(12)
+    seed = 12
+    backend.set_seed(seed)
 
     nshots = 4 * 10**6
-    observable = I(0) * Z(1) + X(0) * Z(1) + Y(0) * X(2) / 2 - Z(0) * (1 - Y(1)) ** 3
+    observable = (
+        3.14 + I(0) * Z(1) + 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)