diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst index c89102799f..8d1178aa6d 100644 --- a/doc/source/api-reference/qibo.rst +++ b/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 diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst index 82c95771a1..81c712f910 100644 --- a/doc/source/code-examples/advancedexamples.rst +++ b/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. diff --git a/src/qibo/hamiltonians/hamiltonians.py b/src/qibo/hamiltonians/hamiltonians.py index 804f9f261d..351a39f559 100644 --- a/src/qibo/hamiltonians/hamiltonians.py +++ b/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, diff --git a/src/qibo/hamiltonians/models.py b/src/qibo/hamiltonians/models.py index f1171e27f5..1fa030fac6 100644 --- a/src/qibo/hamiltonians/models.py +++ b/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 diff --git a/tests/test_hamiltonians_symbolic.py b/tests/test_hamiltonians_symbolic.py index c9d208c047..dea02fc124 100644 --- a/tests/test_hamiltonians_symbolic.py +++ b/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):