Amplitude estimation

Download Notebook - qae.ipynb

Prepare-and-measure (PaM) amplitude estimation.

Given a state-preparation unitary \(A\) that flags a “good” subspace on a distinguished target qubit,

\[A\ket{0} = \sqrt{1-a}\,\ket{\psi_0}\ket{0}_{\text{target}} + \sqrt{a}\,\ket{\psi_1}\ket{1}_{\text{target}},\]

the square of the amplitude \(\sqrt a\) is the probability \(a \in [0, 1] \subset \mathbb{R}\) of measuring the target in \(\ket{1}\). Note that the above equation allows for the target qubit to be entangled with the rest of the system.

The prepare-and-measure estimator is the simplest possible amplitude estimator: it prepares \(A\ket{0}\), measures the target, and repeats. The estimate is just the fraction of \(\ket{1}\) outcomes,

\[\hat{a} = \frac{\text{count}(1)}{\text{repeat}} \;\xrightarrow[\text{repeat}\to\infty]{}\; a.\]

For any repeat \(\geq 1\), \(\hat a\) is a random variable who’s expectation value is \(a\), with standard deviation (= estimation error) \(\sqrt{a(1-a)/\text{repeat}}\). This carries no quantum speedup (it is equivalent to classical Monte Carlo sampling), but it is a useful baseline and a building block for the more advanced amplitude estimators.

guppyalgos.algorithms.amplitude_amplification.qae exposes two composable Guppy functions:

  • prepare_and_measure_once(state_prep) → bool — prepare \(A\ket{0}\) and measure the target once;

  • prepare_and_measure(state_prep, repeat) → float — repeat that repeat times and return the estimate.

The state-preparation unitary \(A\) is supplied as a higher-order argument, and qubit allocation is handled internally, so both functions drop straight into a larger workflow.

import numpy as np
from guppylang import guppy
from guppylang.std.angles import angle
from guppylang.std.builtins import array, comptime, result
from guppylang.std.quantum import cx, h, qubit, ry

from guppyalgos.algorithms.amplitude_amplification.qae import prepare_and_measure, prepare_and_measure_once

Define a state preparation and estimate its amplitude

We use a toy \(A\) that rotates the target with an ry gate so that \(P(\text{target}=1) = \sin^2(\pi t / 2)\), and puts the register into an independent superposition (so the input is genuinely a superposition, while the target marginal is exactly the amplitude). guppylang’s angle is in half-turns, i.e. 1.0 corresponds to \(\pi\) radians.

The runnable program is a zero-argument main that calls the estimator and records the returned float with result.

n_register = 2
angle_half_turns = 1.0 / 3.0  # a = sin^2(pi * t / 2) = 0.25
repeat = 4000


@guppy
def state_prep(register: array[qubit, comptime(n_register)], target: qubit) -> None:
    # A|0>: rotate the target, place the register in a uniform superposition
    ry(target, angle(comptime(angle_half_turns)))
    for i in range(comptime(n_register)):
        h(register[i])


@guppy
def main() -> None:
    result("estimate", prepare_and_measure(state_prep, comptime(repeat)))


shots = main.emulator(n_qubits=n_register + 1).with_seed(1).with_shots(1).run()
estimate = shots.collated_shots()[0]["estimate"][0]
analytic = float(np.sin(np.pi * angle_half_turns / 2) ** 2)
print(f"estimate = {estimate:.4f}   analytic a = {analytic:.4f}")

estimate = 0.2430   analytic a = 0.2500

Entangled target

Nothing changes when the target is entangled with the register rather than in a product state. Here ry rotates register[0] and a cx copies it onto the target, producing \(\sqrt{1-a}\,\ket{00} + \sqrt{a}\,\ket{11}\). Measuring the target still yields \(\ket{1}\) with probability \(a\).

@guppy
def entangled_prep(register: array[qubit, comptime(1)], target: qubit) -> None:
    # Entangle the target with the register: sqrt(1-a)|00> + sqrt(a)|11>
    ry(register[0], angle(comptime(angle_half_turns)))
    cx(register[0], target)


@guppy
def main_entangled() -> None:
    result("estimate", prepare_and_measure(entangled_prep, comptime(repeat)))


shots = main_entangled.emulator(n_qubits=2).with_seed(2).with_shots(1).run()
estimate = shots.collated_shots()[0]["estimate"][0]
print(f"entangled estimate = {estimate:.4f}   analytic a = {analytic:.4f}")

entangled estimate = 0.2485   analytic a = 0.2500

The single-shot primitive

prepare_and_measure is just a loop over prepare_and_measure_once, which returns a single target measurement. Recording that per shot and averaging classically over the emulator’s shots reproduces the same estimate — useful when you want the raw outcomes (e.g. to feed a different estimator).

@guppy
def single_shot() -> None:
    result("outcome", prepare_and_measure_once(state_prep))


n_shots = 4000
shots = (
    single_shot.emulator(n_qubits=n_register + 1).with_seed(3).with_shots(n_shots).run()
)
outcomes = [s["outcome"][0] for s in shots.collated_shots()]
print(f"fraction of |1> over {n_shots} shots = {sum(outcomes) / n_shots:.4f}")

fraction of |1> over 4000 shots = 0.2325

Convergence

The sampling error shrinks like \(\sqrt{a(1-a)/\text{repeat}}\), so more repetitions give a tighter estimate. Note this is the classical Monte Carlo \(1/\sqrt{N}\) scaling — the quantum-accelerated amplitude estimators (which reuse the same \(A\)) improve on it.

for r in [50, 500, 5000]:

    @guppy
    def converge() -> None:
        result("estimate", prepare_and_measure(state_prep, comptime(r)))

    est = (
        converge.emulator(n_qubits=n_register + 1)
        .with_seed(0)
        .with_shots(1)
        .run()
        .collated_shots()[0]["estimate"][0]
    )
    std = float(np.sqrt(analytic * (1 - analytic) / r))
    print(f"repeat={r:5d}   estimate={est:.4f}   |error|={abs(est - analytic):.4f}   ~std={std:.4f}")

repeat=   50   estimate=0.2800   |error|=0.0300   ~std=0.0612
repeat=  500   estimate=0.2560   |error|=0.0060   ~std=0.0194
repeat= 5000   estimate=0.2558   |error|=0.0058   ~std=0.0061