diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst index 623c6c1ce0..bbd592f820 100644 --- a/doc/source/code-examples/advancedexamples.rst +++ b/doc/source/code-examples/advancedexamples.rst @@ -2058,7 +2058,7 @@ diagonalised simultaneously, one has to calculate the expectation value starting # 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_samples(circuit) + 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 diff --git a/src/qibo/hamiltonians/hamiltonians.py b/src/qibo/hamiltonians/hamiltonians.py index a1abcbb72e..710a7b528c 100644 --- a/src/qibo/hamiltonians/hamiltonians.py +++ b/src/qibo/hamiltonians/hamiltonians.py @@ -2,7 +2,7 @@ from itertools import chain from math import prod -from typing import Optional, Union +from typing import Optional import numpy as np import sympy @@ -556,26 +556,34 @@ def calculate_dense(self): def expectation(self, state, normalize=False): return Hamiltonian.expectation(self, state, normalize) - def _exp_from_circuit( - self, circuit: "Circuit", qubit_map: dict, nshots: int = None - ): + def expectation_from_circuit( + self, circuit: "Circuit", qubit_map: dict = None, nshots: int = 1000 + ) -> float: """ Calculate the expectation value from a circuit. - This allows for non diagonal observables. Each term of the observable is - treated separately, by measuring in the correct basis and re-executing the - 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 frequencies. - qubit_map (dict): qubit map. + 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 nshots is None: # pragma: no cover - nshots = 1000 + if qubit_map is None: + qubit_map = list(range(self.nqubits)) + rotated_circuits = [] coefficients = [] Z_observables = [] @@ -614,20 +622,20 @@ def _exp_from_circuit( ] return sum( [ - coeff * obs._exp_from_freq(freq, qmap) + coeff * obs.expectation_from_samples(freq, qmap) for coeff, freq, obs, qmap in zip( coefficients, frequencies, Z_observables, q_maps ) ] ) - def _exp_from_freq(self, freq: dict, qubit_map: dict) -> float: + def expectation_from_samples(self, freq: dict, qubit_map: dict = None) -> float: """ - Calculate the expectation value from some frequencies. + Calculate the expectation value from the samples. The observable has to be diagonal in the computational basis. Args: - freq (dict): input frequencies. + freq (dict): input frequencies of the samples. qubit_map (dict): qubit map. Returns: @@ -641,6 +649,9 @@ def _exp_from_freq(self, freq: dict, qubit_map: dict) -> float: NotImplementedError, "Observable is not a Z Pauli string." ) + 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() @@ -660,38 +671,6 @@ def _exp_from_freq(self, freq: dict, qubit_map: dict) -> float: ) return self.backend.np.sum(expvals @ counts.T) + self.constant.real - def expectation_from_samples( - self, data: Union[dict, "Circuit"], qubit_map: dict = None, nshots: int = None - ): - """ - Calculate the expectation value starting from the samples. 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. - - Args: - data (dict | :class:`qibo.models.Circuit`): either the ``dict`` of the - frequencies of the samples or a :class:`qibo.models.Circuit` where to - extract the samples from. A :class:`qibo.models.Circuit` is needed in - case the observable is not diagonal in the computational basis. - qubit_map (dict): qubit map. - nshots (int, optional): number of shots for the expecation value - calculation. This is used only when a :class:`qibo.models.Circuit` - is passed as input. - - Returns: - (float) the computed expectation value. - """ - if qubit_map is None: - qubit_map = list(range(self.nqubits)) - if isinstance(data, dict): - return self._exp_from_freq(data, qubit_map) - return self._exp_from_circuit(data, qubit_map, nshots) - def __add__(self, o): if isinstance(o, self.__class__): if self.nqubits != o.nqubits: diff --git a/tests/test_hamiltonians.py b/tests/test_hamiltonians.py index bcb5e60f30..5d2e5d57eb 100644 --- a/tests/test_hamiltonians.py +++ b/tests/test_hamiltonians.py @@ -262,38 +262,40 @@ def test_hamiltonian_expectation_errors(backend): h.expectation("test") -@pytest.mark.parametrize( - "observable", - [ - 2 * Z(0) * (1 - Z(1)) ** 2 + Z(0) * Z(2), - X(0) * Z(1) + Y(0) * X(2) / 2 - Z(0) * (1 - Y(1)) ** 3, - ], -) -def test_hamiltonian_expectation_from_samples(backend, observable): - """Test Hamiltonian expectation value calculation.""" - backend.set_seed(12) +def non_exact_expectation_test_setup(backend, observable): nqubits = 3 - nshots = 4 * 10**6 c = Circuit(nqubits) for q in range(nqubits): c.add(gates.RX(q, np.random.rand())) H = hamiltonians.SymbolicHamiltonian(observable, nqubits=nqubits, backend=backend) - matrix = H.matrix - diagonal = ( - backend.np.count_nonzero(matrix - backend.np.diag(backend.np.diagonal(matrix))) - == 0 - ) final_state = backend.execute_circuit(c.copy(True)).state() exp = H.expectation(final_state) - if not diagonal: - exp_from_samples = H.expectation_from_samples(c, nshots=nshots) - else: - c.add(gates.M(*range(nqubits))) - freq = backend.execute_circuit(c, nshots=nshots).frequencies() - exp_from_samples = H.expectation_from_samples(freq) + return exp, H, c + + +def test_hamiltonian_expectation_from_samples(backend): + """Test Hamiltonian expectation value calculation.""" + backend.set_seed(12) + + 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)