Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Expectation value of observables non diagonal in the computational basis #1489

Merged
merged 20 commits into from
Oct 31, 2024
Merged
Show file tree
Hide file tree
Changes from 17 commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
e756691
feat: implemented expectation for non diagonal observbales
BrunoLiegiBastonLiegi Oct 14, 2024
2ced672
feat: making expectation_from_samples take input circuits
BrunoLiegiBastonLiegi Oct 15, 2024
ed8b641
fix: added coefficients
BrunoLiegiBastonLiegi Oct 15, 2024
db41268
feat: removed expecation.py
BrunoLiegiBastonLiegi Oct 15, 2024
6fd2bf1
test: working on a test
BrunoLiegiBastonLiegi Oct 16, 2024
e4b3591
test: finished to modify the test
BrunoLiegiBastonLiegi Oct 16, 2024
761eaa8
fix: fixing the test
BrunoLiegiBastonLiegi Oct 16, 2024
4cdfe6b
feat: splitted the circuit and frequency cases
BrunoLiegiBastonLiegi Oct 16, 2024
63b1ce7
doc: improved Hamiltonian.exp_from_samples error message
BrunoLiegiBastonLiegi Oct 16, 2024
2b34b5c
feat: added support for power elevation + removed outdated test
BrunoLiegiBastonLiegi Oct 16, 2024
68d4e55
fix: typos in the docstrings
BrunoLiegiBastonLiegi Oct 16, 2024
63fba04
fix: minor fixes to docstring and coverage
BrunoLiegiBastonLiegi Oct 18, 2024
84191d9
doc: added example to advanced code examples
BrunoLiegiBastonLiegi Oct 18, 2024
1ddb2c2
doc: small fix
BrunoLiegiBastonLiegi Oct 22, 2024
e8480b5
fix: various fixes to make powers work properly
BrunoLiegiBastonLiegi Oct 26, 2024
0dfb4a8
fix: small fixes for tests
BrunoLiegiBastonLiegi Oct 26, 2024
9a10b11
doc: changed example observable
BrunoLiegiBastonLiegi Oct 27, 2024
e2d1c05
Update src/qibo/hamiltonians/hamiltonians.py
BrunoLiegiBastonLiegi Oct 29, 2024
ea570bf
feat: changed Z import in docs
BrunoLiegiBastonLiegi Oct 29, 2024
f9943d7
feat: separated the two expectation functions
BrunoLiegiBastonLiegi Oct 29, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file added doc/final_result.npy
Binary file not shown.
33 changes: 29 additions & 4 deletions doc/source/code-examples/advancedexamples.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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.hamiltonians import Z as HZ
BrunoLiegiBastonLiegi marked this conversation as resolved.
Show resolved Hide resolved

circuit = Circuit(4)
circuit.add(gates.H(i) for i in range(4))
Expand All @@ -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 = HZ(4)

result = circuit(nshots=1024)
expectation_value = hamiltonian.expectation_from_samples(result.frequencies())
Expand All @@ -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_samples(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:
Expand Down
156 changes: 130 additions & 26 deletions src/qibo/hamiltonians/hamiltonians.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
"""Module defining Hamiltonian classes."""

from itertools import chain
from typing import Optional
from math import prod
from typing import Optional, Union

import numpy as np
import sympy
Expand Down Expand Up @@ -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)))))
Expand Down Expand Up @@ -552,41 +556,141 @@ 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 _exp_from_circuit(
self, circuit: "Circuit", qubit_map: dict, nshots: int = None
):
"""
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.

Args:
circuit (Circuit): input frequencies.
qubit_map (dict): qubit map.

Returns:
(float): the calculated expectation value.
"""
from qibo import Circuit, gates
BrunoLiegiBastonLiegi marked this conversation as resolved.
Show resolved Hide resolved

if nshots is None: # pragma: no cover
nshots = 1000
BrunoLiegiBastonLiegi marked this conversation as resolved.
Show resolved Hide resolved
rotated_circuits = []
coefficients = []
Z_observables = []
q_maps = []
for term in self.terms:
BrunoLiegiBastonLiegi marked this conversation as resolved.
Show resolved Hide resolved
# 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(m) for m in measurements]
BrunoLiegiBastonLiegi marked this conversation as resolved.
Show resolved Hide resolved
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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I understand this is useful only for hardware execution, right? Since in simulation you could just simulate the circuit and use expectation directly, which should also be more efficient (and exact).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes I had exactly in mind hardware execution for this. But even for simulation, you wouldn't want the exact expectation from the final state (is this what you are suggesting? not sure) as this is supposed to be sample dependent.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you wouldn't want the exact expectation from the final state (is this what you are suggesting? not sure) as this is supposed to be sample dependent.

Yes, above I meant to get the exact expectation from the final state, however you can still get a sample dependent estimate by first simulating the circuit with shots, getting the samples and using expectation_from_samples. I would expect this to be faster than using expectation_from_circuit for simulation, since in the first case the "state preparation" circuit is only simulated once, while expectation_from_circuit will repeat the simulation of the same circuit for every basis.

On hardware, or if you are simulating with noise, it is probably a different story.

Copy link
Contributor Author

@BrunoLiegiBastonLiegi BrunoLiegiBastonLiegi Oct 30, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah yeah, I see, I've been thinking about doing the state preparation first and then just measuring for simulation. To do that though, we should probably move (at least some part of) the expectation_from_circuit inside the backend, which I am not against by the way, but for the moment I decided to keep it simple. In any case, I was assuming that since execute_circuits launches several parallel jobs, if I am not mistaken, the overhead of running several copies of the same simulation would have been reduced, assuming you have many cores available.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I also think we don't need to do anything about this, since it is already available by executing the circuit manually and then calling expectation_from_samples. This is a bit more manual for the user, but it is still fairly simple and it is not the best idea to provide multiple ways (almost equally simple) to do the same thing as it would increase maintenance from our side.

]
return sum(
[
coeff * obs._exp_from_freq(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:
"""
Calculate the expectation value from some frequencies.
The observable has to be diagonal in the computational basis.

Args:
freq (dict): input frequencies.
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.")

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)
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 expectation_from_samples(
self, data: Union[dict, "Circuit"], qubit_map: dict = None, nshots: int = None
BrunoLiegiBastonLiegi marked this conversation as resolved.
Show resolved Hide resolved
):
"""
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(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
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__):
Expand Down
14 changes: 12 additions & 2 deletions src/qibo/hamiltonians/terms.py
Original file line number Diff line number Diff line change
@@ -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:
Expand Down Expand Up @@ -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

Expand Down
64 changes: 35 additions & 29 deletions tests/test_hamiltonians.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -261,34 +262,39 @@ def test_hamiltonian_expectation_errors(backend):
h.expectation("test")


def test_hamiltonian_expectation_from_samples(backend):
@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)
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))

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)

backend.assert_allclose(exp, exp_from_samples, atol=1e-2)


def test_hamiltonian_expectation_from_samples_errors(backend):
Expand Down Expand Up @@ -441,8 +447,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):
Expand Down
10 changes: 0 additions & 10 deletions tests/test_hamiltonians_symbolic.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
6 changes: 2 additions & 4 deletions tests/test_hamiltonians_terms.py
Original file line number Diff line number Diff line change
Expand Up @@ -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])


Expand Down