Quantum Amplitude Estimation: From Theory to Finance Applications

Quantum Amplitude Estimation (QAE) is one of the most practically motivated near-term quantum algorithms. It estimates the probability $a = |\langle \chi | \psi \rangle|^2$ that a quantum state $|\psi\rangle$ is in a marked subspace $|\chi\rangle$. Classical Monte Carlo computes such expectations with error $O(1/\sqrt{N})$ using $N$ samples. QAE achieves error $O(1/N)$ using $N$ calls to the state preparation oracle, a quadratic speedup.

The central application in quantum finance is pricing derivatives: the expectation of a payoff function over an underlying asset distribution maps directly to an amplitude estimation problem.

Mathematical Derivation of the QAE Speedup

To understand why QAE offers a quadratic speedup, we need to trace the argument from Grover’s operator through quantum phase estimation to the final accuracy bound.

Classical baseline. Classical Monte Carlo estimates an expectation $\mathbb{E}[f(X)]$ by drawing $N$ independent samples $x_1, \ldots, x_N$ and computing the sample mean $\hat{\mu} = \frac{1}{N}\sum_i f(x_i)$. By the central limit theorem, the standard error is $\sigma / \sqrt{N}$, so achieving error $\epsilon$ requires $N = O(1/\epsilon^2)$ samples.

Quantum setup. Suppose a unitary $\mathcal{A}$ prepares a state that encodes the quantity of interest as an amplitude:

$\mathcal{A}|0\rangle = \sqrt{1-a}\,|\psi_{\text{bad}}\rangle + \sqrt{a}\,|\psi_{\text{good}}\rangle$

where $a \in [0, 1]$ is the amplitude we wish to estimate. The Grover operator $Q$ (defined below) has eigenvalues $e^{\pm 2i\theta_a}$ where $\theta_a = \arcsin(\sqrt{a})$.

Phase estimation precision. Quantum Phase Estimation (QPE) with $m$ counting qubits extracts the phase $\theta_a$ with precision:

$\delta\theta = \frac{\pi}{2^m}$

This requires $O(2^m)$ controlled applications of $Q$, since QPE applies $Q^{2^0}, Q^{2^1}, \ldots, Q^{2^{m-1}}$ for a total of $2^0 + 2^1 + \cdots + 2^{m-1} = 2^m - 1$ oracle calls.

From phase precision to amplitude precision. Since $a = \sin^2(\theta_a)$, the error in $a$ relates to the error in $\theta$ via the chain rule:

$\delta a = \left|\frac{d}{d\theta}\sin^2(\theta)\right| \cdot \delta\theta = |\sin(2\theta_a)| \cdot \frac{\pi}{2^m}$

Since $|\sin(2\theta_a)| \leq 1$, we get $\delta a \leq \pi / 2^m$. Setting $N = 2^m$ oracle queries gives:

$\delta a = O\!\left(\frac{1}{N}\right)$

Compare this to the classical $\delta a = O(1/\sqrt{N})$. For a target accuracy $\epsilon$:

Method	Queries needed	Scaling
Classical Monte Carlo	$O(1/\epsilon^2)$	$N$ queries give $O(1/\sqrt{N})$ error
Quantum Amplitude Estimation	$O(1/\epsilon)$	$N$ queries give $O(1/N)$ error

This is the quadratic speedup. For $\epsilon = 0.001$, classical Monte Carlo needs roughly $10^6$ samples, while QAE needs roughly $10^3$ oracle calls.

The Core Idea: Amplitude and the Grover Operator

Given a circuit $\mathcal{A}$ that prepares $|\psi\rangle = \sqrt{1-a}|0\rangle|\text{good}\rangle^\perp + \sqrt{a}|1\rangle|\text{good}\rangle$, define the Grover operator:

$Q = -\mathcal{A} S_0 \mathcal{A}^\dagger S_\chi$

where $S_0$ reflects about $|0\rangle$ and $S_\chi$ marks the good states. The eigenvalues of $Q$ are $e^{\pm 2i\arcsin(\sqrt{a})}$. Quantum Phase Estimation on $Q$ extracts the phase $\theta_a = \arcsin(\sqrt{a})$, giving $a = \sin^2(\theta_a)$.

The Grover Operator Q in Detail

The Grover operator $Q$ is built from two reflections. Understanding each component is essential for implementing QAE correctly.

$S_\chi$: Reflection about the good subspace. This operator flips the sign of states in the “good” subspace and leaves “bad” states unchanged:

$S_\chi = I - 2|\psi_{\text{good}}\rangle\langle\psi_{\text{good}}|$

In practice, when the good/bad distinction is encoded in a single objective qubit (qubit 0 in Qiskit convention), $S_\chi$ is simply a $Z$ gate on that qubit. The $Z$ gate maps $|1\rangle \to -|1\rangle$ and $|0\rangle \to |0\rangle$, which flips the sign of the “good” component if $|1\rangle$ marks success.

$S_0$: Reflection about $|0\rangle^{\otimes n}$. This operator flips the sign of every computational basis state except $|0\rangle^{\otimes n}$:

$S_0 = 2|0\rangle^{\otimes n}\langle 0|^{\otimes n} - I$

For $n$ qubits, implement it as:

1. Apply $X$ to all $n$ qubits (maps $|0\rangle^{\otimes n}$ to $|1\rangle^{\otimes n}$)
2. Apply a multi-controlled $Z$ gate (flips sign of $|1\rangle^{\otimes n}$ only)
3. Apply $X$ to all $n$ qubits (undo step 1)

This sequence produces $X^{\otimes n} \cdot \text{MCZ} \cdot X^{\otimes n} = 2|0\rangle\langle 0| - I$.

Full operator $Q = -\mathcal{A} S_0 \mathcal{A}^\dagger S_\chi$. The operator applies $S_\chi$ first (mark good states), then $\mathcal{A}^\dagger$ (uncompute state preparation), then $S_0$ (reflect about zero), then $\mathcal{A}$ (re-prepare), with an overall minus sign.

Eigenvalue derivation. The key insight is that $Q$ acts within a two-dimensional subspace spanned by $|\psi_{\text{good}}\rangle$ and $|\psi_{\text{bad}}\rangle$. In this subspace, $Q$ acts as a rotation by angle $2\theta_a$:

$Q = \begin{pmatrix} \cos(2\theta_a) & -\sin(2\theta_a) \\ \sin(2\theta_a) & \cos(2\theta_a) \end{pmatrix}$

where $\theta_a = \arcsin(\sqrt{a})$. The eigenvalues of this rotation matrix are $e^{\pm 2i\theta_a}$, with eigenstates:

$|\psi_\pm\rangle = \frac{1}{\sqrt{2}}\left(|\psi_{\text{good}}\rangle \mp i\,|\psi_{\text{bad}}\rangle\right)$

QPE extracts the phase $2\theta_a / (2\pi)$, from which we recover $a = \sin^2(\theta_a)$.

Canonical QAE with QPE

import numpy as np
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit.library import QFT
from qiskit.primitives import StatevectorSampler as Sampler
from qiskit_algorithms import IterativeAmplitudeEstimation, EstimationProblem
import matplotlib.pyplot as plt
# Note: qiskit-finance examples require: pip install qiskit-finance

def build_state_preparation(a_target):
    """
    Simple A circuit: rotates |0> to sqrt(1-a)|0> + sqrt(a)|1>.
    Objective qubit is qubit 0.
    """
    theta = 2 * np.arcsin(np.sqrt(a_target))
    qc = QuantumCircuit(1)
    qc.ry(theta, 0)
    return qc

def canonical_qae(A_circuit, m_qubits, shots=10000):
    """
    Canonical QAE using QPE with m counting qubits.
    Returns estimated amplitude a_hat.
    """
    n_objective = A_circuit.num_qubits
    counting = QuantumRegister(m_qubits, 'count')
    objective = QuantumRegister(n_objective, 'obj')
    creg = ClassicalRegister(m_qubits)
    qc = QuantumCircuit(counting, objective, creg)

    # Initialize counting qubits
    qc.h(counting)

    # Apply controlled-Q^(2^j) for each counting qubit
    # For this demo, we use A directly and build the Grover operator
    qc.append(A_circuit, objective)

    # Apply A^\dagger for QPE (simplified: just measure A's output)
    # Full QPE-based QAE requires controlled-Q gates; here we sketch structure
    qc.append(QFT(m_qubits, inverse=True), counting)
    qc.measure(counting, creg)

    sampler = Sampler()
    job = sampler.run([(qc, [], shots)])
    result = job.result()
    counts = result[0].data.meas.get_counts()

    # Most frequent bitstring gives estimated phase
    best = max(counts, key=counts.get)
    k = int(best, 2)
    theta_hat = k * np.pi / (2**m_qubits)
    a_hat = np.sin(theta_hat)**2
    return a_hat

# Test with known amplitude
a_true = 0.3
A = build_state_preparation(a_true)
print(f"True amplitude: {a_true}")
print(f"Circuit:\n{A.draw()}")

Iterative QAE (IQAE): No QPE Required

Canonical QAE needs $m$ ancilla qubits and a full QPE circuit with controlled-$Q^{2^j}$ gates, which are deep. Iterative QAE (Suzuki et al., 2020) avoids QPE entirely. It uses a classical outer loop that adaptively chooses the Grover oracle power $k$ and the number of shots, achieving the same $O(1/\epsilon)$ query complexity with shallow circuits.

How IQAE Works: The Algorithm in Detail

The Suzuki et al. (2020) algorithm maintains a confidence interval $[\theta_{\text{lo}}, \theta_{\text{hi}}]$ for the true angle $\theta_a$ and iteratively tightens it. Here is the step-by-step procedure:

Step 1: Initialize. Set $\theta_{\text{lo}} = 0$, $\theta_{\text{hi}} = \pi/2$, and choose the target precision $\epsilon$ and confidence level $\alpha$.

Step 2: Choose the oracle power $k$. In each round $i$, select $k_i$ to be the largest power such that the interval $[\theta_{\text{lo}}, \theta_{\text{hi}}]$ maps injectively under the function $f_k(\theta) = \sin^2((2k+1)\theta)$. Intuitively, applying $Q^k$ amplifies the phase by a factor of $(2k+1)$, which increases sensitivity but introduces aliasing if the interval is too wide. The algorithm picks $k$ to be as large as possible without aliasing.

Step 3: Run the quantum circuit. Prepare $\mathcal{A}|0\rangle$, apply $Q^{k_i}$, and measure the objective qubit. Repeat for $N_i$ shots. Let $h_i$ be the number of times the objective qubit is measured as $|1\rangle$.

Step 4: Update the confidence interval. The probability of measuring $|1\rangle$ after $Q^k$ is $\sin^2((2k+1)\theta_a)$. Using the observed fraction $h_i / N_i$ and a Chernoff-Hoeffding bound, compute a confidence interval for $\sin^2((2k+1)\theta_a)$. Invert through $\arcsin$ to obtain an updated interval for $\theta_a$. Intersect with the previous interval $[\theta_{\text{lo}}, \theta_{\text{hi}}]$.

Step 5: Check convergence. If $\theta_{\text{hi}} - \theta_{\text{lo}} < \epsilon / (2\sin(2\theta_{\text{mid}}))$ (ensuring that the amplitude precision meets the target), stop. Otherwise, return to Step 2.

The schedule of oracle powers roughly follows $k = 1, 2, 4, 8, \ldots, k_{\max}$, doubling each round. Since round $i$ uses $O(N_i)$ shots with $k_i \approx 2^i$, the total number of oracle calls is:

$\sum_i N_i \cdot k_i \approx \sum_i O(1) \cdot 2^i = O(2^{m}) = O(1/\epsilon)$

This achieves the same $O(1/\epsilon)$ total query complexity as canonical QAE, but each individual circuit has depth at most $O(k_i)$ for that round.

from qiskit_algorithms import IterativeAmplitudeEstimation, EstimationProblem

# Build estimation problem
a_true = 0.3
A = build_state_preparation(a_true)

# Mark the |1> state on the objective qubit as "good"
problem = EstimationProblem(
    state_preparation=A,
    objective_qubits=[0],
)

# Run IQAE
iqae = IterativeAmplitudeEstimation(
    epsilon_target=0.01,   # target accuracy
    alpha=0.05,            # confidence level 1 - alpha = 95%
    sampler=Sampler(),
)
result = iqae.estimate(problem)

print(f"True amplitude:      {a_true:.4f}")
print(f"IQAE estimate:       {result.estimation:.4f}")
print(f"95% CI:              [{result.confidence_interval[0]:.4f}, {result.confidence_interval[1]:.4f}]")
print(f"Oracle queries used: {result.num_oracle_queries}")

IQAE’s confidence interval narrows as $O(1/N_{\text{queries}})$, matching the theoretical speedup over Monte Carlo’s $O(1/\sqrt{N})$.

Comparing QAE Variants

Several QAE variants have been developed, each with different trade-offs. The following table summarizes the four main approaches:

Algorithm	Ancilla qubits	Circuit depth	Confidence interval	Key property
Canonical QAE	$m$	$O(2^m)$	$\pi/2^m$	QPE-based, asymptotically optimal
IQAE (Suzuki et al.)	0	$O(1/\epsilon)$	$\epsilon$	No ancillas, adaptive schedule
MLAE (Maximum Likelihood)	0	$O(1/\epsilon)$	$\epsilon$	Uses ML fitting instead of Chernoff bounds
FQAE (Faster QAE)	0	$O(1/\epsilon)$	$\epsilon$	Better constant factors in query count

MLAE (Suzuki et al., 2020) runs circuits with several different Grover powers $k_0, k_1, \ldots, k_L$, collects the measurement outcomes, and fits the amplitude parameter via maximum likelihood estimation. It avoids the Chernoff bound step of IQAE and can achieve tighter estimates for the same number of queries in practice.

from qiskit_algorithms import MaximumLikelihoodAmplitudeEstimation

a_true = 0.3
A = build_state_preparation(a_true)

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

# MLAE with a specified evaluation schedule
# evaluation_schedule is the list of Grover powers k to use
mlae = MaximumLikelihoodAmplitudeEstimation(
    evaluation_schedule=5,   # uses powers k = 0, 1, 2, 4, 8 (geometric schedule)
    sampler=Sampler(),
)
result_mlae = mlae.estimate(problem)

print(f"True amplitude:      {a_true:.4f}")
print(f"MLAE estimate:       {result_mlae.estimation:.4f}")
print(f"95% CI:              [{result_mlae.confidence_interval[0]:.4f}, {result_mlae.confidence_interval[1]:.4f}]")
print(f"Oracle queries used: {result_mlae.num_oracle_queries}")

The choice of variant depends on hardware constraints. IQAE is the most widely used in Qiskit because it adapts to the problem automatically. MLAE can be preferable when you want to fix the circuit schedule in advance (useful for batch execution on cloud hardware). Canonical QAE is mainly of theoretical interest today, since the deep controlled-$Q$ circuits are impractical on near-term devices.

Accuracy vs Circuit Depth Trade-off

The key engineering tension in QAE is between estimation accuracy and circuit depth:

• Canonical QAE with $m$ counting qubits achieves error $\epsilon \approx \pi / 2^m$ but requires circuits of depth $O(2^m)$ due to high-power Grover iterations.
• IQAE achieves error $\epsilon$ with oracle complexity $O(1/\epsilon)$ and maximum circuit depth $O(1/\epsilon)$, matching canonical QAE but spreading the work over many shorter circuits.
• Classical Monte Carlo achieves error $\epsilon$ with $O(1/\epsilon^2)$ circuit evaluations, each of depth $O(1)$.

For near-term hardware with gate error rate $p_{\text{err}}$, circuits deeper than $O(1/p_{\text{err}})$ gates accumulate too much noise. This limits how many Grover iterations can be applied, reducing QAE’s advantage. Error mitigation and fault tolerance are required for the full quantum speedup.

Accuracy vs. Circuit Depth: Printed Comparison Table

import numpy as np

epsilons = [0.1, 0.05, 0.02, 0.01, 0.005]

print(f"{'epsilon':>10} | {'Classical MC':>14} | {'IQAE queries':>14} | {'Max depth (k)':>14} | {'Speedup':>10}")
print("-" * 72)
for eps in epsilons:
    classical = int(np.ceil(1 / eps**2))
    # IQAE total queries: approximately pi/(4*epsilon)
    quantum = int(np.ceil(np.pi / (4 * eps)))
    # Maximum single-circuit depth (largest Grover power)
    max_k = int(np.ceil(1 / (2 * eps)))
    speedup = classical / quantum
    print(f"{eps:>10.3f} | {classical:>14,} | {quantum:>14,} | {max_k:>14,} | {speedup:>10.1f}x")

Expected output:

   epsilon |     Classical MC |   IQAE queries |  Max depth (k) |    Speedup
------------------------------------------------------------------------
     0.100 |              100 |              8 |              5 |      12.5x
     0.050 |              400 |             16 |             10 |      25.0x
     0.020 |            2,500 |             40 |             25 |      62.5x
     0.010 |           10,000 |             79 |             50 |     126.6x
     0.005 |           40,000 |            158 |            100 |     253.2x

The speedup ratio grows as $O(1/\epsilon)$, confirming the quadratic advantage.

import numpy as np

# Illustrate accuracy vs queries
epsilons = [0.1, 0.05, 0.02, 0.01, 0.005]
classical_queries = [int(1 / e**2) for e in epsilons]
quantum_queries = [int(np.pi / (4 * e)) for e in epsilons]  # O(1/epsilon)

plt.figure(figsize=(7, 4))
plt.loglog(epsilons, classical_queries, 'o-', label='Classical MC: O(1/eps^2)')
plt.loglog(epsilons, quantum_queries, 's-', label='QAE: O(1/eps)')
plt.xlabel('Target error epsilon')
plt.ylabel('Number of oracle queries')
plt.title('QAE vs Monte Carlo Query Complexity')
plt.legend()
plt.gca().invert_xaxis()
plt.tight_layout()
plt.savefig('qae_complexity.png', dpi=150)

Oracle Query Count: Empirical Verification

The following code runs IQAE at several target accuracies and records the actual number of oracle queries used. This lets us verify empirically that IQAE scales as $O(1/\epsilon)$.

import numpy as np
from qiskit.primitives import StatevectorSampler as Sampler
from qiskit_algorithms import IterativeAmplitudeEstimation, EstimationProblem

a_true = 0.3
A = build_state_preparation(a_true)

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

epsilons = [0.1, 0.05, 0.01, 0.005]

print(f"{'epsilon':>10} | {'Classical (1/eps^2)':>20} | {'IQAE queries':>14} | {'Speedup':>10}")
print("-" * 64)

for eps in epsilons:
    iqae = IterativeAmplitudeEstimation(
        epsilon_target=eps,
        alpha=0.05,
        sampler=Sampler(),
    )
    result = iqae.estimate(problem)
    classical = int(np.ceil(1 / eps**2))
    speedup = classical / max(result.num_oracle_queries, 1)
    print(f"{eps:>10.3f} | {classical:>20,} | {result.num_oracle_queries:>14,} | {speedup:>10.1f}x")

print("\nIf IQAE scales as O(1/eps), doubling 1/eps should roughly double the query count.")

Confidence Intervals in Depth

IQAE returns a confidence interval $[a_{\text{lo}}, a_{\text{hi}}]$ with the guarantee that:

$\Pr\left[a_{\text{true}} \in [a_{\text{lo}}, a_{\text{hi}}]\right] \geq 1 - \alpha$

The parameter $\alpha$ (default 0.05 for 95% confidence) controls how conservative the interval is. The width of the interval scales as $O(\epsilon_{\text{target}})$, matching the requested precision.

Tightening the Confidence Interval

Reducing $\alpha$ produces a wider interval (higher confidence) for the same $\epsilon_{\text{target}}$. Conversely, reducing $\epsilon_{\text{target}}$ produces a narrower interval at higher query cost.

import numpy as np
from qiskit.primitives import StatevectorSampler as Sampler
from qiskit_algorithms import IterativeAmplitudeEstimation, EstimationProblem

a_true = 0.3
A = build_state_preparation(a_true)
problem = EstimationProblem(state_preparation=A, objective_qubits=[0])

# Compare different confidence levels
for alpha in [0.10, 0.05, 0.01, 0.001]:
    iqae = IterativeAmplitudeEstimation(
        epsilon_target=0.01,
        alpha=alpha,
        sampler=Sampler(),
    )
    result = iqae.estimate(problem)
    ci = result.confidence_interval
    width = ci[1] - ci[0]
    print(f"alpha={alpha:.3f} ({100*(1-alpha):.1f}% CI): "
          f"[{ci[0]:.4f}, {ci[1]:.4f}], width={width:.4f}, "
          f"queries={result.num_oracle_queries}")

Amplitude vs. Expectation Value

An important subtlety: the confidence interval is for the amplitude $a$, not directly for the expectation value of interest. When using QAE for a finance or integration application, the mapping from amplitude to expectation value involves a normalization factor:

$\text{estimate} = a \times \text{normalization\_factor}$

The normalization factor depends on how the payoff function was encoded into the amplitude. For example, in option pricing, Qiskit Finance provides the interpret() method that handles this conversion:

# The raw amplitude from QAE
raw_amplitude = result.estimation

# The interpret() method applies the correct normalization
# to map from amplitude space back to dollar values
estimated_price = european_call.interpret(result.estimation_processed)

The confidence interval must also be mapped through interpret(). Since the mapping is typically linear, the interval transforms straightforwardly:

$\text{CI}_{\text{price}} = [\text{interpret}(a_{\text{lo}}), \text{interpret}(a_{\text{hi}})]$

Note that the confidence interval is not necessarily symmetric around the estimate. The Chernoff bound used in IQAE produces intervals that can be slightly asymmetric, especially near the boundaries $a = 0$ or $a = 1$.

European Call Option Pricing

A European call option pays $\max(S_T - K, 0)$ at maturity, where $S_T$ is the asset price and $K$ is the strike. The fair price is $e^{-rT} \mathbb{E}[\max(S_T - K, 0)]$. Under Black-Scholes, $S_T$ follows a log-normal distribution.

Qiskit Finance encodes the distribution in a quantum register and computes the expected payoff as an amplitude.

# Requires: pip install qiskit-finance
from qiskit_finance.circuit.library import LogNormalDistribution
from qiskit_finance.applications.estimation import EuropeanCallPricing

# Parameters
S0 = 2.0       # initial price
vol = 0.4      # volatility
r = 0.05       # risk-free rate
T = 40/365     # time to maturity
K = 1.896      # strike price
c_approx = 10  # number of qubits for distribution

# Log-normal distribution circuit
mu = (r - 0.5 * vol**2) * T + np.log(S0)
sigma = vol * np.sqrt(T)
mean = np.exp(mu + sigma**2 / 2)
variance = (np.exp(sigma**2) - 1) * np.exp(2 * mu + sigma**2)

low = np.maximum(0, mean - 3 * np.sqrt(variance))
high = mean + 3 * np.sqrt(variance)
num_uncertainty_qubits = 3

uncertainty_model = LogNormalDistribution(
    num_uncertainty_qubits,
    mu=mu,
    sigma=sigma**2,
    bounds=(low, high),
)

# Build the pricing problem
european_call = EuropeanCallPricing(
    num_state_qubits=num_uncertainty_qubits,
    strike_price=K,
    rescaling_factor=0.25,
    bounds=(low, high),
    uncertainty_model=uncertainty_model,
)

problem_finance = european_call.to_estimation_problem()

# Run IQAE on the pricing problem
iqae_finance = IterativeAmplitudeEstimation(
    epsilon_target=0.01,
    alpha=0.05,
    sampler=Sampler(),
)
result_finance = iqae_finance.estimate(problem_finance)

confidence_interval = european_call.interpret(result_finance.confidence_interval_processed)
estimated_price = european_call.interpret(result_finance.estimation_processed)

print(f"Estimated option price:  {estimated_price:.4f}")
print(f"95% confidence interval: [{confidence_interval[0]:.4f}, {confidence_interval[1]:.4f}]")

# Black-Scholes reference
from scipy.stats import norm
d1 = (np.log(S0 / K) + (r + 0.5 * vol**2) * T) / (vol * np.sqrt(T))
d2 = d1 - vol * np.sqrt(T)
bs_price = S0 * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
print(f"Black-Scholes reference: {bs_price:.4f}")

Amplitude Estimation for Numerical Integration

QAE can estimate definite integrals, generalizing the option pricing example to arbitrary functions. The idea is to encode a function $f(x)$ on $[0,1]$ as an amplitude and use QAE to extract the integral $\int_0^1 f(x)\, dx$.

Encoding scheme. Discretize $[0,1]$ into $2^n$ equally spaced points $x_j = j / 2^n$ for $j = 0, \ldots, 2^n - 1$. The state preparation operator $\mathcal{A}$ does two things:

1. Prepares a uniform superposition over $x_j$: $\frac{1}{\sqrt{2^n}} \sum_j |x_j\rangle$
2. For each $|x_j\rangle$, applies a controlled rotation on an ancilla qubit: $|x_j\rangle|0\rangle \to |x_j\rangle\left(\sqrt{1 - f(x_j)}|0\rangle + \sqrt{f(x_j)}|1\rangle\right)$

After these steps, the probability of measuring the ancilla as $|1\rangle$ is:

$a = \frac{1}{2^n} \sum_{j=0}^{2^n - 1} f(x_j) \approx \int_0^1 f(x)\, dx$

QAE estimates this amplitude $a$, giving us the integral estimate.

Important constraint: The function values must satisfy $0 \leq f(x_j) \leq 1$ for the rotation to produce valid amplitudes. For functions with range outside $[0,1]$, normalize by the maximum value and rescale the final answer.

Example: Estimating $\int_0^1 x^2\, dx = 1/3$

import numpy as np
from qiskit import QuantumCircuit
from qiskit.primitives import StatevectorSampler as Sampler
from qiskit_algorithms import IterativeAmplitudeEstimation, EstimationProblem

n_qubits = 4  # 2^4 = 16 discretization points
num_points = 2**n_qubits

def build_integration_circuit(n_qubits):
    """
    Build a circuit that encodes f(x) = x^2 for uniform x in [0,1].
    Uses n_qubits for the x register and 1 ancilla (objective qubit).
    """
    n_total = n_qubits + 1  # x register + objective qubit
    qc = QuantumCircuit(n_total)

    # Step 1: Uniform superposition over x values
    for i in range(n_qubits):
        qc.h(i)

    # Step 2: Controlled rotations to encode f(x) = x^2
    # For each basis state |j>, rotate ancilla by angle 2*arcsin(sqrt(f(x_j)))
    # where x_j = j / 2^n and f(x_j) = x_j^2
    num_points = 2**n_qubits
    for j in range(num_points):
        x_j = j / num_points
        f_val = x_j**2  # f(x) = x^2, already in [0, 1]

        if f_val > 0:
            # Rotation angle for encoding amplitude sqrt(f_val)
            angle = 2 * np.arcsin(np.sqrt(f_val))

            # Build the control pattern for basis state |j>
            # Convert j to binary, apply X gates for 0-bits, then multi-controlled RY
            ctrl_state = format(j, f'0{n_qubits}b')

            # Apply X gates to qubits where ctrl_state has '0'
            for bit_idx, bit in enumerate(reversed(ctrl_state)):
                if bit == '0':
                    qc.x(bit_idx)

            # Multi-controlled RY on the objective qubit (last qubit)
            qc.mcry(angle, list(range(n_qubits)), n_qubits)

            # Undo X gates
            for bit_idx, bit in enumerate(reversed(ctrl_state)):
                if bit == '0':
                    qc.x(bit_idx)

    return qc

A_integration = build_integration_circuit(n_qubits)

# The objective qubit is the last qubit (index n_qubits)
problem_integration = EstimationProblem(
    state_preparation=A_integration,
    objective_qubits=[n_qubits],
)

iqae_integration = IterativeAmplitudeEstimation(
    epsilon_target=0.01,
    alpha=0.05,
    sampler=Sampler(),
)
result_integration = iqae_integration.estimate(problem_integration)

# The estimated amplitude approximates the integral
integral_estimate = result_integration.estimation
exact_value = 1 / 3  # integral of x^2 from 0 to 1

print(f"QAE estimate of integral(x^2, 0, 1): {integral_estimate:.4f}")
print(f"Exact value:                          {exact_value:.4f}")
print(f"Absolute error:                       {abs(integral_estimate - exact_value):.4f}")
print(f"Note: discretization with {num_points} points introduces its own error")

This approach generalizes to any function that can be efficiently encoded as a quantum rotation. The discretization error from using $2^n$ points adds an $O(1/2^n)$ contribution on top of the QAE estimation error.

Portfolio Risk Estimation: VaR and CVaR

Beyond option pricing, QAE applies to portfolio risk metrics. Value at Risk (VaR) and Conditional Value at Risk (CVaR) are standard measures of portfolio downside risk. Both can be formulated as amplitude estimation problems.

Value at Risk (VaR) at confidence level $\beta$ is the loss threshold such that the probability of exceeding it is at most $1 - \beta$:

$\text{VaR}_\beta = \inf\{x : \Pr[\text{Loss} \leq x] \geq \beta\}$

Conditional VaR (CVaR) is the expected loss given that the loss exceeds VaR:

$\text{CVaR}_\beta = \mathbb{E}[\text{Loss} \mid \text{Loss} \geq \text{VaR}_\beta]$

Quantum Approach to Risk Estimation

The quantum encoding works as follows:

1. Encode the portfolio loss distribution as a quantum state, where the amplitudes represent the probability of each discretized loss level.
2. Mark the “bad” outcomes where loss exceeds a threshold $L$. This sets up the objective qubit.
3. Use QAE to estimate the probability of loss exceeding $L$.
4. For VaR: perform a binary search over thresholds $L$ to find the critical value.
5. For CVaR: encode a conditional payoff and estimate its expectation.

import numpy as np
from qiskit import QuantumCircuit
from qiskit.primitives import StatevectorSampler as Sampler
from qiskit_algorithms import IterativeAmplitudeEstimation, EstimationProblem

def build_risk_estimation_circuit(loss_probs, threshold_index):
    """
    Build a circuit for estimating Pr[Loss > threshold].

    Args:
        loss_probs: array of probabilities for each discretized loss level.
                    Must sum to 1. Length must be a power of 2.
        threshold_index: index in loss_probs above which losses are "bad".
    """
    n_qubits = int(np.log2(len(loss_probs)))
    assert 2**n_qubits == len(loss_probs), "Length must be power of 2"

    n_total = n_qubits + 1  # loss register + objective qubit
    qc = QuantumCircuit(n_total)

    # Step 1: Encode the loss distribution
    # Use initialize to prepare the exact probability amplitudes
    amplitudes_full = np.zeros(2**n_total)
    for j, p in enumerate(loss_probs):
        if j >= threshold_index:
            # "Bad" outcome: objective qubit = |1>
            # Index in full register: j (loss) with qubit n_qubits = 1
            idx_bad = j + 2**n_qubits
            amplitudes_full[idx_bad] = np.sqrt(p)
        else:
            # "Good" outcome: objective qubit = |0>
            amplitudes_full[j] = np.sqrt(p)

    # Normalize (should already be normalized if loss_probs sums to 1)
    norm = np.linalg.norm(amplitudes_full)
    if norm > 0:
        amplitudes_full = amplitudes_full / norm

    qc.initialize(amplitudes_full, range(n_total))
    return qc

# Example: simple portfolio with 8 discretized loss levels
n_loss_qubits = 3
num_levels = 2**n_loss_qubits

# Simulated loss distribution (skewed, heavier right tail)
raw_probs = np.array([0.05, 0.15, 0.25, 0.25, 0.15, 0.08, 0.05, 0.02])
loss_probs = raw_probs / raw_probs.sum()  # normalize to valid distribution

# Set threshold: losses at index 5, 6, 7 are "bad"
threshold_index = 5
prob_exceed_classical = loss_probs[threshold_index:].sum()

A_risk = build_risk_estimation_circuit(loss_probs, threshold_index)

problem_risk = EstimationProblem(
    state_preparation=A_risk,
    objective_qubits=[n_loss_qubits],  # objective qubit is the last one
)

iqae_risk = IterativeAmplitudeEstimation(
    epsilon_target=0.01,
    alpha=0.05,
    sampler=Sampler(),
)
result_risk = iqae_risk.estimate(problem_risk)

print(f"Estimated Pr[Loss > threshold]: {result_risk.estimation:.4f}")
print(f"Classical exact probability:     {prob_exceed_classical:.4f}")
print(f"95% CI: [{result_risk.confidence_interval[0]:.4f}, "
      f"{result_risk.confidence_interval[1]:.4f}]")

For a full portfolio risk application, you would replace the simple discrete distribution with one generated from correlated asset models (e.g., using a Gaussian copula or a multivariate normal distribution encoded via quantum circuits). The qiskit-finance library provides utilities for constructing these distributions.

Running on Real Hardware with Error Mitigation

On real quantum hardware, noise degrades QAE estimates. Each application of the Grover operator $Q$ accumulates gate errors, so circuits with large oracle powers $k$ are the most affected. This creates a practical dilemma: higher powers give more information per shot, but noise corrupts the signal.

The Noise-Depth Trade-off

For a gate error rate $p_{\text{err}}$ per two-qubit gate, a circuit with $d$ gates has an overall fidelity of roughly $(1 - p_{\text{err}})^d$. If each Grover iteration uses $g$ gates, then applying $Q^k$ uses $k \cdot g$ gates with fidelity:

$F \approx (1 - p_{\text{err}})^{k \cdot g}$

For current hardware with $p_{\text{err}} \approx 10^{-3}$ and $g \approx 20$ gates per Grover iteration, circuits become unreliable around $k \approx 50$, limiting the achievable accuracy to $\epsilon \approx 1/(2k) \approx 0.01$.

Practical Strategy: Limit Circuit Depth

IQAE’s epsilon_target parameter implicitly controls the maximum circuit depth. On near-term hardware, choose $\epsilon_{\text{target}}$ conservatively to keep circuits shallow enough for the available error rates.

from qiskit.primitives import StatevectorSampler as Sampler
from qiskit_algorithms import IterativeAmplitudeEstimation, EstimationProblem

# For real hardware, use SamplerV2 from qiskit_ibm_runtime:
#
# from qiskit_ibm_runtime import SamplerV2, QiskitRuntimeService
#
# service = QiskitRuntimeService()
# backend = service.least_busy(simulator=False, min_num_qubits=5)
# sampler = SamplerV2(backend)
#
# To limit circuit depth, set epsilon_target to a moderate value.
# On current hardware, epsilon_target=0.05 is often the practical limit.

a_true = 0.3
A = build_state_preparation(a_true)
problem = EstimationProblem(state_preparation=A, objective_qubits=[0])

# Simulate the hardware-limited scenario with a conservative epsilon
iqae_hw = IterativeAmplitudeEstimation(
    epsilon_target=0.05,   # conservative for near-term hardware
    alpha=0.05,
    sampler=Sampler(),     # replace with SamplerV2(backend) for real hardware
)
result_hw = iqae_hw.estimate(problem)

print(f"Estimate (hardware-friendly): {result_hw.estimation:.4f}")
print(f"Oracle queries used:          {result_hw.num_oracle_queries}")
print(f"True value:                   {a_true:.4f}")
print()
print("For real hardware, also consider:")
print("  - Dynamical decoupling to reduce idle-qubit decoherence")
print("  - Twirled readout error extinction (TREX) for measurement error mitigation")
print("  - Zero-noise extrapolation (ZNE) to estimate the noiseless result")

On real devices, additional error mitigation techniques help recover accuracy:

• Measurement error mitigation corrects for readout errors by calibrating the confusion matrix of the measurement qubits.
• Zero-noise extrapolation (ZNE) runs circuits at several artificially increased noise levels and extrapolates back to zero noise.
• Dynamical decoupling inserts identity-equivalent gate sequences during idle periods to suppress decoherence.

These techniques are available through qiskit_ibm_runtime options and can be combined with IQAE to improve results on near-term hardware.

Connection to Quantum Monte Carlo and Beyond

QAE is fundamentally a quantum-accelerated Monte Carlo method. Classical Monte Carlo underpins vast areas of computational science:

• Computational chemistry: estimating molecular energies, reaction rates, and thermodynamic properties via path integral Monte Carlo.
• Materials science: simulating phase transitions, lattice models, and transport properties.
• Financial engineering: pricing exotic derivatives, computing risk metrics, and optimizing portfolios.
• Machine learning: Bayesian inference, sampling from complex distributions, and training generative models.

A quadratic speedup in Monte Carlo would have broad impact across all these fields. However, several challenges remain before QAE delivers practical advantages:

Oracle construction. The hardest part of applying QAE is often building the oracle $\mathcal{A}$ that encodes the distribution of interest. For financial distributions, this requires loading classical probability distributions into quantum states, which can itself be expensive. For physics simulations, encoding the Boltzmann distribution or path integral measure is a major open problem.

Quadratic vs. exponential. The QAE speedup is “only” quadratic. For a problem where classical Monte Carlo takes $10^6$ samples, QAE needs $10^3$ oracle calls, but each oracle call involves a full quantum circuit execution. If the circuit is much more expensive than a classical sample, the crossover point may be very large.

Constant factors. The $O(1/\epsilon)$ scaling of QAE hides constant factors. In practice, IQAE uses roughly $100/\epsilon$ total queries (including shots), which shifts the crossover compared to the idealized $1/\epsilon$.

Despite these challenges, QAE remains one of the strongest candidates for practical quantum advantage in the near to medium term, particularly in finance where the payoff from even modest speedups can be substantial.

Interpreting the Results

QAE encodes the entire distribution in superposition and extracts the expected payoff through interference, rather than sampling individual price paths. The quantum circuit directly encodes the log-normal distribution into amplitudes, and the payoff function is implemented as a controlled rotation.

The result depends critically on the number of qubits used to discretize the distribution. More qubits capture finer price resolution but increase circuit depth. The num_uncertainty_qubits=3 setting above gives a coarse 8-point discretization; real applications would use 10-12 qubits for sufficient accuracy.

The theoretical advantage appears only when the circuit can run with sufficient depth on real hardware. Current error rates require error mitigation or fault-tolerant execution to realize the $O(1/\epsilon)$ query advantage over classical Monte Carlo.

Common Mistakes

Working with QAE involves several subtle pitfalls. Here are the most common mistakes and how to avoid them.

1. Confusing amplitude $a$ with probability $a$. The amplitude $a$ estimated by QAE equals $|\langle 1 | \psi \rangle|^2$, which is the probability of measuring the objective qubit as $|1\rangle$. These are the same quantity only when the state is a simple superposition over computational basis states. In more complex encodings (multi-qubit registers with entanglement), the amplitude estimated by QAE is the probability of the objective qubit being $|1\rangle$, which encodes the quantity of interest through the specific circuit construction. Do not assume the raw amplitude is the final answer without understanding the encoding.

2. Passing a backend directly instead of a Sampler. IQAE and other qiskit_algorithms estimators require a Sampler primitive, not a raw backend. This is a common source of errors when migrating from older Qiskit code:

# WRONG: passing a backend directly
# iqae = IterativeAmplitudeEstimation(epsilon_target=0.01, alpha=0.05, backend=backend)

# CORRECT: wrap in a Sampler
from qiskit.primitives import StatevectorSampler as Sampler
sampler = Sampler()
iqae = IterativeAmplitudeEstimation(epsilon_target=0.01, alpha=0.05, sampler=sampler)

# For real hardware, use SamplerV2:
# from qiskit_ibm_runtime import SamplerV2
# sampler = SamplerV2(backend)

3. Setting epsilon_target too small for near-term hardware. A target accuracy of $\epsilon = 0.001$ requires Grover circuits with $k \approx 500$ iterations. On current hardware with error rates around $10^{-3}$ per gate, such circuits produce noise-dominated results. Start with $\epsilon = 0.05$ on real devices and decrease only if results are stable.

4. Not normalizing the payoff function to $[0, 1]$. Amplitude encoding requires that the function being estimated maps to valid probabilities in $[0, 1]$. If your payoff function has range $[0, C]$ for some $C > 1$, divide by $C$ before encoding and multiply the final result by $C$. Failing to normalize produces invalid quantum states (the rotation angles exceed $\pi/2$, wrapping around and giving incorrect results).

# Example: if payoff ranges from 0 to 100
max_payoff = 100.0
normalized_payoff = payoff / max_payoff  # now in [0, 1]
# ... run QAE on normalized_payoff ...
# Rescale the result
estimated_value = qae_result.estimation * max_payoff

5. Forgetting to call interpret() for finance problems. The raw QAE result is an amplitude in $[0, 1]$. For finance applications built with qiskit-finance, the interpret() method applies the rescaling factor and offset to convert from amplitude space back to the original units (e.g., dollars). Skipping this step gives a number that looks plausible but has the wrong scale.

# WRONG: using the raw amplitude as the price
wrong_price = result_finance.estimation

# CORRECT: apply the interpret() method
correct_price = european_call.interpret(result_finance.estimation_processed)

6. Treating the confidence interval as symmetric. IQAE’s confidence interval $[a_{\text{lo}}, a_{\text{hi}}]$ is generally not symmetric around the point estimate. Near the boundaries ($a$ close to 0 or 1), the interval can be noticeably skewed. Always report both bounds rather than using $\text{estimate} \pm \text{half-width}$. This asymmetry arises because the Chernoff bound constrains each tail independently, and the mapping through $\arcsin$ introduces additional asymmetry.