Adaptive QPE¶

Download Notebook - adaptive-qpe.ipynb Open in Playpond

This example illustrates a possible implementation of a quantum algorithm in Guppy, showcasing the possible integration, usage and pitfalls of features including:

Control flow
Implicit capture of Python values, especially for referencing global constants

The implemented algorithm is the basic adaptive random walk phase estimation algorithm found at arXiv:2208.04526 (Alg. 1 with sampling circuit from Fig. 1), with an example Hamiltonian and numbers taken from arXiv:2206.12950 (Section V.B).

Setup¶

We begin by specifying the example Hamiltonian in use, in this case \(H = 0.5 * Z\), in two variants:

The unitary oracle \(U(t) = R_Z[-0.5t]\) generated by the Hamiltonian, and
An eigenstate \(\ket{0}\) of the Hamiltonian.

Although this may seem redundant, we will need both for executing the algorithm.

from guppylang import guppy
from guppylang.std.quantum import qubit, crz


@guppy
def oracle(ctrl: qubit, q: qubit, t: float) -> None:
    """Applies a controlled e^-iπHt/2 gate for the example Hamiltonian H = 0.5 * Z."""
    crz(ctrl, q, angle(0.5 * t))


@guppy
def eigenstate() -> qubit:
    """Prepares eigenstate of the example Hamiltonian H = 0.5 * Z."""
    q = qubit()
    x(q)
    return q

Algorithm¶

The algorithm roughly performs the following steps on each iteration, which are outlined below:

Sample from a distribution generated by a circuit involving the eigenstate, an iteration-specific ancilla and the oracle (see Fig. 1 of arXiv:2208.04526).
Based on the outcome, step the phase estimate up or down by a fixed amount.
Decrease the step size by a fixed amount.

Following arXiv:2206.12950 (Fig. 4), we use 24 such iterations, and finally record the phase estimate using result.

import math
from guppylang.std.angles import angle
from guppylang.std.builtins import output
from guppylang.std.quantum import discard, measure, h, rz, x

sqrt_e = math.sqrt(math.e)
sqrt_e_div = math.sqrt((math.e - 1) / math.e)


@guppy
def main() -> None:
    # Pick initial estimate of phase mean and stdv
    # and prepare eigenstate
    mu, sigma = sqrt_e_div, 1 / sqrt_e
    tgt = eigenstate()

    for _ in range(24):
        t = 1 / sigma

        # Sample from the phase estimate distribution
        aux = qubit()
        h(aux)
        rz(aux, angle((sigma - mu) * t))
        oracle(aux, tgt, t)
        h(aux)
        outcome = measure(aux)

        # Adjust estimate and step size
        if outcome:
            mu += sigma / sqrt_e
        else:
            mu -= sigma / sqrt_e
        sigma *= sqrt_e_div

    discard(tgt)
    output("eigenvalue", 2 * mu)

main.check()

The constants sqrt_e and sqrt_e_div are defined in Python outside the function annotated with @guppy. Guppy captures these values implicitly, so they can be referenced directly without wrapping them in comptime.

The arithmetic involving sigma still runs as part of the Guppy program because sigma changes at runtime. A comptime expression may only depend on values known at compile time, so it cannot include the runtime value sigma.

Testing and Evaluation¶

We test the implementation by creating an emulator and allocating the maximum number of qubits: One holding the eigenstate and one ancilla that will be reused. Running the emulator once yields an example result:

emulator = main.emulator(n_qubits=2).with_seed(2)
emulator.run()

EmulatorResult(results=[QsysShot(entries=[('eigenvalue', 1.0948862183455086)])])

However, this is not guaranteed to be the actual phase of our example Hamiltonian, as it is only one example drawn from a distribution. We thus run a larger number of shots, extract the recorded phase from each and visualize them using matplotlib. Doing this reveals the expected phase of the Hamiltonian as the highest peak in the histogram.

# Ensure that matplotlib does not print warnings. Added to avoid failing on CI.
import matplotlib; matplotlib.set_loglevel("critical")
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker

shots = emulator.with_shots(500).run()

fig, ax = plt.subplots(1, 1)
ax.hist([shot.as_dict()["eigenvalue"] for shot in shots], bins=100)
ax.xaxis.set_major_locator(ticker.MultipleLocator(0.5))
ax.set_title("Phase estimate [frequency / phase]")
ax.set_xlabel("Phase estimate")
ax.set_ylabel("Frequency")
plt.show()

../../_images/4b8a91604715f0dff3b6aaa388b20d6abb3621cf3b315563eac26f704a8d1a0a.png