Chapter 16: Quantum Phase Estimation — The Subroutine That Powers Shor's, Simulation, and Half of Quantum Computing

Learning Objectives

By the end of this chapter, you will be able to:

• Formulate the eigenvalue estimation problem and explain why it is classically hard.
• Construct the quantum phase estimation (QPE) circuit from controlled-$U$ operations and the inverse QFT.
• Analyze the precision and success probability of QPE as a function of the number of ancilla qubits.
• Derive the measurement probability distribution for both exact and inexact phases.
• Implement QPE in Qiskit for a simple unitary and extract its eigenvalues.
• Explain how QPE serves as the engine for Shor's algorithm, the HHL linear systems solver, and quantum simulation.
• Connect QPE to Hamiltonian simulation via $U = e^{-iHt}$.
• Implement iterative phase estimation (IPE) and compare it with standard QPE.
• Analyze error sources and mitigation strategies for QPE on noisy hardware.
• Understand the role of QPE in the broader landscape of quantum algorithms.

16.1 The Eigenvalue Estimation Problem

16.1.1 Problem Statement

Let $U$ be a unitary operator acting on $m$ qubits. Let $|\psi\rangle$ be an eigenvector of $U$ with eigenvalue $e^{2\pi i \phi}$:

$$U |\psi\rangle = e^{2\pi i \phi} |\psi\rangle, \quad \phi \in [0, 1)$$

The quantum phase estimation problem is: given a quantum circuit that implements controlled-$U^{2^j}$ for $j = 0, 1, \ldots, t-1$, and given the ability to prepare the eigenstate $|\psi\rangle$, estimate the phase $\phi$ to $t$ bits of precision.

Why this chapter matters. Quantum phase estimation is arguably the single most important subroutine in quantum computing. Shor's algorithm, quantum simulation, the HHL algorithm for linear systems, amplitude estimation, and many other algorithms all reduce to phase estimation. If you understand QPE deeply, you understand the engine that drives most quantum algorithms with provable speedups. This is not an exaggeration — QPE appears in roughly half of all known quantum algorithms with exponential or polynomial speedups.

16.1.2 Why Is This Hard Classically?

Classically, to estimate an eigenvalue of a $2^m \times 2^m$ unitary matrix, one must either:

• Diagonalize the matrix: This requires $O(2^{3m})$ operations — exponential in the number of qubits.
• Use the power method: Repeatedly apply $U$ to a trial vector and extract the dominant eigenvalue. This converges in $O(1/\Delta)$ iterations where $\Delta$ is the eigenvalue gap, and each iteration costs $O(2^{2m})$ for dense matrices.
• Use the Lanczos algorithm: More efficient for sparse matrices, but still exponential in $m$ for general unitaries.

Quantum phase estimation achieves this in $O(t^2)$ gates plus the cost of implementing controlled-$U$, which for many physically relevant unitaries (e.g., $U = e^{-iHt}$) is polynomial in $m$. This is an exponential speedup in the dimension of the Hilbert space.

Recurring Theme: Quantum is Linear Algebra, Not Magic

QPE exploits a fundamental property of linear algebra on high-dimensional spaces. An $m$-qubit system lives in a $2^m$-dimensional Hilbert space. The quantum computer naturally manipulates vectors in this space — it doesn't need to store $2^m$ amplitudes explicitly. The QFT operates on this exponentially large space using only $O(m^2)$ gates. The speedup comes from this dimensionality mismatch: the quantum computer operates on the full Hilbert space, while classical algorithms must enumerate its dimensions.

16.1.3 The Phase Kickback Trick

The core mechanism of QPE is phase kickback. Consider a single control qubit in state $|+\rangle = (|0\rangle + |1\rangle)/\sqrt{2}$ and a target eigenstate $|\psi\rangle$:

$$\text{controlled-}U \left( \frac{|0\rangle + |1\rangle}{\sqrt{2}} \otimes |\psi\rangle \right) = \frac{1}{\sqrt{2}} \left( |0\rangle|\psi\rangle + |1\rangle e^{2\pi i \phi} |\psi\rangle \right)$$

$$= \frac{1}{\sqrt{2}} \left( |0\rangle + e^{2\pi i \phi} |1\rangle \right) \otimes |\psi\rangle$$

The phase $e^{2\pi i \phi}$ has been "kicked back" from the target to the control qubit. The target eigenstate $|\psi\rangle$ is unchanged — it factors out of the expression. This is a special property of eigenstates: the controlled-$U$ gate shifts the phase from the target register to the control register.

Try It Yourself: Phase Kickback with Z

Consider $U = Z$ (the Pauli-Z gate). Its eigenstates are $|0\rangle$ (eigenvalue $+1$, phase $\phi = 0$) and $|1\rangle$ (eigenvalue $-1$, phase $\phi = 1/2$). Apply controlled-$Z$ with control in $|+\rangle$ and target in $|1\rangle$:

$$\text{CZ}(|+\rangle \otimes |1\rangle) = \frac{1}{\sqrt{2}}(|0\rangle|1\rangle + |1\rangle(-1)|1\rangle) = |-\rangle \otimes |1\rangle$$

The phase $-1 = e^{2\pi i \cdot 1/2}$ has been kicked back to the control qubit, flipping it from $|+\rangle$ to $|-\rangle$. Measuring the control qubit in the X basis reveals the phase $\phi = 1/2$.

By using $t$ control qubits with powers $U^{2^j}$, we encode $\phi$ into a $t$-qubit Fourier state. The key insight is that different powers of $U$ kick back different phase factors:

• $U^{2^0}$ kicks back $e^{2\pi i \phi \cdot 2^0} = e^{2\pi i \phi}$
• $U^{2^1}$ kicks back $e^{2\pi i \phi \cdot 2^1} = e^{2\pi i \cdot 2\phi}$
• $U^{2^j}$ kicks back $e^{2\pi i \phi \cdot 2^j} = e^{2\pi i \cdot 2^j \phi}$

Each qubit encodes one binary digit of $\phi$: the $j$-th qubit (from the least significant) encodes the $j$-th bit of the binary expansion of $\phi$. This is exactly what the QFT converts from the Fourier basis to the computational basis.

16.2 The Phase Estimation Circuit

16.2.1 Circuit Structure

The QPE circuit uses two registers:

• Register 1 (estimation register): $t$ qubits, initialized to $|0\rangle^{\otimes t}$.
• Register 2 (eigenstate register): $m$ qubits, prepared in $|\psi\rangle$.

ASCII Circuit Diagram: Quantum Phase Estimation
================================================

Register 1 (t qubits)              Register 2 (m qubits)
|0⟩ ──H───●─────────────────────────────────────────────[QFT†]── Measure
            │
|0⟩ ──H───│────●────────────────────────────────────────[QFT†]── Measure
            │    │
|0⟩ ──H───│────│────●───────────────────────────────────[QFT†]── Measure
            │    │    │
  ...      │    │    │         ...                           │
|0⟩ ──H───│────│────│──────────●────────────────────────[QFT†]── Measure
            │    │    │          │
            │    │    │          │
|ψ⟩ ──────[U]──[U²]──[U⁴]──...──[U^(2^{t-1})]──────────────────── (unchanged)

Each controlled-U^(2^j) applies U^(2^j) to Register 2 when the
j-th qubit of Register 1 is |1⟩.

16.2.2 Mathematical Derivation — Step by Step

Let's trace the quantum state through each stage of the circuit with complete mathematical detail.

Step 1: Superposition. Apply $H^{\otimes t}$ to Register 1:

$$|\Psi_1\rangle = \frac{1}{\sqrt{2^t}} \sum_{x=0}^{2^t-1} |x\rangle \otimes |\psi\rangle$$

where $|x\rangle = |x_{t-1} x_{t-2} \cdots x_0\rangle$ is the binary representation of $x$.

Step 2: Controlled-$U$ operations. For each $j = 0, 1, \ldots, t-1$, apply controlled-$U^{2^j}$ with control qubit $j$:

$$|\Psi_2\rangle = \frac{1}{\sqrt{2^t}} \sum_{x=0}^{2^t-1} |x\rangle \otimes U^x |\psi\rangle$$

Since $|\psi\rangle$ is an eigenstate, $U^x |\psi\rangle = (e^{2\pi i \phi})^x |\psi\rangle = e^{2\pi i \phi x} |\psi\rangle$:

$$|\Psi_2\rangle = \frac{1}{\sqrt{2^t}} \sum_{x=0}^{2^t-1} e^{2\pi i \phi x} |x\rangle \otimes |\psi\rangle$$

This is the crucial step: the phase factors $e^{2\pi i \phi x}$ are entangled with the computational basis states $|x\rangle$.

Try It Yourself: Verify Phase Kickback

Consider $t = 3$ estimation qubits and $\phi = 0.375 = 3/8$. After Step 2:

$$|\Psi_2\rangle = \frac{1}{\sqrt{8}} \sum_{x=0}^{7} e^{2\pi i \cdot 0.375 \cdot x} |x\rangle |\psi\rangle$$

Computing the phase factors: - $x = 0$: $e^0 = 1$ - $x = 1$: $e^{2\pi i \cdot 3/8} = e^{3\pi i/4}$ - $x = 2$: $e^{2\pi i \cdot 3/4} = e^{3\pi i/2}$ - $x = 3$: $e^{2\pi i \cdot 9/8} = e^{9\pi i/4}$ - etc.

This state is the QFT of $|3\rangle$ (since $2^3 \times 0.375 = 3$), as we'll see in the next step.

Step 3: Inverse QFT. The state of Register 1 is exactly the definition of the quantum Fourier transform of the state $|2^t \phi\rangle$ (where $2^t \phi$ is interpreted as an integer if $\phi$ is a $t$-bit fraction). Applying $\text{QFT}^\dagger$:

$$|\Psi_3\rangle = \text{QFT}^\dagger_{2^t} \left( \frac{1}{\sqrt{2^t}} \sum_{x=0}^{2^t-1} e^{2\pi i \phi x} |x\rangle \right) \otimes |\psi\rangle$$

The inverse QFT acts as:

$$\text{QFT}^\dagger |k\rangle = \frac{1}{\sqrt{2^t}} \sum_{x=0}^{2^t-1} e^{-2\pi i k x / 2^t} |x\rangle$$

So applying it to the state in Register 1:

$$|\Psi_3\rangle = \frac{1}{2^t} \sum_{x=0}^{2^t-1} \sum_{k=0}^{2^t-1} e^{2\pi i x (\phi - k/2^t)} |k\rangle \otimes |\psi\rangle$$

The amplitude for $|k\rangle$ is:

$$\alpha_k = \frac{1}{2^t} \sum_{x=0}^{2^t-1} e^{2\pi i x (\phi - k/2^t)}$$

Step 4: Measurement. We measure Register 1 in the computational basis. The probability of outcome $k$ is:

$$P(k) = |\alpha_k|^2$$

Let's evaluate this sum explicitly. Using the geometric series formula:

$$\alpha_k = \frac{1}{2^t} \cdot \frac{1 - e^{2\pi i \cdot 2^t (\phi - k/2^t)}}{1 - e^{2\pi i (\phi - k/2^t)}}$$

This simplifies to:

$$\alpha_k = \frac{1}{2^t} \cdot \frac{1 - e^{2\pi i (2^t \phi - k)}}{1 - e^{2\pi i (\phi - k/2^t)}}$$

And the probability is:

$$P(k) = |\alpha_k|^2 = \frac{1}{2^{2t}} \cdot \frac{|1 - e^{2\pi i (2^t \phi - k)}|^2}{|1 - e^{2\pi i (\phi - k/2^t)}|^2}$$

Using the identity $|1 - e^{i\theta}|^2 = 4\sin^2(\theta/2)$:

$$P(k) = \frac{1}{2^{2t}} \cdot \frac{\sin^2(\pi(2^t \phi - k))}{\sin^2(\pi(\phi - k/2^t))}$$

This is sharply peaked around $k \approx 2^t \phi$.

16.2.3 Detailed Analysis of the Probability Distribution

Let $\delta = \phi - k/2^t$ be the error in the phase estimate. The probability of measuring outcome $k$ is:

$$P(k) = \frac{1}{2^{2t}} \cdot \frac{\sin^2(\pi \cdot 2^t \delta)}{\sin^2(\pi \delta)}$$

For $\delta$ close to 0, we can approximate $\sin(\pi\delta) \approx \pi\delta$:

$$P(k) \approx \frac{1}{2^{2t}} \cdot \frac{\sin^2(\pi \cdot 2^t \delta)}{(\pi \delta)^2} = \text{sinc}^2(2^t \delta)$$

This is a sinc-squared function centered at $\delta = 0$, with main lobe width $\sim 1/2^t$.

ASCII Diagram: QPE Measurement Probability Distribution
=========================================================

P(k)
  |         /\
  |        /  \
  |       /    \
  |      /      \
  |     /        \
  |    /          \
  |   /            \
  |  /              \
  | /                \
  |/__________________\_______________→ k
  |                   |
  0         k ≈ 2^t φ       2^t

  The peak is at k ≈ 2^t φ, with width ~1/2^t.
  For exact t-bit phases, P(k = 2^t φ) = 1.
  For inexact phases, the probability spreads over
  neighboring integers.

16.2.4 Precision and Success Probability

Case 1: $\phi$ is an exact $t$-bit fraction. If $\phi = k/2^t$ for some integer $k$, then $\delta = 0$ and $P(k) = 1$. The measurement always yields the exact phase.

Case 2: $\phi$ is not a $t$-bit fraction. Let $b$ be the nearest integer to $2^t \phi$, so that $2^t \phi = b + \delta'$ where $|\delta'| \leq 1/2$. Then:

$$P(b) \geq \frac{4}{\pi^2} \approx 0.405$$

This is the minimum probability of measuring the best approximation. The probability of measuring a value within $\Delta$ of the true value is:

$$P(|k - 2^t \phi| \leq \Delta) \geq 1 - \frac{1}{2(\Delta - 1)}$$

for $\Delta \geq 2$.

To achieve $n$ bits of precision with success probability $1 - \epsilon$, we need:

$$t = n + \left\lceil \log_2\left(2 + \frac{1}{2\epsilon}\right) \right\rceil$$

additional qubits. For example: - 10 bits of precision, 99% success: $t = 10 + \lceil \log_2(2 + 50) \rceil = 10 + 6 = 16$ qubits - 10 bits of precision, 99.9% success: $t = 10 + \lceil \log_2(2 + 500) \rceil = 10 + 9 = 19$ qubits - 20 bits of precision, 90% success: $t = 20 + \lceil \log_2(2 + 5) \rceil = 20 + 3 = 23$ qubits

Worked Example: QPE for $\phi = 1/3$ with $t = 4$ qubits

$\phi = 1/3 = 0.010101..._2$ (repeating). With $t = 4$ qubits, $2^t = 16$.

The nearest integers to $2^4 \phi = 16/3 \approx 5.33$ are $b = 5$ ($\phi_b = 5/16 = 0.3125$) and $b = 6$ ($\phi_b = 6/16 = 0.375$).

Probability of measuring $k = 5$: $$P(5) = \frac{1}{256} \cdot \frac{\sin^2(\pi(16/3 - 5))}{\sin^2(\pi(1/3 - 5/16))} = \frac{1}{256} \cdot \frac{\sin^2(\pi/3)}{\sin^2(\pi/48)}$$

$\sin^2(\pi/3) = 3/4$, $\sin^2(\pi/48) \approx (\pi/48)^2 \approx 0.00428$.

$P(5) \approx \frac{3/4}{256 \times 0.00428} \approx \frac{0.75}{1.096} \approx 0.685$

Similarly, $P(6)$ will also be significant. The probability is spread between $k = 5$ and $k = 6$, with no single outcome dominating.

16.2.5 What Happens with Non-Eigenstate Input

If the input state $|\psi_{\text{in}}\rangle$ is not an eigenstate of $U$ but a superposition $\sum_j c_j |\psi_j\rangle$ of eigenstates with eigenvalues $e^{2\pi i \phi_j}$, then QPE projects into one of the eigenstates and returns the corresponding phase:

$$|\psi_{\text{in}}\rangle = \sum_j c_j |\psi_j\rangle \xrightarrow{\text{QPE}} |2^t \phi_j\rangle \otimes |\psi_j\rangle$$

with probability $|c_j|^2$. This is a measurement in the eigenbasis of $U$ — a useful property for algorithms like HHL and quantum simulation.

ASCII Diagram: QPE with Superposition Input
=============================================

Input: |ψ_in⟩ = c₀|ψ₀⟩ + c₁|ψ₁⟩ + c₂|ψ₂⟩ + ...

            ┌──────┐
c₀|ψ₀⟩ ────│      │──── c₀|φ₀⟩ |ψ₀⟩
c₁|ψ₁⟩ ────│ QPE  │──── c₁|φ₁⟩ |ψ₁⟩
c₂|ψ₂⟩ ────│      │──── c₂|φ₂⟩ |ψ₂⟩
            └──────┘

Measurement outcome: |φ_j⟩ with probability |c_j|²
The eigenstate |ψ_j⟩ collapses onto its corresponding eigenvalue.

16.2.6 Detailed Derivation of the QPE Output for Exact Phases

When $\phi = m/2^t$ is an exact $t$-bit fraction, the output of QPE is deterministic. Let us trace through the full derivation.

After the controlled-$U$ operations (Step 2), the state of Register 1 is:

$$|\Psi_2\rangle = \frac{1}{\sqrt{2^t}} \sum_{x=0}^{2^t-1} e^{2\pi i \phi x} |x\rangle \otimes |\psi\rangle$$

When $\phi = m/2^t$, the phase factors become $e^{2\pi i m x / 2^t}$, and the state is exactly the QFT of $|m\rangle$:

$$|\Psi_2\rangle = \text{QFT}|m\rangle \otimes |\psi\rangle$$

Applying $\text{QFT}^\dagger$ to Register 1:

$$|\Psi_3\rangle = \text{QFT}^\dagger \cdot \text{QFT}|m\rangle \otimes |\psi\rangle = |m\rangle \otimes |\psi\rangle$$

The measurement always yields $m$, giving the exact phase $\phi = m/2^t$ with probability 1.

Numerical example: $\phi = 3/8 = 0.011_2$ with $t = 3$ estimation qubits. The controlled-$U$ operations encode:

• Qubit 0 ($j=0$): phase $2^0 \cdot 3/8 = 3/8$
• Qubit 1 ($j=1$): phase $2^1 \cdot 3/8 = 3/4$
• Qubit 2 ($j=2$): phase $2^2 \cdot 3/8 = 3/2$

After the inverse QFT, the measurement yields $|011\rangle = |3\rangle$, confirming $\phi = 3/8$.

Try It Yourself: Verify that for $\phi = 5/16$ with $t = 4$ estimation qubits, QPE yields $|0101\rangle = |5\rangle$ with certainty. Compute the phase factors for each qubit and show that they correspond to the QFT of $|5\rangle$.

Common Misconception: "QPE always gives the exact phase."

QPE gives an exact result only when $\phi$ is exactly representable as a $t$-bit binary fraction. For general $\phi$, the probability is spread over several outcomes centered around the best approximation. For example, $\phi = 1/3 = 0.010101..._2$ cannot be exactly represented with any finite number of bits, and QPE with $t = 4$ gives two major outcomes: $k = 5$ ($\phi \approx 0.3125$) and $k = 6$ ($\phi \approx 0.3750$). The more estimation qubits you use, the more precisely you can approximate $\phi$, but you never get an exact answer for irrational phases.

16.3 Qiskit Implementation: Phase Estimation for a Simple Unitary

We implement QPE for the single-qubit unitary:

$$U = R_z(2\pi \cdot 0.625) = \begin{pmatrix} e^{-i\pi \cdot 0.625} & 0 \\ 0 & e^{i\pi \cdot 0.625} \end{pmatrix}$$

The eigenstate $|1\rangle$ has eigenvalue $e^{i\pi \cdot 0.625} = e^{2\pi i \cdot 0.3125}$. We expect to estimate $\phi = 0.3125 = 5/16$.

import numpy as np
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit_aer import AerSimulator
from qiskit.circuit.library import QFT

phi_true = 0.3125          # true phase (5/16)
t = 4                      # estimation qubits (4 bits of precision)
m = 1                      # eigenstate register size

# ============================================================
# Build the QPE circuit
# ============================================================
est_reg = QuantumRegister(t, 'est')
eig_reg = QuantumRegister(m, 'eig')
meas_reg = ClassicalRegister(t, 'meas')
qc = QuantumCircuit(est_reg, eig_reg, meas_reg)

# Step 1: Superposition on estimation register
for i in range(t):
    qc.h(est_reg[i])

# Step 2: Prepare eigenstate |1⟩
qc.x(eig_reg[0])

# Step 3: Controlled-U^(2^j) operations
for j in range(t):
    power = 2 ** j
    angle = 2 * np.pi * phi_true * power
    qc.cp(angle, est_reg[j], eig_reg[0])

# Step 4: Inverse QFT on estimation register
iqft = QFT(t, inverse=True, do_swaps=True)
qc.append(iqft, est_reg[:])

# Step 5: Measure
qc.measure(est_reg, meas_reg)

# ============================================================
# Simulate
# ============================================================
simulator = AerSimulator()
from qiskit import transpile
qc_compiled = transpile(qc, simulator)
job = simulator.run(qc_compiled, shots=4096)
result = job.result()
counts = result.get_counts()

# ============================================================
# Analyze results
# ============================================================
print("Measurement results (top 8):")
for outcome, count in sorted(counts.items(),
                              key=lambda x: -x[1])[:8]:
    phase_est = int(outcome, 2) / (2**t)
    print(f"  |{outcome}⟩ → φ_est = {phase_est:.4f} "
          f"(true: {phi_true:.4f}), count = {count}")

# Find most likely phase
best_outcome = max(counts, key=counts.get)
phi_estimated = int(best_outcome, 2) / (2**t)
error = abs(phi_estimated - phi_true)
print(f"\nBest estimate: φ = {phi_estimated:.4f}")
print(f"True phase:    φ = {phi_true:.4f}")
print(f"Absolute error: {error:.6f}")
print(f"Success: {error < 1/(2**t)}")

Expected output: The measurement |0101⟩ (binary for 5) should dominate, giving $\phi_{\text{est}} = 5/16 = 0.3125$, matching the true phase exactly since it is a 4-bit fraction.

16.3.1 Non-Exact Phase Example

When $\phi$ is not an exact $t$-bit fraction, the probability distribution broadens. Let's test with $\phi = 1/3 \approx 0.3333$:

phi_true = 1/3  # not exactly representable in 4 bits
t = 4

# Rebuild circuit with new phase
qc2 = QuantumCircuit(QuantumRegister(t, 'est'),
                     QuantumRegister(m, 'eig'),
                     ClassicalRegister(t, 'meas'))
for i in range(t):
    qc2.h(i)
qc2.x(t)  # eigenstate |1⟩
for j in range(t):
    power = 2 ** j
    angle = 2 * np.pi * phi_true * power
    qc2.cp(angle, j, t)
iqft2 = QFT(t, inverse=True, do_swaps=True)
qc2.append(iqft2, range(t))
qc2.measure(range(t), range(t))

qc2_compiled = transpile(qc2, simulator)
counts2 = simulator.run(qc2_compiled, shots=4096).result().get_counts()

print("\nNon-exact phase (φ = 1/3 ≈ 0.3333):")
for outcome, count in sorted(counts2.items(),
                              key=lambda x: -x[1])[:8]:
    phase_est = int(outcome, 2) / (2**t)
    print(f"  |{outcome}⟩ → φ_est = {phase_est:.4f}, count = {count}")

# The distribution should peak around 5/16 = 0.3125 and 6/16 = 0.3750

16.3.2 Increasing Precision with More Qubits

Let's see how the precision improves with more estimation qubits:

def run_qpe(phi, t, shots=4096):
    """Run QPE with t estimation qubits and return the best phase estimate."""
    qc = QuantumCircuit(QuantumRegister(t, 'est'),
                       QuantumRegister(1, 'eig'),
                       ClassicalRegister(t, 'meas'))
    for i in range(t):
        qc.h(i)
    qc.x(t)
    for j in range(t):
        angle = 2 * np.pi * phi * (2 ** j)
        qc.cp(angle, j, t)
    iqft = QFT(t, inverse=True, do_swaps=True)
    qc.append(iqft, range(t))
    qc.measure(range(t), range(t))
    qc_compiled = transpile(qc, simulator)
    counts = simulator.run(qc_compiled, shots=shots).result().get_counts()
    best = max(counts, key=counts.get)
    return int(best, 2) / (2**t), counts[best] / shots

phi = 1/3
print("Precision vs. estimation qubits for φ = 1/3:")
print(f"{'t':>3} {'φ_est':>10} {'error':>12} {'prob':>8}")
for t in [3, 4, 5, 6, 8, 10]:
    est, prob = run_qpe(phi, t)
    error = abs(est - phi)
    print(f"{t:3d} {est:10.6f} {error:12.2e} {prob:8.4f}")

16.3.3 QPE for a Two-Qubit Unitary: The CNOT-H Gate

Let us demonstrate QPE for a two-qubit unitary $U = \text{CNOT} \cdot (H \otimes I)$. This unitary has eigenvalues that are not simple fractions, making it an instructive test case.

import numpy as np
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile
from qiskit_aer import AerSimulator
from qiskit.circuit.library import QFT
from qiskit.quantum_info import Operator, Statevector

# Build the unitary U = CNOT * (H ⊗ I)
U_circuit = QuantumCircuit(2)
U_circuit.h(0)
U_circuit.cx(0, 1)
U_matrix = Operator(U_circuit).data

# Find eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(U_matrix)
print("Eigenvalues of U:")
for i, (val, vec) in enumerate(zip(eigenvalues, eigenvectors.T)):
    phase = np.angle(val) / (2 * np.pi) % 1
    print(f"  λ_{i} = {val:.6f}, phase φ = {phase:.6f}")

# Choose an eigenstate for QPE
target_idx = 0  # use the first eigenstate
eigvec = eigenvectors[:, target_idx]
true_phase = np.angle(eigenvalues[target_idx]) / (2 * np.pi) % 1
print(f"\nTarget eigenstate index: {target_idx}")
print(f"True phase: φ = {true_phase:.6f}")

# QPE circuit
t = 5  # estimation qubits
m = 2  # eigenstate register size
est_reg = QuantumRegister(t, 'est')
eig_reg = QuantumRegister(m, 'eig')
meas_reg = ClassicalRegister(t, 'meas')
qc = QuantumCircuit(est_reg, eig_reg, meas_reg)

# Prepare eigenstate
# We need to prepare the eigenstate on eig_reg
# For a general eigenstate, we use state preparation
eigvec_normalized = eigvec / np.linalg.norm(eigvec)
# Use initialize to prepare the eigenstate
init_gate = Statevector(eigvec_normalized).to_instruction()
# We'll use statevector simulation instead
qc.eig_reg = eig_reg

# Simpler approach: use statevector simulation
from qiskit.quantum_info import Statevector

# Build QPE circuit for statevector simulation
qc_sv = QuantumCircuit(t + m)
# Prepare superposition on estimation register
for i in range(t):
    qc_sv.h(i)
# Prepare eigenstate on eigenstate register
# Use the init instruction
qc_sv.initialize(eigvec_normalized, range(t, t + m))

# Apply controlled-U^(2^j) gates
for j in range(t):
    power = 2 ** j
    # Apply U^power using matrix exponentiation
    U_power = np.linalg.matrix_power(U_matrix, power)
    # Convert to gate
    U_power_gate = Operator(U_power)
    # Apply as controlled operation
    for _ in range(power):
        qc_sv.compose(U_circuit.to_instruction(), 
                      qubits=[t, t+1], inplace=True)

# Apply inverse QFT
iqft = QFT(t, inverse=True, do_swaps=True)
qc_sv.append(iqft, range(t))

# Statevector simulation
state = Statevector(qc_sv)
# Trace out eigenstate register
probs = np.zeros(2**t)
for k in range(2**t):
    prob = 0
    for j in range(2**m):
        idx = k * 2**m + j
        prob += abs(state.data[idx])**2
    probs[k] = prob

best_k = np.argmax(probs)
est_phase = best_k / 2**t
print(f"\nQPE with t = {t} estimation qubits:")
print(f"  Best measurement: k = {best_k} = {bin(best_k)[2:].zfill(t)}")
print(f"  Estimated phase: φ = {est_phase:.6f}")
print(f"  True phase:      φ = {true_phase:.6f}")
print(f"  Error: {abs(est_phase - true_phase):.6f}")
print(f"  Probability of best outcome: {probs[best_k]:.4f}")

This example demonstrates that QPE works for multi-qubit unitaries and that the eigenstate preparation step is crucial — if we cannot prepare a good eigenstate, the phase estimate will be a superposition of different phases.

Common Misconception: "QPE requires knowing the eigenstate in advance."

For QPE to give a clean phase estimate, the input must be an eigenstate of $U$. If the input is a superposition $\sum_j c_j |\psi_j\rangle$, the measurement collapses to a random eigenstate and returns its phase. This is actually useful — it's how Shor's algorithm and quantum simulation work — but it means you get only one eigenvalue per measurement. To get a specific eigenvalue (e.g., the ground state energy), you need to prepare a state with large overlap with the desired eigenstate, which is itself a hard problem.

16.4 Applications of Quantum Phase Estimation

Before iterative phase estimation became popular, Kitaev proposed an elegant alternative that uses only single-qubit measurements but requires classical post-processing to reconstruct the phase. The idea is to measure the phase bit by bit using Hadamard tests:

For each bit $j$ (from most significant to least): 1. Prepare an ancilla qubit in $|+\rangle$ and the eigenstate register in $|\psi\rangle$. 2. Apply controlled-$U^{2^j}$ with the ancilla as control. 3. Measure the ancilla in the $X$ basis to obtain bit $j$ of the phase.

The probability of measuring $|+\rangle$ is $P(+) = (1 + \cos(2\pi \cdot 2^j \phi + \omega_j))/2$, where $\omega_j$ is a correction based on previously measured bits.

Kitaev's approach is mathematically equivalent to IPE but was developed independently and has slightly different noise properties. It's particularly useful when the controlled-$U^{2^j}$ operations are noisy, because each bit is measured independently and errors in one bit don't propagate to others.

16.5.1 Detailed Derivation of QPE Precision Bounds

Let us derive the precision bounds for QPE more carefully.

Theorem (QPE Success Probability): Let $\phi$ be the phase to be estimated, $t$ be the number of estimation qubits, and $b$ be the nearest integer to $2^t \phi$. Then:

$$P(|c - 2^t \phi| \leq \Delta) \geq 1 - \frac{1}{2(\Delta - 1)}$$

for $\Delta \geq 2$, where $c$ is the measurement outcome.

Proof sketch: The probability of measuring outcome $k$ is:

$$P(k) = \frac{1}{2^{2t}} \frac{\sin^2(\pi(2^t \phi - k))}{\sin^2(\pi(\phi - k/2^t))}$$

For $|k - 2^t \phi| > \Delta$, we can bound:

$$P(k) \leq \frac{1}{4\Delta^2}$$

by using $|\sin(\pi(\phi - k/2^t))| \geq |\phi - k/2^t| \geq \Delta/2^t$ and $|\sin(\pi(2^t \phi - k))| \leq 1$.

Summing over all $|k - 2^t \phi| > \Delta$:

$$\sum_{|k - 2^t \phi| > \Delta} P(k) \leq \sum_{|k - 2^t \phi| > \Delta} \frac{1}{4\Delta^2} \leq \frac{2^t}{4\Delta^2} \leq \frac{1}{2\Delta}$$

Wait, this bound is too loose. The tighter bound uses:

$$\sum_{|k - 2^t \phi| > \Delta} P(k) \leq \frac{1}{2(\Delta - 1)}$$

which follows from the bound $|\sin(\pi x)| \geq 2|x|$ for $|x| \leq 1/2$.

Corollary (Bits of precision): To achieve $n$ bits of precision with success probability $1 - \epsilon$:

$$t = n + \left\lceil \log_2\left(2 + \frac{1}{2\epsilon}\right) \right\rceil$$

This gives the "extra qubits" needed beyond the desired precision to boost the success probability.

Examples: - $n = 10$ bits, $\epsilon = 0.01$ (99% confidence): $t = 10 + \lceil \log_2(52) \rceil = 10 + 6 = 16$ qubits - $n = 10$ bits, $\epsilon = 0.001$ (99.9% confidence): $t = 10 + \lceil \log_2(502) \rceil = 10 + 9 = 19$ qubits - $n = 20$ bits, $\epsilon = 0.01$ (99% confidence): $t = 20 + 6 = 26$ qubits

16.5.2 Energy Resolution and Chemical Accuracy

For quantum chemistry applications, QPE is used to find eigenvalues of a molecular Hamiltonian $H$. The unitary is $U = e^{-iH\tau}$ for some evolution time $\tau$.

Energy resolution: If QPE estimates the phase $\phi$ to $t$ bits, the energy resolution is:

$$\Delta E = \frac{2\pi}{\tau \cdot 2^t}$$

For chemical accuracy ($\Delta E \leq 1.6 \times 10^{-3}$ Hartree) with $\tau = 1$ atomic unit:

$$2^t \geq \frac{2\pi}{1.6 \times 10^{-3}} \approx 3927$$

$$t \geq \lceil \log_2(3927) \rceil = 12 \text{ qubits}$$

With the overhead for success probability ($\epsilon = 0.01$), we need $t = 12 + 6 = 18$ estimation qubits.

Evolution time $\tau$: The choice of $\tau$ affects both resolution and the depth of controlled-$U^{2^j}$ operations. Larger $\tau$ gives better resolution for the same number of qubits, but requires deeper circuits.

The controlled-$U^{2^j}$ problem: For $j = t-1$ (the most significant bit), we need to implement $U^{2^{t-1}} = e^{-iH \tau 2^{t-1}}$. This requires $2^{t-1}$ repetitions of the Trotter circuit for $e^{-iH\tau}$, making the circuit exponentially deeper in $t$.

This is why QPE-based algorithms for chemistry require fault-tolerant quantum computers — the circuit depth is $O(2^t \cdot \text{Trotter\_depth})$, which for $t = 18$ and a Trotter circuit depth of 1000, gives a total depth of $2^{18} \times 1000 \approx 2.6 \times 10^8$ — far beyond NISQ capabilities.

Worked Example: QPE for H$_2$ Ground State Energy

The H$_2$ molecule in minimal basis has Hamiltonian eigenvalues (in Hartree): - Ground state: $E_0 \approx -1.85$ Hartree - First excited state: $E_1 \approx -0.57$ Hartree - Gap: $\Delta E \approx 1.28$ Hartree

To resolve the ground state energy to chemical accuracy ($1.6 \times 10^{-3}$ Hartree): - Choose $\tau = 1$ atomic unit - Need $t = 12$ estimation qubits (plus 6 for 99% confidence) - Total: 18 estimation qubits + 4 qubits for H$_2$ + ancilla $\approx$ 25 qubits - Circuit depth: $O(2^{18}) \approx 260{,}000$ Trotter steps

This is feasible on a fault-tolerant quantum computer but impossible on NISQ hardware.

16.5 Applications of Quantum Phase Estimation

QPE is arguably the most important subroutine in quantum computing. Here are its major applications:

16.5.1 Shor's Algorithm (Chapter 15)

In Shor's algorithm, $U_a |y\rangle = |ay \bmod N\rangle$. The eigenstates of $U_a$ are:

$$|\psi_k\rangle = \frac{1}{\sqrt{r}} \sum_{j=0}^{r-1} e^{-2\pi i k j / r} |a^j \bmod N\rangle$$

with eigenvalues $e^{2\pi i k / r}$. QPE extracts $k/r$, and continued fractions recover $r$.

The connection is precise: Shor's period-finding algorithm IS quantum phase estimation applied to the modular multiplication unitary. The only difference is that we don't know the eigenstates of $U_a$, so we use the trick of preparing a superposition of eigenstates by applying a Hadamard on the control register and modular exponentiation on the target register. The QPE circuit naturally "kicks back" the phase information into the control register.

16.5.2 The HHL Algorithm for Linear Systems

The Harrow-Hassidim-Lloyd (HHL) algorithm solves $A\mathbf{x} = \mathbf{b}$ for sparse, well-conditioned $A$. It uses QPE to decompose $|\mathbf{b}\rangle$ in the eigenbasis of $A$:

1. QPE on $e^{iAt}$: Decompose $|\mathbf{b}\rangle = \sum_j \beta_j |\lambda_j\rangle$ into eigenstates of $A$ with eigenvalues $\lambda_j$.
2. Eigenvalue inversion: Apply a controlled-rotation that maps $|\lambda_j\rangle$ to $C/\lambda_j |\lambda_j\rangle$.
3. Inverse QPE: Uncompute the eigenvalue register.
4. Result: The state $\sum_j (C/\lambda_j) \beta_j |\lambda_j\rangle \propto |A^{-1}\mathbf{b}\rangle = |\mathbf{x}\rangle$.

The runtime is $O(\log N \cdot s^2 \kappa^2 / \epsilon)$, exponentially faster than classical $O(N \sqrt{\kappa})$ for dense matrices, where $N$ is the matrix dimension, $s$ is the sparsity, $\kappa$ is the condition number, and $\epsilon$ is the precision.

Important caveat: HHL outputs a quantum state $|\mathbf{x}\rangle$, not a classical vector. Reading out all $N$ components would require $O(N)$ measurements, destroying the quantum speedup. HHL is useful when we only need specific properties of the solution (e.g., $\langle \mathbf{x} | M | \mathbf{x} \rangle$ for some observable $M$), not the full vector.

16.5.3 Quantum Simulation (Chapter 17)

Given a Hamiltonian $H$, the time evolution operator is $U(t) = e^{-iHt}$. QPE applied to $U(t)$ extracts the eigenvalues of $H$ — i.e., the energy spectrum. This is the foundation of quantum chemistry on quantum computers.

For a molecular Hamiltonian with eigenstates $|E_j\rangle$ and energies $E_j$:

$$U(t) |E_j\rangle = e^{-iE_j t} |E_j\rangle$$

QPE with $U(t)$ and input state $|E_j\rangle$ yields phase $\phi_j = -E_j t / (2\pi)$, from which we recover $E_j = -2\pi \phi_j / t$.

To resolve energy differences $\Delta E$, we need:

$$t \geq \log_2\left( \frac{2\pi}{\tau \cdot \Delta E} \right)$$

estimation qubits. For chemical accuracy ($\sim 1.6 \times 10^{-3}$ Hartree $\approx 0.04$ eV), with $\tau = 1$ (atomic units), this requires $t \approx 12$–$16$ qubits.

16.5.4 Quantum Counting and Amplitude Estimation

QPE can estimate the amplitude of a marked state in Grover's algorithm. If there are $M$ marked items in a database of $N$ items, Grover's search applies the unitary:

$$Q = -A S_0 A^{-1} S_f$$

where $A$ is the diffusion operator and $S_f$ marks the solutions. The eigenvalues of $Q$ are $e^{\pm 2i\theta}$ where $\sin^2\theta = M/N$. QPE on $Q$ yields $\theta$, from which we extract $M = N\sin^2\theta$.

This gives a quadratic speedup over classical counting: $O(\sqrt{N/M})$ quantum queries vs. $O(N)$ classical queries.

16.5.5 Summary of QPE Applications

16.6 Connection to Hamiltonian Simulation

The most important application of QPE is extracting eigenvalues of a Hamiltonian $H$. If we can implement $U = e^{-iH\tau}$ for some small time step $\tau$, then:

$$U |E_j\rangle = e^{-iE_j \tau} |E_j\rangle$$

QPE on $U$ with eigenstate $|E_j\rangle$ yields the phase $\phi_j = -E_j \tau / (2\pi)$ (mod $2\pi$), from which we recover the energy $E_j$.

Choosing $\tau$: The phase must satisfy $0 \leq \phi < 1$, so we need $0 \leq E_j \tau / (2\pi) < 1$ for all eigenvalues. This requires $\tau < 2\pi / E_{\max}$. In practice, we shift the Hamiltonian so that all eigenvalues are positive ($H \to H - E_{\min} I$) and choose $\tau$ appropriately.

Implementing controlled-$U^{2^j}$: For a local Hamiltonian $H = \sum_k H_k$, we use Trotterization (Chapter 17) to approximate $e^{-iH\tau}$ as a product of $e^{-iH_k\tau}$ terms. Then $e^{-iH \cdot 2^j \tau}$ is implemented by running the Trotter circuit with time step $2^j \tau$, or equivalently, by repeating the Trotter circuit $2^j$ times.

The gate cost of QPE for chemistry scales as:

$$\text{Gate count} \sim O\left( \frac{t \cdot \eta^4}{\epsilon} \right)$$

where $\eta$ is the number of spin-orbitals and $\epsilon$ is the target precision. For a molecule with 100 orbitals, 16 bits of precision, and chemical accuracy, this amounts to roughly $10^{10}$ gates — far beyond current NISQ capabilities but feasible for fault-tolerant quantum computers.

Recurring Theme: Noise is the Enemy

QPE requires deep circuits — controlled-$U^{2^j}$ for $j$ up to $t-1$ means applying $U$ up to $2^{t-1}$ times. For $t = 16$, that's over 32,000 applications of $U$. Each application introduces errors that accumulate. On a NISQ device with 99.5% gate fidelity, 32,000 gates would produce essentially random noise. This is why QPE-based algorithms (Shor's, HHL, simulation) require fault-tolerant quantum computers, and why NISQ-era algorithms (VQE, QAOA) avoid deep circuits.

16.7 Iterative Phase Estimation

16.7.1 The IPE Algorithm

An alternative to the multi-qubit QFT-based approach is iterative phase estimation (IPE), which uses a single ancilla qubit and classical feedback:

ASCII Circuit Diagram: Iterative Phase Estimation (One Bit per Round)
====================================================================

Round k (determining bit φ_{t-k}):

|0⟩ ──H───R_z(ω_k)───●───R_z(-π/2)───H─── Measure → bit φ_{t-k}
                       │
|ψ⟩ ─────────────────[U^(2^{t-k-1})]─────────────────────

where ω_k = -2π Σ_{j=0}^{k-1} φ_{t-j} / 2^{k-j}

After measuring bit k, the classical feedback angle ω_k is updated
based on previously measured bits. This eliminates the need for
the multi-qubit QFT and reduces qubit count at the cost of
sequential measurements.

The IPE algorithm works as follows:

1. Start with $k = t-1$ (most significant bit).
2. In round $k$, prepare the ancilla in $|+\rangle$, apply controlled-$U^{2^k}$ with a phase correction $\omega_k$ based on previously measured bits, measure the ancilla in the X basis, and record the outcome as bit $\phi_{t-1-k}$.
3. Repeat for $k = t-2, t-3, \ldots, 0$.
4. The measured bits $\phi_{t-1}, \phi_{t-2}, \ldots, \phi_0$ give the binary expansion of $\phi = 0.\phi_{t-1}\phi_{t-2}\cdots\phi_0$.

Advantages of IPE over standard QPE: - Qubit count: IPE needs $m + 1$ qubits (one ancilla plus the eigenstate register), compared to $t + m$ for standard QPE. - Circuit depth: Each round is shallow (one controlled-$U^{2^k}$), but $t$ rounds are needed sequentially. - Robustness: Since each round uses a shallow circuit, IPE is more resilient to decoherence on NISQ devices.

16.7.2 Qiskit Implementation of IPE

import numpy as np
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit_aer import AerSimulator
from qiskit import transpile

def iterative_phase_estimation(phi_true, t, num_shots=1024):
    """
    Iterative phase estimation for U = R_z(2πφ) on eigenstate |1⟩.

    Parameters:
        phi_true: true phase φ ∈ [0, 1)
        t: number of bits of precision
        num_shots: shots per measurement

    Returns:
        estimated phase φ_est
    """
    measured_bits = []

    for k in range(t - 1, -1, -1):
        # Phase correction from previously measured bits
        omega = 0.0
        for j, bit in enumerate(measured_bits):
            omega += bit / (2 ** (j + 1))
        omega *= 2 * np.pi

        # Build single-round circuit
        anc = QuantumRegister(1, 'anc')
        eig = QuantumRegister(1, 'eig')
        meas = ClassicalRegister(1, 'meas')
        qc = QuantumCircuit(anc, eig, meas)

        # Prepare eigenstate |1⟩
        qc.x(eig[0])

        # Apply H on ancilla
        qc.h(anc[0])

        # Apply phase correction
        if abs(omega) > 1e-10:
            qc.rz(omega, anc[0])

        # Controlled-U^(2^k)
        angle = 2 * np.pi * phi_true * (2 ** k)
        qc.cp(angle, anc[0], eig[0])

        # Final rotation and H
        qc.rz(-np.pi / 2, anc[0])
        qc.h(anc[0])

        # Measure
        qc.measure(anc[0], meas[0])

        # Run
        simulator = AerSimulator()
        qc_compiled = transpile(qc, simulator)
        result = simulator.run(qc_compiled, shots=num_shots).result()
        counts = result.get_counts()

        # Determine bit (majority vote)
        bit = 1 if counts.get('1', 0) > counts.get('0', 0) else 0
        measured_bits.append(bit)

    # Reconstruct phase from bits
    phi_est = 0.0
    for j, bit in enumerate(measured_bits):
        phi_est += bit / (2 ** (j + 1))

    return phi_est

# Test IPE
phi_true = 0.3125  # 5/16
for t in [3, 4, 5, 6, 8]:
    phi_est = iterative_phase_estimation(phi_true, t)
    error = abs(phi_est - phi_true)
    print(f"t = {t}: φ_est = {phi_est:.6f}, error = {error:.2e}")

Comparison of QPE vs. IPE:

Recurring Theme: We're at the Beginning

IPE exemplifies the creative algorithmic workarounds that the NISQ era demands. By trading parallel qubits for sequential measurements, IPE reduces the quantum resource requirements at the cost of more classical processing. This is a pattern we'll see repeatedly: hybrid classical-quantum algorithms that squeeze useful computation out of imperfect hardware. As hardware improves, we'll transition back to the more elegant but resource-intensive standard QPE.

16.8 Error Analysis and Mitigation

16.8.1 Sources of Error in QPE

1. Finite precision: Using $t$ estimation qubits limits resolution to $2^{-t}$. This is a fundamental limitation — to get more bits of the phase, you need more qubits.

2. Non-eigenstate input: If the input state $|\psi_{\text{in}}\rangle = \sum_j c_j |\psi_j\rangle$ is a superposition of eigenstates, the measurement yields $\phi_j$ with probability $|c_j|^2$. This is actually useful — it's how Shor's algorithm and quantum simulation work — but it means a single measurement only gives one eigenvalue.

3. Gate errors: Imperfect controlled-$U^{2^j}$ operations introduce coherent and incoherent errors. The controlled-$U^{2^j}$ for large $j$ requires many repetitions of $U$, amplifying gate errors.

4. Decoherence: Loss of phase coherence during the computation smears the QFT peaks. The estimation register must maintain coherence for the entire circuit depth, which scales as $O(t^2 + t \cdot \text{depth}(U^{2^{t-1}}))$.

16.8.2 Error Mitigation Strategies

Readout error mitigation: The measurement error matrix $M_{ij} = P(\text{measure } i \mid \text{prepared } j)$ can be characterized by preparing all $2^t$ basis states and measuring. Inverting $M$ (or using pseudoinverse) corrects the measurement statistics:

$$\vec{p}_{\text{corrected}} = M^{-1} \vec{p}_{\text{measured}}$$

Zero-noise extrapolation (ZNE): Run QPE at multiple noise levels (by stretching gate durations or inserting identity-equivalent gate sequences), then extrapolate to the zero-noise limit. The key insight is that the phase estimate $\phi(\epsilon)$ is a smooth function of the noise rate $\epsilon$.

Symmetry verification: If the unitary $U$ has a symmetry (e.g., particle number conservation for molecular Hamiltonians), check whether the output state satisfies the symmetry. Discard runs that violate the symmetry — this is a form of post-selection that can dramatically improve results.

Randomized compiling: Replace deterministic gate sequences with randomly chosen but logically equivalent alternatives. This converts coherent errors (which accumulate) into stochastic errors (which average out), making ZNE more effective.

16.9 The Kitaev Phase Estimation Algorithm

An intermediate approach between standard QPE and IPE is the Kitaev algorithm (also called semi-classical phase estimation), which uses multiple rounds with different numbers of repetitions:

1. For bit $j$ (from most significant to least), prepare a single ancilla qubit.
2. Apply controlled-$U^{2^j}$.
3. Measure the ancilla with a phase correction based on previously measured bits.
4. Repeat steps 1-3 for each bit.

The key difference from IPE is that Kitaev's original formulation uses the semiclassical approximation — the measurement result is used classically, and the ancilla qubit is reinitialized each round. This is mathematically equivalent to IPE but was discovered independently.

The advantage over IPE is that each round can be repeated many times to improve the confidence in each bit, at the cost of more total circuit executions.

16.10 Advanced Topic: Robust Phase Estimation

Standard QPE assumes access to perfect controlled-$U^{2^j}$ operations. In practice, these operations are noisy. Robust phase estimation protocols modify QPE to handle noise:

Bayesian phase estimation: Instead of fixed phases $2^j$, choose phases $\phi_k$ adaptively based on previous measurements. After each measurement, update a Bayesian posterior over $\phi$:

$$P(\phi \mid \text{data}) \propto P(\text{data} \mid \phi) \cdot P_{\text{prior}}(\phi)$$

The next phase $\phi_{k+1}$ is chosen to maximize the expected information gain. This approach achieves the Heisenberg limit ($1/n$ scaling of uncertainty with number of measurements) and is more robust to noise than standard QPE.

Bayesian QPE algorithm:

def bayesian_phase_estimation(U_implementation, prior_samples=1000, 
                                max_iterations=20, confidence=0.99):
    """
    Bayesian phase estimation: adaptively choose measurement phases
    to maximize information gain about the unknown phase φ.
    """
    from scipy.stats import norm

    # Initialize prior: uniform on [0, 1)
    phi_samples = np.random.uniform(0, 1, prior_samples)
    weights = np.ones(prior_samples) / prior_samples

    for iteration in range(max_iterations):
        # Choose next evolution time to maximize expected information
        # Heuristic: choose t_k such that the prior variance is maximally reduced
        t_k = 2 ** iteration  # or use more sophisticated adaptive scheme

        # Run the circuit and get measurement outcome
        outcome = run_qpe_round(U_implementation, t_k)

        # Bayesian update
        for i in range(prior_samples):
            # P(measurement | φ) for each φ sample
            prob = 0.5 * (1 + np.cos(2 * np.pi * t_k * phi_samples[i] - outcome * np.pi))
            weights[i] *= prob

        # Normalize
        weights /= weights.sum()

        # Check convergence
        phi_est = np.average(phi_samples, weights=weights)
        variance = np.average((phi_samples - phi_est) ** 2, weights=weights)

        if variance < (1 - confidence) ** 2:
            break

    return np.average(phi_samples, weights=weights)

16.11 The QFT Circuit in Detail

16.11.1 Decomposing the QFT into Elementary Gates

The quantum Fourier transform on $n$ qubits, $\text{QFT}_{2^n}$, can be decomposed into $O(n^2)$ elementary gates. The circuit consists of:

1. Hadamard on qubit $j$: $H_j$
2. Controlled rotations: For each pair $(j, k)$ with $j < k$, apply $R_k = \text{CP}(2\pi/2^{k-j+1})$ with qubit $k$ as control and qubit $j$ as target.
3. SWAP gates: Reverse the qubit order.

The circuit for $\text{QFT}_8$ (3 qubits) is:

ASCII Circuit: QFT on 3 Qubits
===============================

|j₂⟩ ──H──R₂──R₃──×──
          │    │   │
|j₁⟩ ─────H──R₂──×─×──
               │   │ │
|j₀⟩ ──────────H──×─×──

R_k = controlled phase rotation by angle 2π/2^k
× = SWAP gates (to reverse qubit order)

Gate count analysis: - Hadamard gates: $n$ - Controlled rotations: $n(n-1)/2$ - SWAP gates: $n/2$ (each SWAP requires 3 CNOTs) - Total: $n + n(n-1)/2 + 3n/2 = O(n^2)$

For Shor's algorithm factoring an $n$-bit number, the QFT requires $O(n^2) = O((\log N)^2)$ gates — exponentially fewer than the classical FFT's $O(N \log N) = O(2^n n)$.

16.11.2 The Inverse QFT

The inverse QFT ($\text{QFT}^\dagger$) is simply the reverse of the QFT circuit with conjugated rotations:

$$\text{QFT}^\dagger |k\rangle = \frac{1}{\sqrt{2^n}} \sum_{j=0}^{2^n-1} e^{-2\pi i jk/2^n} |j\rangle$$

In practice, $\text{QFT}^\dagger$ is implemented by running the QFT circuit backward with negated rotation angles. The gate count is identical to the QFT.

Try It Yourself: Implement the QFT from Scratch

Instead of using Qiskit's built-in QFT, try implementing it manually:

python def qft_circuit(n): """Build QFT on n qubits from elementary gates.""" qc = QuantumCircuit(n) for j in range(n): qc.h(j) for k in range(j+1, n): angle = 2 * np.pi / (2 ** (k - j + 1)) qc.cp(angle, k, j) # SWAP to reverse qubit order for i in range(n // 2): qc.swap(i, n - 1 - i) return qc

Verify that your implementation produces the same results as Qiskit's built-in QFT on a small example (e.g., 3 qubits).

16.11.3 Approximate QFT

The controlled-phase rotations $R_k$ for large $k$ have very small angles. Since these rotations become exponentially small, they can be omitted with only a small reduction in fidelity. The approximate QFT discards all rotations $R_k$ with $k > m$ for some cutoff $m$, reducing the gate count from $O(n^2)$ to $O(nm)$.

Error bound: If we drop all rotations with $|angle| < 2\pi/2^{m+1}$, the error in the QFT output is:

$$\|\text{QFT} - \text{QFT}_{\text{approx}}\| \leq \frac{n}{2^{m+1}}$$

For QPE with $t$ estimation qubits, we need $m \geq t + \log_2(n) + \text{constant}$ to ensure the approximation error doesn't affect the phase estimate.

For Shor's algorithm factoring an $n$-bit number, the approximate QFT with $m = O(\log n)$ controlled rotations per qubit gives a total gate count of $O(n \log n)$ for the QFT — a significant saving over the exact $O(n^2)$.

16.12 Phase Estimation in the Presence of Noise

16.12.1 Depolarizing Noise

Depolarizing noise replaces the quantum state $\rho$ with a mixture of the original state and the maximally mixed state:

$$\mathcal{E}_p(\rho) = (1-p)\rho + \frac{p}{2^n} I$$

For QPE with $t$ estimation qubits and depolarizing rate $p$ per gate, the effective noise on the phase estimate is:

$$\text{Phase error} \approx O(p \cdot t \cdot \text{depth}(U))$$

For a circuit with $D$ total gates and depolarizing rate $p$, the output distribution is a mixture:

$$P_{\text{noisy}}(k) = (1-p)^D P_{\text{ideal}}(k) + \left(1 - (1-p)^D\right) P_{\text{uniform}}(k)$$

where $P_{\text{uniform}}(k) = 1/2^t$ is the uniform distribution. The signal-to-noise ratio is:

$$\text{SNR} = \frac{(1-p)^D \max_k P_{\text{ideal}}(k)}{(1-(1-p)^D)/2^t} \approx \frac{(1-p)^D \cdot 4/\pi^2}{(1-(1-p)^D)/2^t}$$

For QPE to work, we need $\text{SNR} \gg 1$, which requires $p \ll 1/D$. This is the fundamental noise threshold for QPE.

16.12.2 Coherent Errors

Coherent errors (systematic over-rotations or under-rotations) are more pernicious than stochastic errors because they accumulate constructively rather than averaging out. A systematic over-rotation of $\epsilon$ per gate accumulates as $O(\epsilon D)$, while a stochastic error of $\epsilon$ per gate accumulates as $O(\epsilon\sqrt{D})$.

For QPE, coherent errors in the controlled-$U^{2^j}$ gates can shift the phase estimate by:

$$\Delta\phi \approx \sum_{j=0}^{t-1} 2^j \epsilon_j$$

where $\epsilon_j$ is the error in the $j$-th controlled rotation. The $2^j$ factor means that errors in the most significant bits (largest $j$) are amplified exponentially.

Mitigation: Randomized compiling converts coherent errors into stochastic errors by randomly selecting equivalent gate sequences. This replaces the $O(\epsilon D)$ coherent accumulation with $O(\epsilon\sqrt{D})$ stochastic accumulation, significantly improving QPE performance.

16.13 QPE and the Quantum Linear Systems Algorithm (HHL)

The HHL algorithm (Harrow, Hassidim, Lloyd, 2009) solves linear systems $A\mathbf{x} = \mathbf{b}$ by using QPE as a subroutine. Here's how it works:

1. Encode the vector $\mathbf{b}$ as a quantum state $|b\rangle = \sum_i b_i |i\rangle$.
2. Apply QPE to $U = e^{iAt}$ (or $U = e^{-iAt}$) to decompose $|b\rangle$ in the eigenbasis of $A$:

$$|b\rangle = \sum_j \beta_j |u_j\rangle \xrightarrow{\text{QPE}} \sum_j \beta_j |u_j\rangle |\lambda_j\rangle$$

3. Invert the eigenvalues using a controlled rotation: $|\lambda_j\rangle \to C/\lambda_j |\lambda_j\rangle$, where $C$ is a normalization constant.
4. Apply inverse QPE to uncompute the eigenvalue register.
5. Result: The state $\sum_j (C/\lambda_j) \beta_j |u_j\rangle \propto |x\rangle = A^{-1}|b\rangle$.

Key properties of HHL: - Runtime: $O(\log N \cdot s^2 \kappa^2 / \epsilon)$ where $N$ is the matrix dimension, $s$ is the sparsity, $\kappa$ is the condition number, and $\epsilon$ is the desired precision. - Exponential speedup over classical $O(N\sqrt{\kappa})$ for dense matrices. - But: only produces a quantum state $|x\rangle$, not a classical vector. Reading out all components would require $O(N)$ measurements. - Useful when we only need specific properties of $x$ (e.g., $\|x\|$, $x^\top M x$), not the full vector.

Common Misconception: "HHL solves linear systems exponentially faster than classical algorithms."

HHL produces a quantum state encoding the solution, not a classical vector. To extract a single component of $x$, you need $O(1)$ measurements; to extract all components, you need $O(N)$ measurements, destroying the quantum speedup. HHL is useful only when we need global properties of the solution (e.g., its norm, or an inner product with another vector), not the full solution vector. This is analogous to how QPE produces a phase estimate, not the full eigenvector.

16.14 QPE for Quantum Chemistry: A Complete Example

16.14.1 The Chemistry QPE Protocol

Quantum phase estimation is the gold-standard algorithm for computing molecular energies on a fault-tolerant quantum computer. Here's the complete protocol:

1. Hamiltonian preparation: Compute the molecular Hamiltonian $H$ in second quantization using PySCF or a similar package.
2. Fermion-to-qubit mapping: Transform $H$ to a qubit Hamiltonian using Jordan-Wigner or Bravyi-Kitaev.
3. Initial state preparation: Prepare a Hartree-Fock state or an ADAPT-VQE state with good overlap with the ground state.
4. Controlled time evolution: Implement controlled-$e^{-iHt}$ using Trotterization or qubitization.
5. QPE: Apply the QPE circuit to estimate the ground state energy.
6. Classical post-processing: Convert the measured phase to an energy using $E_0 = -2\pi\phi/t + E_{\text{shift}}$.

16.14.2 Resource Estimates for Chemistry QPE

The cost of QPE for chemistry depends on several factors:

Estimation qubits: To resolve the ground state energy to chemical accuracy ($\epsilon = 1.6$ mHartree) with $t$ estimation qubits:

$$t \geq \log_2\left(\frac{2\pi}{\tau \epsilon}\right) + \log_2\left(\frac{2 + 1/(2\epsilon_{\text{prob}})}{1}\right)$$

For $\tau = 1$ and $\epsilon = 1.6 \times 10^{-3}$ Hartree, this gives $t \approx 12 + 6 = 18$ estimation qubits.

Trotter steps per controlled-$U$: For second-order Trotter with $L$ Hamiltonian terms and time step $\Delta t$:

$$\text{Trotter steps per } U^{2^j} = O\left(\frac{2^j \tau \lambda^{3/2}}{\epsilon^{1/2}}\right)$$

where $\lambda = \sum_j \|H_j\|$ is the Hamiltonian norm.

Total gate count: For the FeMoco active space (113 orbitals):

This is within reach of a fault-tolerant quantum computer but far beyond NISQ capabilities.

16.14.3 Comparison: QPE vs. VQE for Chemistry

Recurring Theme: We're at the Beginning

QPE is the "gold standard" algorithm for quantum chemistry, but it requires fault-tolerant quantum computers that don't yet exist. VQE is the "bronze standard" — it works on NISQ devices but doesn't have proven speedups. As hardware improves, we'll transition from VQE to QPE, but this transition is likely 10-15 years away for molecules beyond the reach of classical methods.

16.14.4 QPE Simulation for H$_2$ (Educational)

For small systems like H$_2$, we can simulate QPE classically to verify the algorithm:

import numpy as np
from scipy.linalg import expm

# H2 Hamiltonian in 2-qubit reduced form
# H = f0*I + f1*Z0 + f2*Z1 + f3*Z0*Z1 + f4*X0*X1 + f5*Y0*Y1
# (Coefficients from Qiskit Nature for bond length 0.735 Å)
f0, f1, f2, f3, f4, f5 = -1.0524, 0.3979, -0.3979, -0.0113, 0.1809, 0.1809

# Build the 4x4 matrix
I = np.eye(2)
Z = np.array([[1, 0], [0, -1]])
X = np.array([[0, 1], [1, 0]])
Y = np.array([[0, -1j], [1j, 0]])

H = (f0 * np.kron(I, I) + f1 * np.kron(Z, I) + f2 * np.kron(I, Z) +
     f3 * np.kron(Z, Z) + f4 * np.kron(X, X) + f5 * np.kron(Y, Y))

# Exact eigenvalues
eigenvalues = np.linalg.eigvalsh(H)
print("Exact eigenvalues:", eigenvalues)
print("Ground state energy:", eigenvalues[0])

# Simulate QPE with t estimation qubits
t = 10  # estimation qubits
tau = 1.0  # evolution time
U = expm(-1j * H * tau)

# Ground state (eigenstate with lowest energy)
_, eigvecs = np.linalg.eigh(H)
psi_0 = eigvecs[:, 0]  # ground state

# Compute phase for ground state
phase_0 = eigenvalues[0] * tau / (2 * np.pi)  # should be negative
print(f"\nGround state phase: {phase_0:.6f}")
print(f"Ground state energy: {eigenvalues[0]:.6f} Hartree")
print(f"Ground state energy: {eigenvalues[0] * 27.2114:.6f} eV")

# Simulate measurement probabilities
N = 2**t
probs = np.zeros(N)
for k in range(N):
    # Probability of measuring k given phase phi
    phi = -phase_0  # QPE measures positive phase
    alpha_k = (1/N) * np.sum([np.exp(2j * np.pi * (phi - k/N) * x) for x in range(N)])
    # More precise formula
    delta = phi - k/N
    if abs(np.sin(np.pi * delta)) < 1e-10:
        probs[k] = 1.0 if abs(delta) < 1/(2*N) else 0.0
    else:
        probs[k] = abs(np.sin(np.pi * N * delta) / (N * np.sin(np.pi * delta)))**2

# Normalize
probs /= probs.sum()

best_k = np.argmax(probs)
est_phase = best_k / N
est_energy = -2 * np.pi * est_phase / tau
print(f"\nQPE with {t} estimation qubits:")
print(f"  Best measurement: k = {best_k}")
print(f"  Estimated phase: {est_phase:.6f}")
print(f"  Estimated energy: {est_energy:.6f} Hartree")
print(f"  Exact energy: {eigenvalues[0]:.6f} Hartree")
print(f"  Error: {abs(est_energy - eigenvalues[0]):.2e} Hartree")
print(f"  Chemical accuracy achieved: {abs(est_energy - eigenvalues[0]) < 1.6e-3}")