amplitude-estimation

name	amplitude-estimation
description	A quantum algorithm for estimating the amplitude of a specific state in a quantum superposition, which can be used for various applications such as Monte Carlo simulations and optimization problems. Provides efficient implementations and educational resources for understanding and utilizing amplitude estimation in quantum algorithm development.

Amplitude Estimation (QAE)

Purpose

Quantum Amplitude Estimation (QAE) estimates the success probability $a$ of a quantum circuit $U$ (i.e., $U|0\rangle = \sqrt{a}|\text{good}\rangle + \sqrt{1-a}|\text{bad}\rangle$) with $O(1/\epsilon)$ circuit evaluations — a quadratic speedup over the classical $O(1/\epsilon^2)$ Monte Carlo sampling.

Use this skill when you need to:

• Precisely estimate a probability amplitude output by a quantum circuit.
• Apply QAE as a subroutine in quantum finance, quantum Monte Carlo, or other estimation tasks.

Overview

QAE combines Grover/Amplitude Amplification with Quantum Phase Estimation (QPE):

1. The Grover operator $Q = S_\psi S_f$ acts as a rotation in the 2D space spanned by $|\text{good}\rangle$ and $|\text{bad}\rangle$, with eigenvalue phases $\pm 2\theta$ where $\sin\theta = \sqrt{a}$.
2. A QPE circuit with $d$ phase-register qubits measures $\phi \approx 2\theta/(2\pi)$.
3. The amplitude is recovered via $\hat{a} = \sin^2(\pi \phi)$.

Prerequisites

• Understanding of Grover iteration (oracle + diffuser).
• Understanding of Quantum Phase Estimation (QPE) and the inverse QFT.
• Python: numpy, Circuit, Register.

Using the Provided Implementation

from unitarylab_algorithms import AmplitudeEstimationAlgorithm
from unitarylab.core import Circuit

# Build state preparation U (data register only, no ancilla)
U = Circuit(2, name="PrepU")
U.ry(1.1, 0)
U.cx(0, 1)

algo = AmplitudeEstimationAlgorithm()
result = algo.run(
    U=U,
    good_zero_qubits=[0],   # Good state: qubit 0 in |0>
    d=6,                    # Phase register bits (higher d = better precision)
    backend='torch'
)

print(result['Target amplitude'])     # Estimated probability p
print(result['Phase'])                # Estimated phase phi
print(result['Most likely phase (bits)'])  # Best phase register bit-string
print(result['circuit_path'])         # SVG circuit diagram path
print(result['plot'])                 # List of saved output file dicts

Core Parameters Explained

Parameter	Type	Default	Description
U	Circuit	required	State preparation circuit on data qubits (no ancilla).
good_zero_qubits	List[int]	required	Indices of qubits that must be $
d	int	6	Number of qubits in the phase register. Precision scales as $\delta a \approx \pi/2^d$.
backend	str	'torch'	Simulation backend.
device	str	'cpu'	Device for simulation (e.g. 'cpu' or 'cuda').
dtype	dtype	np.complex128	Numerical dtype for the statevector simulation.

AmplitudeEstimationAlgorithm(text_mode='plain', algo_dir=None) — constructor parameters:

• text_mode: Output formatting style ('plain' or 'legacy').
• algo_dir: Directory for saving results; auto-derived from cwd if None.

Common misunderstandings:

• d controls precision, not the number of Grover iterations. Larger d means more controlled-$Q$ applications in QPE.
• The total qubit count is d + n_data + 1 (phase + data + Grover ancilla).
• result['Target amplitude'] is $\sin^2(\pi\phi)$; do not confuse the phase $\phi$ with the amplitude.

Return Fields

The result dict is built by _build_return_dict(success, circuit_path, filepath, circuit) and contains:

Key	Type	Description
status	str	'ok' on success, 'failed' otherwise.
circuit_path	str	Path to the saved SVG circuit diagram.
plot	list[dict]	List of saved output files: each entry has 'format' (file extension) and 'filename' (full path).
circuit	Circuit	The QPE circuit object.
Target amplitude	float	Estimated success probability $\hat{a} = \sin^2(\pi\phi)$.
Phase	float	Estimated phase $\phi \in [0, 0.5]$ after folding.
Most likely phase (bits)	str	Bit-string of the highest-probability phase register measurement.
Computation time (s)	float	Wall-clock seconds for the quantum simulation step.
Total qubits	int	Total qubit count: d + n_data + 1.

Implementation Architecture

AmplitudeEstimationAlgorithm in algorithm.py organizes the QAE algorithm into five ordered stages within run(), plus a set of self-contained quantum-circuit building helpers.

run(U, good_zero_qubits, d, backend, device, dtype) — Five Stages:

Stage	Code Action	Algorithmic Role
1 — Parameter Setup	Computes n_data, sets total_qubits = d + n_data + 1	Establishes register layout: phase (d qubits) + data (n_data qubits) + Grover ancilla (1 qubit)
2 — Circuit Construction	Builds prepare (data init), calls _grover_operator_from_zero_oracle(U, good_zero_qubits) to get G, then calls _qpe_circuit(G, d, prepare)	Assembles the complete QPE circuit around the Grover operator
3 — Simulation	qpe_circ.execute(backend=backend, device=device, dtype=dtype) → .state → np.asarray(state) → _phase_histogram(statevector, d)	Runs statevector simulation; extracts phase register probability histogram
4 — Classical Post-Processing	Reads top bit-string; computes phi_raw = int(bits, 2) / 2^d; folds to [0, 0.5]; computes est_amp = sin²(π·φ); calls update_output()	Converts QPE phase readout back to amplitude estimate; stores output fields
5 — Export	save_circuit(qpe_circ) + save_txt() → _build_return_dict(True, circuit_path, filename, qpe_circ)	Saves SVG circuit diagram and text report; assembles return dict

Helper Methods:

• _grover_operator_from_zero_oracle(U, good_zero_qubits) — Builds a single Grover iteration G as a Circuit(n_data+1). Calls _phase_oracle_all_zeros (good-state phase flip), then _diffusion_about_prepared_state (reflection about $|\psi\rangle$). Appends a global phase correction (X + kickback ancilla) so that controlled-G functions correctly in QPE.
• _qpe_circuit(G, d, prepare_target) — Constructs the full QPE circuit. Appends state prep to target register, applies H to all phase qubits, applies controlled-G^(2^k) for each phase qubit k, then appends the iQFT via _iqft_circuit(d).
• _iqft_circuit(n, do_swaps) — Builds the inverse QFT circuit directly from Hadamard gates and controlled-phase rotations mcp. Qubit-swap reversal is applied first.
• _phase_oracle_all_zeros(qc, zero_qubits, ancilla) — Phase-kickback oracle: prepares ancilla in $|-\rangle$, X-flips zero_qubits, applies MCX, unflips, unprepares ancilla. Same as amplitude amplification's oracle.
• _diffusion_about_prepared_state(qc, U, data_qubits, ancilla) — Applies U†, calls _phase_oracle_all_zeros on all data qubits, applies U. Same as amplitude amplification's diffuser.
• _phase_histogram(statevector, d) — Marginalizes the full statevector over the non-phase qubits by grouping indices idx % 2^d and summing probabilities. Returns a sorted dict of bit-string → probability.
• _build_return_dict(success, circuit_path, filepath, circuit) — Assembles the final result dict. Sets status to 'ok'/'failed', wraps file paths into plot (list of {'format': ext, 'filename': path} dicts), attaches the circuit object, and merges all self.output fields via update_output().

Data flow: U + good_zero_qubits → _grover_operator_from_zero_oracle() → G → _qpe_circuit() → qpe_circ → execute() → _phase_histogram() → est_amp → update_output() → _build_return_dict().

Understanding the Key Quantum Components

1. State Preparation Operator

The input circuit U prepares the state

$$ U|0\rangle^n = \sqrt{a},|\text{good}\rangle + \sqrt{1-a},|\text{bad}\rangle $$

where the good state is determined by good_zero_qubits all being $|0\rangle$..

2. Grover Operator $Q$

The Grover operator is constructed internally as:

• Phase oracle $S_f$: Flips the phase of the good state using phase kickback from an ancilla qubit prepared in $|-\rangle$.
• Diffuser $S_\psi = 2|\psi\rangle\langle\psi| - I$: Reflects about the prepared state $|\psi\rangle = U|0\rangle$.

In the 2D invariant subspace, $Q$ acts as a rotation by angle $2\theta$: $$Q|\psi_\pm\rangle = e^{\pm i2\theta}|\psi_\pm\rangle, \quad \sin\theta = \sqrt{a}$$

A global phase correction is applied to cancel the $(-1)$ introduced by the diffuser sign convention, so controlled-$Q$ works correctly in QPE.

3. Quantum Phase Estimation (QPE)

A $d$-qubit phase register is prepared in uniform superposition. Controlled powers of $Q$ ($Q^{2^0}, Q^{2^1}, \ldots, Q^{2^{d-1}}$) encode the phase $2\theta$ into the register. The inverse QFT (iQFT) converts phase to binary.

4. Phase Folding and Amplitude Recovery

QPE produces two symmetric peaks at $\phi \approx 2\theta/(2\pi)$ and $1 - 2\theta/(2\pi)$. The code folds to $[0, 0.5]$ via: $$\phi = \min(\phi_{\text{raw}}, 1 - \phi_{\text{raw}})$$ Then: $$\hat{a} = \sin^2(\pi\phi)$$

Theory-to-Code Mapping

README / Theory Concept	Code Object or Location
State $	\psi\rangle = U
Good state: all good_zero_qubits in $	0\rangle$
Grover operator $Q$	_grover_operator_from_zero_oracle(U, good_zero_qubits) → returns G as a Circuit
Global phase correction $(-1)$	X-kick on ancilla appended at end of _grover_operator_from_zero_oracle
QPE phase register ($d$ bits)	phase = list(range(d)) in _qpe_circuit(); Hadamard applied to all d bits
Controlled $Q^{2^k}$ applications	qc.append(U.repeat(2**k), target=target, control=phase[k], control_state='1') in _qpe_circuit()
Inverse QFT (iQFT)	_iqft_circuit(d) — appended to phase qubits in _qpe_circuit()
Phase readout $\phi$	phi_raw = int(best_bits, 2) / 2^d; folded to [0, 0.5]
Amplitude formula $\hat{a} = \sin^2(\pi\phi)$	est_amp = float(np.sin(np.pi * phi) ** 2)
Phase histogram	_phase_histogram(statevector, d) — marginalizes all non-phase qubits
Query complexity $O(1/\epsilon)$	Implicit: d phase bits → $2^d$ controlled-G total applications → accuracy $\pi/2^d$

Notes on encapsulation: The iQFT is built directly from mcp and h gates in _iqft_circuit() rather than calling an external library. The Grover operator and QPE circuit are both fully self-contained Circuit objects that can be inspected or drawn independently. The phase folding and $\sin^2$ recovery are purely classical post-processing steps after simulation.

Mathematical Deep Dive

Grover eigenvalues: $$Q|\psi_\pm\rangle = e^{\pm i2\theta}|\psi_\pm\rangle, \quad |\psi_\pm\rangle = \frac{1}{\sqrt{2}}(|\text{good}\rangle \pm i|\text{bad}\rangle)$$

QPE readout: With $d$ phase bits, the resolution is $\delta\phi = 1/2^d$, leading to amplitude accuracy: $$\delta a \approx \pi \cdot \delta\phi = \frac{\pi}{2^d}$$

Query complexity: QPE with $d$ bits requires $O(2^d)$ queries to $Q$, which equals $O(1/\epsilon)$ for target accuracy $\epsilon$. Classical Monte Carlo requires $O(1/\epsilon^2)$.

Peak success probability: The main QPE peak appears with probability $\geq 4/\pi^2 \approx 0.405$.

Hands-On Example (UnitaryLab)

from unitarylab_algorithms import AmplitudeEstimationAlgorithm
from unitarylab.core import Circuit
import numpy as np

# True success probability: p = cos^2(pi/8) ≈ 0.854
# Single qubit state: |psi> = cos(theta)|0> + sin(theta)|1>
# Good state: qubit 0 in |0>
theta = np.pi / 8
U = Circuit(1, name="prep")
U.ry(2 * theta, 0)  # cos(theta)|0> + sin(theta)|1>; good = |0>, p = cos^2(theta)

algo = AmplitudeEstimationAlgorithm()
result = algo.run(U=U, good_zero_qubits=[0], d=8, backend='torch')

print(f"True p = {np.cos(theta)**2:.6f}")
print(f"QAE estimate = {result['Target amplitude']:.6f}")
print(f"Phase phi   = {result['Phase']:.6f}")
print(f"Best bits   = {result['Most likely phase (bits)']}")
print(f"Status      = {result['status']}")
print(result['plot'])   # list of saved output file dicts

Reference Implementation (Qiskit)

In addition to the UnitaryLab implementation above, the same amplitude estimation idea can also be expressed using Qiskit’s EstimationProblem and AmplitudeEstimation workflow. This section is provided only as a reference example for users who want to compare different software ecosystems. The main implementation path of this skill remains the UnitaryLab version described above.

Example A: Minimal Qiskit Amplitude Estimation Run

from qiskit.circuit import QuantumCircuit
from qiskit.primitives import StatevectorSampler
from qiskit_algorithms.amplitude_estimators import AmplitudeEstimation, EstimationProblem

# State preparation A
# A|0> = cos(theta)|0> + sin(theta)|1>
theta = 0.3
A = QuantumCircuit(1)
A.ry(2 * theta, 0)

problem = EstimationProblem(
    state_preparation=A,
    objective_qubits=0,
)

ae = AmplitudeEstimation(
    num_eval_qubits=4,
    sampler=StatevectorSampler()
)

result = ae.estimate(problem)

print("Grid estimate:", result.estimation)
print("MLE estimate:", result.mle)
print("Confidence interval:", result.confidence_interval)
print("Oracle queries:", result.num_oracle_queries)

Example B: Qiskit with Custom Grover Operator

from qiskit.circuit import QuantumCircuit
from qiskit.circuit.library import GroverOperator
from qiskit.primitives import StatevectorSampler
from qiskit_algorithms.amplitude_estimators import AmplitudeEstimation, EstimationProblem

theta = 0.3
A = QuantumCircuit(1)
A.ry(2 * theta, 0)

custom_Q = GroverOperator(oracle=A)

problem = EstimationProblem(
    state_preparation=A,
    objective_qubits=0,
    grover_operator=custom_Q,
)

ae = AmplitudeEstimation(
    num_eval_qubits=4,
    sampler=StatevectorSampler()
)

result = ae.estimate(problem)

print("Grid estimate:", result.estimation)
print("MLE estimate:", result.mle)
print("Samples:", result.samples)

Other Qiskit AE Variants (Reference)

Qiskit also includes several amplitude-estimation variants beyond the standard QPE-based AmplitudeEstimation:

• IterativeAmplitudeEstimation
  An iterative AE method that does not rely on Quantum Phase Estimation. It uses selected Grover powers to produce an estimate with target error epsilon_target and confidence level 1 - alpha.
  Official reference:
  https://qiskit-community.github.io/qiskit-algorithms/stubs/qiskit_algorithms.IterativeAmplitudeEstimation.html

• MaximumLikelihoodAmplitudeEstimation
  An AE method without phase estimation. It evaluates multiple Grover powers and determines the final result through maximum-likelihood estimation, without requiring additional evaluation qubits.
  Official reference:
  https://qiskit-community.github.io/qiskit-algorithms/stubs/qiskit_algorithms.MaximumLikelihoodAmplitudeEstimation.html

• FasterAmplitudeEstimation
  A faster AE variant that replaces the QPE component with an iterative Grover-search-style procedure. Compared with canonical QAE, it does not require the additional evaluation qubits and yields less complex circuits.
  Official reference:
  https://qiskit-community.github.io/qiskit-algorithms/stubs/qiskit_algorithms.FasterAmplitudeEstimation.html#qiskit_algorithms.FasterAmplitudeEstimation

Minimal Manual Implementation (UnitaryLab)

The following Python skeleton reconstructs the key components of the QAE algorithm and can be adapted into a compatible implementation.

Step 1: Build the iQFT circuit

# Simplified reconstruction — mirrors AmplitudeEstimationAlgorithm._iqft_circuit()
import numpy as np
from unitarylab.core import Circuit

def build_iqft(n: int, do_swaps: bool = True, backend: str = 'torch') -> Circuit:
    qc = Circuit(n, name=f"iQFT_{n}")
    if do_swaps:
        for i in range(n // 2):
            qc.swap(i, n - 1 - i)
    for j in range(n):
        for k in range(j):
            angle = -np.pi / (2 ** (j - k))
            qc.mcp(angle, k, j)   # controlled-phase rotation
        qc.h(j)
    return qc

Step 2: Build the Grover operator with global phase correction

# Simplified reconstruction — mirrors AmplitudeEstimationAlgorithm._grover_operator_from_zero_oracle()

def build_phase_oracle(n_data: int, good_zero_qubits: list, backend: str = 'torch') -> Circuit:
    """Phase-kickback oracle: flips phase of states where good_zero_qubits are all |0>."""
    qc = Circuit(n_data + 1)
    ancilla = n_data
    # 1. Prepare ancilla in |->
    qc.x(ancilla); qc.h(ancilla)
    # 2. Flip controlled qubits; MCX on ancilla
    for q in good_zero_qubits: qc.x(q)
    if len(good_zero_qubits) == 1:
        qc.cx(good_zero_qubits[0], ancilla)
    else:
        qc.mcx(good_zero_qubits, ancilla)
    for q in good_zero_qubits: qc.x(q)
    # 3. Unprepare ancilla
    qc.h(ancilla); qc.x(ancilla)
    return qc

def build_grover_operator(U: Circuit, good_zero_qubits: list, backend: str = 'torch') -> Circuit:
    """Single Grover iteration G = Diffuser ∘ Oracle, with global phase correction."""
    n_data = U.get_num_qubits()
    ancilla = n_data
    data = list(range(n_data))
    qc = Circuit(n_data + 1, name="GroverIter")
    # Oracle
    oracle = build_phase_oracle(n_data, good_zero_qubits, backend)
    qc.append(oracle, list(range(n_data + 1)))
    # Diffuser: U† → all-zeros oracle → U
    qc.append(U.dagger(), data)
    all_zeros_oracle = build_phase_oracle(n_data, data, backend)
    qc.append(all_zeros_oracle, list(range(n_data + 1)))
    qc.append(U, data)
    # Global phase correction (-1) so controlled-G works in QPE
    qc.x(ancilla); qc.h(ancilla)
    qc.x(ancilla)
    qc.h(ancilla); qc.x(ancilla)
    return qc

Step 3: Build the QPE circuit around G

# Simplified reconstruction — mirrors AmplitudeEstimationAlgorithm._qpe_circuit()

def build_qpe_circuit(G: Circuit, d: int,
                      prepare_target=None, backend: str = 'torch') -> Circuit:
    """QPE circuit: controlled-G^(2^k) powers + iQFT on phase register."""
    n_target = G.get_num_qubits()
    qc = Circuit(d + n_target, name=f"QPE_d{d}")
    phase = list(range(d))
    target = list(range(d, d + n_target))

    if prepare_target is not None:
        qc.append(prepare_target, target)

    for q in phase:
        qc.h(q)

    cG = G.control(1)
    for k in range(d):
        power = 2 ** k
        for _ in range(power):
            qc.append(cG, target + [phase[k]])

    iqft = build_iqft(d, do_swaps=True, backend=backend)
    qc.append(iqft, phase)
    return qc

Step 4: Extract amplitude from phase histogram

def phase_histogram(statevector: np.ndarray, d: int) -> dict:
    """Marginalizes all non-phase qubits; groups by phase register index."""
    probs = np.abs(statevector) ** 2
    counts = {}
    modulus = 2 ** d
    for idx, p in enumerate(probs):
        if p < 1e-12: continue
        k = idx % modulus
        bits = format(k, f'0{d}b')
        counts[bits] = counts.get(bits, 0.0) + float(p)
    return dict(sorted(counts.items(), key=lambda kv: kv[1], reverse=True))

def qae_full(U, good_zero_qubits, d, backend='torch'):
    """End-to-end QAE: returns estimated amplitude."""
    G = build_grover_operator(U, good_zero_qubits, backend)
    n_data = U.get_num_qubits()
    prepare = Circuit(n_data + 1)
    prepare.append(U, list(range(n_data)))
    qpe_circ = build_qpe_circuit(G, d, prepare_target=prepare, backend=backend)
    state = qpe_circ.execute()
    sv = np.asarray(state, dtype=complex).reshape(-1)
    histogram = phase_histogram(sv, d)
    best_bits = next(iter(histogram))
    phi_raw = int(best_bits, 2) / (2 ** d)
    phi = min(phi_raw, 1.0 - phi_raw)
    return float(np.sin(np.pi * phi) ** 2), phi, histogram

Debugging Tips

1. Large estimation error: Increase d. Each additional bit halves the phase resolution.
2. Histogram has many small peaks instead of one dominant peak: This is expected for small d; increase d to sharpen the peak.
3. result['Target amplitude'] misses a factor: Ensure the good state condition (all good_zero_qubits in $|0\rangle$) matches your circuit's actual target. A common mistake is using wrong qubit indices.
4. Symmetric peaks: The histogram may show two peaks at $\phi$ and $1-\phi$. The code folds both to the same estimate — this is correct behavior.
5. Circuit size scaling: Total circuit depth scales as $O(2^d)$ due to QPE. For d=10, expect a deep circuit.